dragonfly_sim.core.reduced_order_model

POD reduced-order model (ROM) of a steady flow solve: least-squares Petrov-Galerkin (LSPG) on a POD basis of converged solutions.

The ROM works with any model built on CompressibleEulerModel – the Euler model itself and CompressibleRANSModel, which forwards to it – through this interface, and nothing model-specific:

u_vec                     the state (its layout is the basis's row layout)
set_up_solver(), problem  the Newton problem: problem.F(x, b) and
                          problem.J(x, problem._A) assemble the residual
                          and its exact Jacobian at the current u_vec
                          (for RANS the condensed-gradient problem, which
                          updates the reconstructed gradient on the way)
steady_residual()         context in which problem assembles the steady
                          residual (pseudo/physical time terms off)
positivity_step_length    step length keeping density and pressure positive
state_scales()            per-component magnitudes for the inner product
_assemble_physical_residual(), relative_residual_norm()

Inner product. Snapshots and basis use (u, v)_W = u^T W v with, by default, the state-scaled Euclidean weight W = D^-2, D the per-component state scales s_k: every component on the same footing (in absolute terms the SA variable rho*nu_tilde is ~1e-4 against rho*E ~ 2, and a truncated SVD would ignore it). inner_product=’l2’ adds the (diagonal) DG0 mass matrix, W = M D^-2, i.e. int sum_k u_k v_k / s_k^2 dx; on the NACA0012 alpha sweeps it was 10-100x less accurate than the Euclidean weight (M emphasizes the large far-field cells). W is built once on the undeformed mesh: the incremental SVD keeps Q^T W Q = I only for a fixed W. The LSPG residual is measured in the matching norm ||S r||, S = D^-1 (D^-1 M^-1/2 with ‘l2’, on the current mesh). state_scaling=None with ‘euclidean’ is the plain Euclidean POD/LSPG.

Solve. With V the basis (V^T W V = I), u = V a (+ mean), LSPG is the Gauss-Newton iteration on min ||S R(V a)||:

(S A V)^T (S A V) da = -(S A V)^T S r, a <- a + theta da,

theta from positivity_step_length, optionally backtracked (Armijo) on ||S r||^2 / 2. The reduced system is small and solved redundantly on every rank. Convergence follows project-sandbox benchmark_pod.py (its PODSolver and outer loop): passes of at most newton_max_it iterations, each stopping on ||g|| = ||(S A V)^T S r|| < max(atol, rtol ||g_0||), with g_0 from the first iteration of the first pass (not reset by a restart), on a negligible step, or on a stall of ||g||; up to max_outer_iterations passes, converged when a pass ends with ||g|| < max(newton_atol, g_conv_limit ||g_0||), g_0 the same reference.

Derivatives. The ROM solves the LSPG stationarity condition, the projected residual equations

Phi^T R(V a; x) = 0, Phi = S^T S A V,

so when its solution is accepted the adjoint is that of these k equations (ReducedOrderModel.adjoint, an LSPGAdjoint). With V, S and Phi frozen at the solution (the basis is not differentiated, and the Gauss-Newton terms, which are proportional to S R, are dropped), for an output f with g = (df/du)^T:

(S A V)^T (S A V) lam = V^T g, psi = Phi lam = S^T S A V lam,

and df/dx = df/dx|_u - psi^T dR/dx with the full-order dR/dx: psi takes the place of the full-order adjoint A^-T g. It is the Petrov-Galerkin solution of A^T psi = g in the span of Phi, tested with V: V^T (A^T psi - g) = 0. psi is computed as S Q y with S A V = Q R (a distributed thin QR) and R^T y = V^T g, without forming the normal matrix, whose condition number is squared.

Classes

`LSPGAdjoint`	The adjoint of an accepted LSPG solution (Derivatives in the module
`LSPGSolver`	Reduced Newton (Gauss-Newton) loop on a POD basis, driving a model's
`ROMResult`
`ReducedOrderModel`	POD snapshots, basis and LSPG solve for DG_windtunnel_model

Functions

scaled_jacobian_basis(model, V, S[, AV])

S A V at the model's current u_vec, A the exact Jacobian assembled by its

Module Contents

class dragonfly_sim.core.reduced_order_model.LSPGAdjoint(V, S, comm, parameter_vector, solution)

The adjoint of an accepted LSPG solution (Derivatives in the module docstring): psi = S Q y, with S A V = Q R and R^T y = V^T g, for the basis V and residual weight S of that solve. The factorization is formed by the first solve() and reused for every later right-hand side (CSDL passes one per output row). It belongs to one design and state; matches() checks both.

Parameters:

VPETSc.Mat: The basis of the solve (owned by the POD, not by this object).
SPETSc.Vec or None: The diagonal residual weight of the solve; destroyed with this object.
parameter_vector, solution: The design and this rank’s owned entries of the accepted solution.

_factor(model): Q, R of S A V, A the steady Jacobian at the model’s current u_vec.

_galerkin_residual(model, psi, rhs, b): |V^T (A^T psi - g)| / |V^T g|, with the Jacobian _factor assembled (zero up to round-off).

destroy()

matches(parameter_vector, u_local): Whether this is the adjoint at design parameter_vector and state u_local (owned entries): both equal to the accepted solve’s on every rank. Collective.

solve(model, rhs_local): psi (this rank’s owned entries) for the right-hand side g, rhs_local being its owned entries. The model’s u_vec, mesh and angle of attack must be those of the accepted solution (the first call assembles the Jacobian there). Collective.

S

V

_Q = None

_R = None

comm

parameter_vector

property rb_size

solution

class dragonfly_sim.core.reduced_order_model.LSPGSolver(V, comm, W=None, S=None, mean=None, rtol=1e-10, atol=1e-12, max_it=30, stall_window=3, stall_rtol=1e-10, line_search=True, armijo_c=0.0001, max_backtracks=10)

Reduced Newton (Gauss-Newton) loop on a POD basis, driving a model’s Newton problem (see the module docstring).

Parameters:

VPETSc.Mat: The basis, a distributed MATDENSE in the state layout.
WPETSc.Vec or None: The diagonal inner-product weight, for the restriction a = V^T W (u - mean); None: Euclidean (V^T V = I).
SPETSc.Vec or None: The diagonal residual weight; None: the plain residual norm.
meanPETSc.Vec or None: Optional affine offset (u = V a + mean).

_backtrack(model, a, da, theta, phi, slope): Largest theta * 2^-k (k <= max_backtracks) with phi(a + theta da) <= phi + armijo_c theta slope; the given theta if there is none.

_reduced_system(model): (H, g, ||S r||) with H = (S A V)^T (S A V), g = (S A V)^T S r at the current u_vec.

_weighted_residual(model): Assemble r at the current u_vec; returns (S r local entries, ||S r||).

destroy()

expand(a, u): u <- V a (+ mean), ghosts updated.

restrict(u_vec): a = V^T W (u - mean), replicated on every rank.

solve(model)

One pass of the LSPG Gauss-Newton loop: advance model.u_vec, starting from its restriction onto the basis, towards the LSPG solution in the span of V. Call inside model.steady_residual() for a steady solve. Returns (n_iterations, converged). ReducedOrderModel.solve runs up to max_outer_iterations passes.

Stopping rules (those of the project-sandbox PODSolver, which benchmark_pod.py uses), checked at every iteration from the reduced right-hand side g = (S A V)^T S r of the current state:

converged   ||g|| < max(atol, rtol ||g_0||), g_0 from the first iteration of
            the solver's first pass (kept over restarts); or, after a step, ||da|| < max(atol,
            rtol max(||a||, 1))
stalled     ||g|| > (1 - stall_rtol) times its value stall_window
            iterations earlier
stop        non-finite ||g||, a singular reduced Jacobian, or max_it

Step: the Gauss-Newton step da limited to the positivity-preserving length, then (line_search) halved until the merit ||S r||^2 / 2 decreases sufficiently (Armijo), at most max_backtracks times; if no halving does, the positivity-limited step is taken as is. A full Gauss-Newton step can overshoot when the minimum residual in the basis is large.

S = None

V

W = None

_AV = None

_r = None

armijo_c = 0.0001

atol = 1e-12

comm

converged = False

g0 = None

last_g_norm = None

line_search = True

max_backtracks = 10

max_it = 30

mean = None

n_iterations = 0

property rb_size

rtol = 1e-10

stall_rtol = 1e-10

stall_window = 3

stalled = False

class dragonfly_sim.core.reduced_order_model.ROMResult

accepted: bool

basis_width: int

converged: bool

eta: float

iterations: int

residual_norm: float

solution: numpy.ndarray

walltime: float

class dragonfly_sim.core.reduced_order_model.ReducedOrderModel(rb_size=20, basis_mode='global', max_rank=None, qr_rank_tol=1e-10, subtract_mean=False, mean_weighting='snapshot_weights', weightfunction='Cubic', n_nonzero_weights=None, cubic_cutoff=None, parameter_scales=None, distance_weighting='total_derivative', inner_product='euclidean', state_scaling='freestream', newton_max_it=15, newton_rtol=1e-10, newton_atol=1e-12, stall_window=3, stall_rtol=1e-10, max_outer_iterations=3, g_conv_limit=1e-10, line_search=True, eta_threshold=1e-10, min_snapshots=2, run_fom_every_evaluation=False)

POD snapshots, basis and LSPG solve for DG_windtunnel_model (reduced_order_model=…). One instance per windtunnel model.

Basis:

rb_size            basis width
basis_mode         'global' or 'local_weighted' (core/pod.py)
max_rank           cap on the stored snapshot SVD rank (default 3 rb_size)
subtract_mean      local_weighted only: centre the snapshots on their weighted mean,
                   u = V a + mean (default False)
mean_weighting     'snapshot_weights' (the local distance weights, default) or
                   'uniform' (1/n per snapshot); used with subtract_mean
qr_rank_tol        relative weighted residual (outside the span of the stored
                   snapshots) below which a snapshot adds no direction
weightfunction, n_nonzero_weights, cubic_cutoff
                   local_weighted weighting (PODBasis.construct_basis);
                   cubic_cutoff = c sets the cubic support radius to
                   c min_dist + (1 - c) max_dist, widened if needed so
                   that the rb_size nearest snapshots keep a weight
parameter_scales   fixed per-parameter factors on the snapshot distance
distance_weighting output-sensitivity weights on the snapshot distance, from the
                   total derivatives of the objective and constraints
                   (utils/dv_sensitivity.py): None (off), 'total_derivative'
                   (default; outputs divided by their own magnitude) or
                   'total_derivative_raw'. They are multiplied with
                   parameter_scales. Without a derivative probe feeding
                   set_distance_sensitivities there is no information, and the
                   distance stays unweighted

Inner product (default: Euclidean with free-stream state scaling):

inner_product      'euclidean' or 'l2' (DG0 mass matrix)
state_scaling      'freestream' (model.state_scales()), None, or one
                   scale per state component

Solve and acceptance (the defaults are project-sandbox benchmark_pod.py’s):

newton_max_it, newton_rtol, newton_atol, stall_window, stall_rtol
                   per-pass LSPGSolver stopping rules
max_outer_iterations
                   LSPG passes, each restarted from the current state
g_conv_limit       the ROM solve has converged when the last ||g|| of a pass is
                   below max(newton_atol, g_conv_limit ||g_0||), g_0 the first
                   ||g|| of the first pass (the reference of newton_rtol too);
                   0 never accepts. Keep it >= newton_rtol, or the passes stop
                   before this target
line_search        Armijo damping of the Gauss-Newton step
eta_threshold      the ROM solution is accepted when ||R||/||u|| is below it
min_snapshots      snapshots needed before the ROM is tried
run_fom_every_evaluation
                   also run the FOM when the ROM is accepted (benchmarking;
                   the accepted output stays the ROM's)

After an accepted solve, adjoint is its LSPGAdjoint (None otherwise), valid until the next solve.

_check_parameter_size(parameter_vector)

_diagonal(model, mass_power, scale_power): diag(M^mass_power D^scale_power) as a Vec in the state layout, on the model’s current mesh (None when it is the identity).

_distance_weights(n_parameters): The per-parameter factors on the snapshot distance: parameter_scales times the output-sensitivity weights. None where there is neither, and the sensitivity part is dropped when it does not fit the parameter vector (e.g. alpha is not a design variable).

add_snapshot(model, parameter_vector): Add the model’s current u_vec as the snapshot at parameter_vector. Returns whether it added a basis direction. Collective.

destroy()

discard_adjoint(): Drop the adjoint of the last accepted solve. Collective.

set_distance_sensitivities(norms): Take the weights of the snapshot distance from the total-derivative norms of dv_sensitivity (the norms dict its probe passes to on_norms): ||dF/dp_j|| over the objective and constraints, normalized to mean 1. Unusable norms (non-finite, or all zero) keep the previous weights.

set_up(model): Build the inner-product weight and the POD on the model’s current (undeformed) mesh. Call once the model’s boundary conditions, function spaces and state exist. Collective.

snapshot_singular_values()

solve(model, parameter_vector): The ROM solution at parameter_vector, starting from the model’s current u_vec (projected onto the basis); u_vec holds it afterwards. Returns a ROMResult, or None when there are fewer than min_snapshots snapshots. When the solution is accepted, adjoint is its adjoint; any earlier one is discarded first. Collective.

W = None

_warned_weight_size = False

adjoint = None

basis_mode = 'global'

cubic_cutoff = None

distance_sensitivities = None

distance_weighting = 'total_derivative'

eta_threshold = 1e-10

g_conv_limit = 1e-10

inner_product = 'euclidean'

last_solver = None

line_search = True

max_outer_iterations = 3

max_rank = 60

mean_weighting = 'snapshot_weights'

min_snapshots = 2

n_nonzero_weights = None

property n_snapshots

newton_atol = 1e-12

newton_max_it = 15

newton_rtol = 1e-10

parameter_scales = None

pod = None

qr_rank_tol = 1e-10

rb_size = 20

run_fom_every_evaluation = False

scales = None

stall_rtol = 1e-10

stall_window = 3

state_scaling = 'freestream'

subtract_mean = False

weightfunction = 'Cubic'

dragonfly_sim.core.reduced_order_model.scaled_jacobian_basis(model, V, S, AV=None): S A V at the model’s current u_vec, A the exact Jacobian assembled by its Newton problem (problem.J). Returns (this rank’s rows of S A V, the PETSc product A V); AV, if given, is a product of the same shape to reuse. Without S the rows are a view of AV’s array.