Coverage for dragonfly_sim/core/reduced_order_model.py: 92%
371 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""
2POD reduced-order model (ROM) of a steady flow solve: least-squares
3Petrov-Galerkin (LSPG) on a POD basis of converged solutions.
5The ROM works with any model built on CompressibleEulerModel -- the Euler
6model itself and CompressibleRANSModel, which forwards to it -- through this
7interface, and nothing model-specific::
9 u_vec the state (its layout is the basis's row layout)
10 set_up_solver(), problem the Newton problem: problem.F(x, b) and
11 problem.J(x, problem._A) assemble the residual
12 and its exact Jacobian at the current u_vec
13 (for RANS the condensed-gradient problem, which
14 updates the reconstructed gradient on the way)
15 steady_residual() context in which problem assembles the steady
16 residual (pseudo/physical time terms off)
17 positivity_step_length step length keeping density and pressure positive
18 state_scales() per-component magnitudes for the inner product
19 _assemble_physical_residual(), relative_residual_norm()
21Inner product. Snapshots and basis use (u, v)_W = u^T W v with, by default,
22the state-scaled Euclidean weight W = D^-2, D the per-component state scales
23s_k: every component on the same footing (in absolute terms the SA variable
24rho*nu_tilde is ~1e-4 against rho*E ~ 2, and a truncated SVD would ignore
25it). inner_product='l2' adds the (diagonal) DG0 mass matrix, W = M D^-2,
26i.e. int sum_k u_k v_k / s_k^2 dx; on the NACA0012 alpha sweeps it was 10-100x
27less accurate than the Euclidean weight (M emphasizes the large far-field
28cells). W is built once on the undeformed mesh: the incremental SVD keeps
29Q^T W Q = I only for a fixed W. The LSPG residual is measured in the matching
30norm ||S r||, S = D^-1 (D^-1 M^-1/2 with 'l2', on the current mesh).
31state_scaling=None with 'euclidean' is the plain Euclidean POD/LSPG.
33Solve. With V the basis (V^T W V = I), u = V a (+ mean), LSPG is the
34Gauss-Newton iteration on min ||S R(V a)||:
36 (S A V)^T (S A V) da = -(S A V)^T S r, a <- a + theta da,
38theta from positivity_step_length, optionally backtracked (Armijo) on
39||S r||^2 / 2. The reduced system is small and solved redundantly on every
40rank. Convergence follows project-sandbox benchmark_pod.py (its PODSolver and
41outer loop): passes of at most newton_max_it iterations, each stopping on
42||g|| = ||(S A V)^T S r|| < max(atol, rtol ||g_0||), on a negligible step, or on
43a stall of ||g||; up to max_outer_iterations passes, converged when a pass
44ends with ||g|| < g_conv_limit.
46Derivatives. The ROM solves the LSPG stationarity condition, the projected
47residual equations
49 Phi^T R(V a; x) = 0, Phi = S^T S A V,
51so when its solution is accepted the adjoint is that of these k equations
52(ReducedOrderModel.adjoint, an LSPGAdjoint). With V, S and Phi frozen at the
53solution (the basis is not differentiated, and the Gauss-Newton terms, which
54are proportional to S R, are dropped), for an output f with g = (df/du)^T:
56 (S A V)^T (S A V) lam = V^T g, psi = Phi lam = S^T S A V lam,
58and df/dx = df/dx|_u - psi^T dR/dx with the full-order dR/dx: psi takes the
59place of the full-order adjoint A^-T g. It is the Petrov-Galerkin solution of
60A^T psi = g in the span of Phi, tested with V: V^T (A^T psi - g) = 0. psi is
61computed as S Q y with S A V = Q R (a distributed thin QR) and R^T y = V^T g,
62without forming the normal matrix, whose condition number is squared.
63"""
64from dataclasses import dataclass
65from time import perf_counter
67import numpy as np
68import scipy.linalg as sp_la
69from mpi4py import MPI
70from petsc4py import PETSc
72from dragonfly_sim.core.pod import PODBasis
73from dragonfly_sim.utils.pod_utils import DiagonalWeight, tsqr
74from dragonfly_sim.utils.derived_fields import diagonal_mass_inverse
77def scaled_jacobian_basis(model, V, S, AV=None):
78 """
79 S A V at the model's current u_vec, A the exact Jacobian assembled by its
80 Newton problem (problem.J). Returns (this rank's rows of S A V, the PETSc
81 product A V); `AV`, if given, is a product of the same shape to reuse.
82 Without S the rows are a view of AV's array.
83 """
84 problem = model.problem
85 problem.J(model.u_vec.x.petsc_vec, problem._A)
86 A = problem._A
87 if AV is None:
88 AV = A.matMult(V)
89 else:
90 A.matMult(V, result=AV)
91 SAV = AV.getDenseArray()
92 if S is not None:
93 SAV = S.getArray(readonly=True)[:, None] * SAV
94 return SAV, AV
97class LSPGSolver:
98 """
99 Reduced Newton (Gauss-Newton) loop on a POD basis, driving a model's
100 Newton problem (see the module docstring).
102 Parameters
103 ----------
104 V : PETSc.Mat
105 The basis, a distributed MATDENSE in the state layout.
106 W : PETSc.Vec or None
107 The diagonal inner-product weight, for the restriction
108 a = V^T W (u - mean); None: Euclidean (V^T V = I).
109 S : PETSc.Vec or None
110 The diagonal residual weight; None: the plain residual norm.
111 mean : PETSc.Vec or None
112 Optional affine offset (u = V a + mean).
113 """
115 def __init__(self, V, comm, W=None, S=None, mean=None, rtol=1e-10, atol=1e-12, max_it=30,
116 stall_window=3, stall_rtol=1e-10, line_search=True, armijo_c=1e-4, max_backtracks=10):
117 self.V = V
118 self.comm = comm
119 self.W = W
120 self.S = S
121 self.mean = mean
122 self.rtol = rtol
123 self.atol = atol
124 self.max_it = max_it
125 self.stall_window = stall_window
126 self.stall_rtol = stall_rtol
127 self.line_search = line_search
128 self.armijo_c = armijo_c
129 self.max_backtracks = max_backtracks
130 self.converged = False
131 self.stalled = False
132 self.n_iterations = 0
133 self.last_g_norm = None
134 self._AV = None
135 self._r = None
137 @property
138 def rb_size(self):
139 return self.V.getSize()[1]
141 def destroy(self):
142 for name in ("_AV", "_r"):
143 obj = getattr(self, name)
144 if obj is not None:
145 obj.destroy()
146 setattr(self, name, None)
148 # -- maps between the full and reduced spaces --------------------------
149 def restrict(self, u_vec):
150 """a = V^T W (u - mean), replicated on every rank."""
151 x = u_vec.getArray(readonly=True).copy()
152 if self.mean is not None:
153 x -= self.mean.getArray(readonly=True)
154 if self.W is not None:
155 x *= self.W.getArray(readonly=True)
156 return self.comm.allreduce(self.V.getDenseArray().T @ x, op=MPI.SUM)
158 def expand(self, a, u):
159 """u <- V a (+ mean), ghosts updated."""
160 x = u.x.petsc_vec
161 x.array[:] = self.V.getDenseArray() @ a
162 if self.mean is not None:
163 x.axpy(1.0, self.mean)
164 u.x.scatter_forward()
166 # -- the reduced Newton loop ---------------------------------------------
167 def _weighted_residual(self, model):
168 """Assemble r at the current u_vec; returns (S r local entries, ||S r||)."""
169 problem = model.problem
170 if self._r is None:
171 self._r = problem._b.duplicate()
172 problem.F(model.u_vec.x.petsc_vec, self._r)
173 r = self._r.getArray(readonly=True)
174 if self.S is not None:
175 r = self.S.getArray(readonly=True) * r
176 return r, float(np.sqrt(self.comm.allreduce(float(r @ r), op=MPI.SUM)))
178 def _reduced_system(self, model):
179 """(H, g, ||S r||) with H = (S A V)^T (S A V), g = (S A V)^T S r at the
180 current u_vec."""
181 r, r_norm = self._weighted_residual(model)
182 if self._AV is not None and self._AV.getSize()[1] != self.V.getSize()[1]:
183 self._AV.destroy()
184 self._AV = None
185 SAV, self._AV = scaled_jacobian_basis(model, self.V, self.S, self._AV)
186 H = self.comm.allreduce(SAV.T @ SAV, op=MPI.SUM)
187 g = self.comm.allreduce(SAV.T @ r, op=MPI.SUM)
188 return H, g, r_norm
190 def solve(self, model):
191 """
192 One pass of the LSPG Gauss-Newton loop: advance model.u_vec, starting
193 from its restriction onto the basis, towards the LSPG solution in the
194 span of V. Call inside model.steady_residual() for a steady solve.
195 Returns (n_iterations, converged). ReducedOrderModel.solve runs up to
196 max_outer_iterations passes.
198 Stopping rules (those of the project-sandbox PODSolver, which
199 benchmark_pod.py uses), checked at every iteration from the reduced
200 right-hand side g = (S A V)^T S r of the current state::
202 converged ||g|| < max(atol, rtol ||g_0||), g_0 from this pass's first
203 iteration; or, after a step, ||da|| < max(atol,
204 rtol max(||a||, 1))
205 stalled ||g|| > (1 - stall_rtol) times its value stall_window
206 iterations earlier
207 stop non-finite ||g||, a singular reduced Jacobian, or max_it
209 Step: the Gauss-Newton step da limited to the positivity-preserving
210 length, then (line_search) halved until the merit ||S r||^2 / 2
211 decreases sufficiently (Armijo), at most max_backtracks times; if no
212 halving does, the positivity-limited step is taken as is. A full
213 Gauss-Newton step can overshoot when the minimum residual in the
214 basis is large.
215 """
216 u = model.u_vec
217 a = self.restrict(u.x.petsc_vec)
218 dx = u.x.petsc_vec.duplicate()
219 self.converged = self.stalled = False
220 g0 = None
221 history = []
222 n = 0
223 PETSc.Sys.Print("LSPG ROM solve, basis width {}".format(self.rb_size))
224 for n in range(self.max_it):
225 self.expand(a, u)
226 H, g, r_norm = self._reduced_system(model)
227 g_norm = float(np.linalg.norm(g))
228 self.last_g_norm = g_norm
229 g0 = g_norm if g0 is None else g0
230 if not np.isfinite(g_norm):
231 PETSc.Sys.Print("ROM: non-finite reduced residual at iteration {}; stopping".format(n))
232 break
233 if g_norm < max(self.atol, self.rtol * g0):
234 self.converged = True
235 PETSc.Sys.Print("ROM converged at iteration {}: |g| {:.6e} (from {:.6e})".format(n, g_norm, g0))
236 break
237 history.append(g_norm)
238 if len(history) > self.stall_window and \
239 g_norm > (1.0 - self.stall_rtol) * history[-(self.stall_window + 1)]:
240 self.stalled = True
241 PETSc.Sys.Print("ROM stalled at iteration {}: |g| {:.6e}, vs {:.6e} {} iterations "
242 "earlier".format(n, g_norm, history[-(self.stall_window + 1)],
243 self.stall_window))
244 break
245 try:
246 da = np.linalg.solve(H, -g)
247 except np.linalg.LinAlgError:
248 PETSc.Sys.Print("ROM: singular reduced Jacobian at iteration {}".format(n))
249 break
250 # the full-space update is +V da; positivity_step_length limits x - theta dx
251 dx.array[:] = -(self.V.getDenseArray() @ da)
252 theta, _ = model.positivity_step_length(u.x.petsc_vec, dx)
253 if self.line_search:
254 theta = self._backtrack(model, a, da, theta, 0.5 * r_norm**2, float(g @ da))
255 a = a + theta * da
256 da_norm = float(np.linalg.norm(da))
257 PETSc.Sys.Print("ROM iter {}: |g| = {:.6e}, ||S r|| = {:.6e}, |da| = {:.6e}, theta = {:.4f}".format(
258 n, g_norm, r_norm, da_norm, theta))
259 if da_norm < max(self.atol, self.rtol * max(float(np.linalg.norm(a)), 1.0)):
260 self.converged = True
261 PETSc.Sys.Print("ROM converged at iteration {}: |da| {:.6e}".format(n, da_norm))
262 break
263 self.expand(a, u)
264 dx.destroy()
265 if not self.converged:
266 PETSc.Sys.Print("[warning] ROM pass stopped after {} iterations without converging{}".format(
267 n + 1, " (stalled)" if self.stalled else ""))
268 self.n_iterations = n + 1
269 return self.n_iterations, self.converged
271 def _backtrack(self, model, a, da, theta, phi, slope):
272 """Largest theta * 2^-k (k <= max_backtracks) with phi(a + theta da)
273 <= phi + armijo_c theta slope; the given theta if there is none."""
274 t = theta
275 for _ in range(self.max_backtracks + 1):
276 self.expand(a + t * da, model.u_vec)
277 _, r_trial = self._weighted_residual(model)
278 if np.isfinite(r_trial) and 0.5 * r_trial**2 <= phi + self.armijo_c * t * slope:
279 return t
280 t *= 0.5
281 PETSc.Sys.Print("ROM: no Armijo step along da; taking theta = {:.4f}".format(theta))
282 return theta
285class LSPGAdjoint:
286 """
287 The adjoint of an accepted LSPG solution (Derivatives in the module
288 docstring): psi = S Q y, with S A V = Q R and R^T y = V^T g, for the basis
289 V and residual weight S of that solve. The factorization is formed by the
290 first solve() and reused for every later right-hand side (CSDL passes one
291 per output row). It belongs to one design and state; matches() checks
292 both.
294 Parameters
295 ----------
296 V : PETSc.Mat
297 The basis of the solve (owned by the POD, not by this object).
298 S : PETSc.Vec or None
299 The diagonal residual weight of the solve; destroyed with this object.
300 parameter_vector, solution
301 The design and this rank's owned entries of the accepted solution.
302 """
304 def __init__(self, V, S, comm, parameter_vector, solution):
305 self.V = V
306 self.S = S
307 self.comm = comm
308 self.parameter_vector = np.array(parameter_vector, dtype=float).ravel()
309 self.solution = np.array(solution, dtype=float)
310 self._Q = None
311 self._R = None
313 @property
314 def rb_size(self):
315 return self.V.getSize()[1]
317 def matches(self, parameter_vector, u_local):
318 """Whether this is the adjoint at design `parameter_vector` and state
319 `u_local` (owned entries): both equal to the accepted solve's on every
320 rank. Collective."""
321 pv = np.asarray(parameter_vector, dtype=float).ravel()
322 local = bool(np.array_equal(pv, self.parameter_vector)
323 and np.array_equal(np.asarray(u_local, dtype=float), self.solution))
324 return self.comm.allreduce(local, op=MPI.LAND)
326 def _factor(self, model):
327 """Q, R of S A V, A the steady Jacobian at the model's current u_vec."""
328 model.set_up_solver()
329 with model.steady_residual():
330 SAV, AV = scaled_jacobian_basis(model, self.V, self.S)
331 self._Q, self._R = tsqr(self.comm, SAV)
332 AV.destroy()
333 diag = np.abs(np.diag(self._R))
334 if not np.all(np.isfinite(self._R)) or diag.min() <= np.finfo(float).eps * diag.max():
335 raise ValueError("reduced adjoint: S A V is singular to working precision "
336 "(|R_ii| from {:.3e} to {:.3e})".format(diag.min(), diag.max()))
338 def solve(self, model, rhs_local):
339 """
340 psi (this rank's owned entries) for the right-hand side g, rhs_local
341 being its owned entries. The model's u_vec, mesh and angle of attack
342 must be those of the accepted solution (the first call assembles the
343 Jacobian there). Collective.
344 """
345 first = self._Q is None
346 if first:
347 self._factor(model)
348 rhs = np.asarray(rhs_local, dtype=float)
349 b = self.comm.allreduce(self.V.getDenseArray().T @ rhs, op=MPI.SUM)
350 y = sp_la.solve_triangular(self._R, b, trans='T', lower=False)
351 psi = self._Q @ y
352 if self.S is not None:
353 psi = self.S.getArray(readonly=True) * psi
354 psi_norm = np.sqrt(self.comm.allreduce(float(psi @ psi), op=MPI.SUM))
355 if not np.isfinite(psi_norm):
356 raise ValueError("reduced adjoint: non-finite solution")
357 if first:
358 diag = np.abs(np.diag(self._R))
359 PETSc.Sys.Print("reduced (LSPG) adjoint: basis width {}, |R_ii| from {:.3e} to {:.3e}, "
360 "Petrov-Galerkin residual |V^T (A^T psi - g)| / |V^T g| = {:.3e}".format(
361 self.rb_size, diag.min(), diag.max(),
362 self._galerkin_residual(model, psi, rhs, b)))
363 else:
364 PETSc.Sys.Print("reduced (LSPG) adjoint: reused the factorization, |psi| = {:.6e}".format(psi_norm))
365 return psi
367 def _galerkin_residual(self, model, psi, rhs, b):
368 """|V^T (A^T psi - g)| / |V^T g|, with the Jacobian _factor assembled
369 (zero up to round-off)."""
370 A = model.problem._A
371 x, Atx = A.createVecLeft(), A.createVecRight()
372 x.array[:] = psi
373 A.multTranspose(x, Atx)
374 res = self.comm.allreduce(self.V.getDenseArray().T @ (Atx.array - rhs), op=MPI.SUM)
375 x.destroy()
376 Atx.destroy()
377 return float(np.linalg.norm(res) / max(np.linalg.norm(b), 1e-300))
379 def destroy(self):
380 if self.S is not None:
381 self.S.destroy()
382 self.S = None
383 self._Q = self._R = None
386@dataclass
387class ROMResult:
388 solution: np.ndarray # this rank's owned entries of the ROM state
389 eta: float # ||R|| / ||u|| of the full steady residual
390 residual_norm: float
391 converged: bool # the reduced Newton loop converged
392 accepted: bool # eta < eta_threshold
393 iterations: int
394 basis_width: int
395 walltime: float
398class ReducedOrderModel:
399 """
400 POD snapshots, basis and LSPG solve for DG_windtunnel_model
401 (reduced_order_model=...). One instance per windtunnel model.
403 Basis::
405 rb_size basis width
406 basis_mode 'global' or 'local_weighted' (core/pod.py)
407 max_rank cap on the stored snapshot SVD rank (default 3 rb_size)
408 qr_rank_tol relative weighted residual (outside the span of the stored
409 snapshots) below which a snapshot adds no direction
410 weightfunction, n_nonzero_weights, cubic_cutoff
411 local_weighted weighting (PODBasis.construct_basis);
412 cubic_cutoff = c sets the cubic support radius to
413 c min_dist + (1 - c) max_dist, widened if needed so
414 that the rb_size nearest snapshots keep a weight
415 parameter_scales fixed per-parameter factors on the snapshot distance
416 distance_weighting output-sensitivity weights on the snapshot distance, from the
417 total derivatives of the objective and constraints
418 (utils/dv_sensitivity.py): None (off), 'total_derivative'
419 (default; outputs divided by their own magnitude) or
420 'total_derivative_raw'. They are multiplied with
421 parameter_scales. Without a derivative probe feeding
422 set_distance_sensitivities there is no information, and the
423 distance stays unweighted
425 Inner product (default: Euclidean with free-stream state scaling)::
427 inner_product 'euclidean' or 'l2' (DG0 mass matrix)
428 state_scaling 'freestream' (model.state_scales()), None, or one
429 scale per state component
431 Solve and acceptance (the defaults are project-sandbox benchmark_pod.py's)::
433 newton_max_it, newton_rtol, newton_atol, stall_window, stall_rtol
434 per-pass LSPGSolver stopping rules
435 max_outer_iterations
436 LSPG passes, each restarted from the current state
437 g_conv_limit the ROM solve has converged when the last ||g|| of a
438 pass is below this
439 line_search Armijo damping of the Gauss-Newton step
440 eta_threshold the ROM solution is accepted when ||R||/||u|| is below it
441 min_snapshots snapshots needed before the ROM is tried
442 run_fom_every_evaluation
443 also run the FOM when the ROM is accepted (benchmarking;
444 the accepted output stays the ROM's)
446 After an accepted solve, `adjoint` is its LSPGAdjoint (None otherwise),
447 valid until the next solve.
448 """
450 def __init__(self, rb_size=20, basis_mode='global', max_rank=None, qr_rank_tol=1e-10,
451 weightfunction='Cubic', n_nonzero_weights=None, cubic_cutoff=None,
452 parameter_scales=None, distance_weighting='total_derivative',
453 inner_product='euclidean', state_scaling='freestream',
454 newton_max_it=15, newton_rtol=1e-10, newton_atol=1e-12, stall_window=3,
455 stall_rtol=1e-10, max_outer_iterations=3, g_conv_limit=1e-12, line_search=True,
456 eta_threshold=1e-10, min_snapshots=2, run_fom_every_evaluation=False):
457 if inner_product not in ('l2', 'euclidean'):
458 raise ValueError("inner_product must be 'l2' or 'euclidean', got {!r}".format(inner_product))
459 if not (state_scaling is None or isinstance(state_scaling, str) and state_scaling == 'freestream'
460 or np.ndim(state_scaling) == 1):
461 raise ValueError("state_scaling must be 'freestream', None or one scale per state component")
462 self.rb_size = rb_size
463 self.basis_mode = basis_mode
464 self.max_rank = 3 * rb_size if max_rank is None else max_rank
465 self.qr_rank_tol = qr_rank_tol
466 self.weightfunction = weightfunction
467 self.n_nonzero_weights = n_nonzero_weights
468 self.cubic_cutoff = cubic_cutoff
469 self.parameter_scales = parameter_scales
470 if distance_weighting not in (None, 'total_derivative', 'total_derivative_raw'):
471 raise ValueError("distance_weighting must be None, 'total_derivative' or "
472 "'total_derivative_raw', got {!r}".format(distance_weighting))
473 self.distance_weighting = distance_weighting
474 # Latest mean-normalized weights, written by set_distance_sensitivities;
475 # None until a gradient evaluation has been probed. They lag the basis
476 # they weight by one design point.
477 self.distance_sensitivities = None
478 self._warned_weight_size = False
479 self.inner_product = inner_product
480 self.state_scaling = state_scaling
481 self.newton_max_it = newton_max_it
482 self.newton_rtol = newton_rtol
483 self.newton_atol = newton_atol
484 self.stall_window = stall_window
485 self.stall_rtol = stall_rtol
486 self.max_outer_iterations = max_outer_iterations
487 self.g_conv_limit = g_conv_limit
488 self.line_search = line_search
489 self.eta_threshold = eta_threshold
490 self.min_snapshots = min_snapshots
491 self.run_fom_every_evaluation = run_fom_every_evaluation
493 self.pod = None
494 self.W = None
495 self.scales = None
496 self.last_solver = None
497 self.adjoint = None
499 # ------------------------------------------------------------------
500 # Set-up and weights
501 # ------------------------------------------------------------------
502 def set_up(self, model):
503 """Build the inner-product weight and the POD on the model's
504 current (undeformed) mesh. Call once the model's boundary conditions,
505 function spaces and state exist. Collective."""
506 layout = model.u_vec.x.petsc_vec
507 if self.state_scaling is None:
508 self.scales = None
509 elif isinstance(self.state_scaling, str):
510 self.scales = np.asarray(model.state_scales(), dtype=float)
511 else:
512 self.scales = np.asarray(self.state_scaling, dtype=float)
513 if self.scales is not None and self.scales.shape != (model.n_state,):
514 raise ValueError("{} state scales for {} state components".format(
515 self.scales.shape[0], model.n_state))
516 if self.scales is not None and not np.all(self.scales > 0):
517 raise ValueError("state scales must be positive, got {}".format(self.scales))
519 if self.W is not None:
520 self.W.destroy()
521 w = self._diagonal(model, mass_power=1.0, scale_power=-2.0)
522 self.W = None if w is None else DiagonalWeight(w)
523 if self.pod is not None:
524 self.pod.destroy()
525 self.pod = PODBasis(model.mesh_obj.mesh.comm, layout.getSize(), basis_mode=self.basis_mode,
526 IP_matrix=self.W, max_rank=self.max_rank, qr_rank_tol=self.qr_rank_tol)
527 PETSc.Sys.Print("Initialized POD ROM: basis_mode={}, rb_size={}, max_rank={}, n_dofs={}, "
528 "inner product {}{}".format(
529 self.basis_mode, self.rb_size, self.max_rank, layout.getSize(),
530 self.inner_product,
531 "" if self.scales is None else ", state scales {}".format(self.scales)))
533 def set_distance_sensitivities(self, norms):
534 """Take the weights of the snapshot distance from the total-derivative
535 norms of dv_sensitivity (the `norms` dict its probe passes to
536 `on_norms`): ||dF/dp_j|| over the objective and constraints, normalized
537 to mean 1. Unusable norms (non-finite, or all zero) keep the previous
538 weights."""
539 if self.distance_weighting is None:
540 return
541 weights = norms['column_norms_relative'] if self.distance_weighting == 'total_derivative' \
542 else norms['column_norms']
543 if weights is None:
544 # no output values in the derivative result
545 PETSc.Sys.Print("ROM distance weighting: output values unavailable, using the raw "
546 "total-derivative norms for this evaluation")
547 weights = norms['column_norms']
548 weights = np.asarray(weights, dtype=float).ravel()
549 if not np.all(np.isfinite(weights)) or not np.any(weights > 0.):
550 PETSc.Sys.Print("ROM distance weighting: the total-derivative norms are all zero or "
551 "not finite; keeping the previous weights")
552 return
553 self.distance_sensitivities = weights / weights.mean()
554 inert = "" if self.basis_mode == 'local_weighted' else \
555 " -- NOT IN USE, basis_mode is '{}'".format(self.basis_mode)
556 PETSc.Sys.Print("ROM distance weights ({}, {} parameters, mean-normalized): min {:.4e}, "
557 "median {:.4e}, max {:.4e}{}".format(
558 self.distance_weighting, weights.shape[0],
559 self.distance_sensitivities.min(), np.median(self.distance_sensitivities),
560 self.distance_sensitivities.max(), inert))
562 def _distance_weights(self, n_parameters):
563 """The per-parameter factors on the snapshot distance: parameter_scales
564 times the output-sensitivity weights. None where there is neither, and
565 the sensitivity part is dropped when it does not fit the parameter
566 vector (e.g. alpha is not a design variable)."""
567 scales = None if self.parameter_scales is None else np.asarray(self.parameter_scales, dtype=float)
568 weights = self.distance_sensitivities
569 if self.distance_weighting is None or weights is None:
570 return scales
571 if weights.shape[0] != n_parameters:
572 if not self._warned_weight_size:
573 PETSc.Sys.Print("ROM distance weighting: {} weights do not fit the {}-entry parameter "
574 "vector; using the distance without them. They come from the design "
575 "variables of the CSDL problem, which must be the snapshot "
576 "parameters.".format(weights.shape[0], n_parameters))
577 self._warned_weight_size = True
578 return scales
579 return weights if scales is None else scales * weights
581 def _diagonal(self, model, mass_power, scale_power):
582 """diag(M^mass_power D^scale_power) as a Vec in the state layout, on
583 the model's current mesh (None when it is the identity)."""
584 if self.inner_product == 'euclidean' and self.scales is None:
585 return None
586 V = model.functionspaces["V"]
587 if self.inner_product == 'l2':
588 d = diagonal_mass_inverse(V, model.forms.dx)
589 arr = d.array
590 arr[:] = arr ** (-mass_power)
591 else:
592 d = model.u_vec.x.petsc_vec.duplicate()
593 d.set(1.0)
594 arr = d.array
595 if self.scales is not None:
596 # owned entries come in blocks of n_state, one per cell
597 arr *= np.tile(self.scales ** scale_power, arr.shape[0] // model.n_state)
598 return d
600 @property
601 def n_snapshots(self):
602 return 0 if self.pod is None else self.pod.snapshot_matrix.n_snapshots
604 def snapshot_singular_values(self):
605 sigma = None if self.pod is None else self.pod.snapshot_matrix.Sigma
606 return None if sigma is None else sigma.getArray().copy()
608 # ------------------------------------------------------------------
609 # Snapshots and solve
610 # ------------------------------------------------------------------
611 def add_snapshot(self, model, parameter_vector):
612 """Add the model's current u_vec as the snapshot at parameter_vector.
613 Returns whether it added a basis direction. Collective."""
614 added = self.pod.snapshot_matrix.process_and_append_to_snapshot_lists(
615 model.u_vec.x.petsc_vec, np.asarray(parameter_vector, dtype=float))
616 self.pod.snapshot_matrix.print_info()
617 return added
619 def discard_adjoint(self):
620 """Drop the adjoint of the last accepted solve. Collective."""
621 if self.adjoint is not None:
622 self.adjoint.destroy()
623 self.adjoint = None
625 def solve(self, model, parameter_vector):
626 """
627 The ROM solution at parameter_vector, starting from the model's
628 current u_vec (projected onto the basis); u_vec holds it afterwards.
629 Returns a ROMResult, or None when there are fewer than min_snapshots
630 snapshots. When the solution is accepted, `adjoint` is its adjoint;
631 any earlier one is discarded first. Collective.
632 """
633 self.discard_adjoint()
634 if self.n_snapshots < self.min_snapshots:
635 PETSc.Sys.Print("ROM: {} snapshot(s) stored, {} needed; skipping the ROM solve".format(
636 self.n_snapshots, self.min_snapshots))
637 return None
638 t0 = perf_counter()
639 comm = model.mesh_obj.mesh.comm
640 self.pod.construct_basis(rb_size=self.rb_size,
641 local_parametervector=np.asarray(parameter_vector, dtype=float),
642 parameter_weights=self._distance_weights(len(parameter_vector)),
643 weightfunction=self.weightfunction,
644 n_nonzero_weights=self.n_nonzero_weights,
645 cubic_cutoff=self.cubic_cutoff,
646 layout_vec=model.u_vec.x.petsc_vec)
647 if self.pod.basis is None or self.pod.basis.getSize()[1] == 0:
648 PETSc.Sys.Print("ROM: empty POD basis; skipping the ROM solve")
649 return None
651 model.set_up_solver()
652 S = self._diagonal(model, mass_power=-0.5, scale_power=-1.0)
653 solver = LSPGSolver(self.pod.basis, comm, W=None if self.W is None else self.W.w, S=S,
654 mean=self.pod.mean, rtol=self.newton_rtol, atol=self.newton_atol,
655 max_it=self.newton_max_it, stall_window=self.stall_window,
656 stall_rtol=self.stall_rtol, line_search=self.line_search)
657 # Up to max_outer_iterations passes, each restarted from the current
658 # state; converged when the pass's last |g| is below g_conv_limit (the
659 # LSPG stationarity residual -- the projected residual V^T r does not
660 # vanish at the LSPG solution, so it is not tested)
661 converged = False
662 iterations = 0
663 residual_norm = np.nan
664 try:
665 for _ in range(self.max_outer_iterations):
666 with model.steady_residual():
667 its, _ = solver.solve(model)
668 iterations += its
669 residual_norm = model._assemble_physical_residual()
670 if not np.isfinite(residual_norm):
671 PETSc.Sys.Print("ROM: non-finite residual; stopping")
672 break
673 if solver.last_g_norm is not None and solver.last_g_norm < self.g_conv_limit:
674 converged = True
675 PETSc.Sys.Print("LSPG converged: |g| = {:.6e} < {:.1e}".format(
676 solver.last_g_norm, self.g_conv_limit))
677 break
678 else:
679 PETSc.Sys.Print("ROM: {} LSPG passes without |g| < {:.1e} (last |g| {:.6e})".format(
680 self.max_outer_iterations, self.g_conv_limit, solver.last_g_norm))
681 except BaseException:
682 if S is not None:
683 S.destroy()
684 raise
685 finally:
686 solver.destroy()
687 self.last_solver = solver
688 solver.n_iterations = iterations
690 eta = model.relative_residual_norm(residual_norm)
691 model.u_vec.x.scatter_forward()
692 result = ROMResult(solution=model.u_vec.x.petsc_vec.getArray().copy(), eta=eta,
693 residual_norm=residual_norm, converged=converged,
694 accepted=bool(np.isfinite(eta) and eta < self.eta_threshold),
695 iterations=iterations, basis_width=solver.rb_size,
696 walltime=perf_counter() - t0)
697 PETSc.Sys.Print("ROM relative residual norm eta: {} (||r||: {}); {} iterations, wall time {}".format(
698 eta, residual_norm, iterations, result.walltime))
699 if result.accepted:
700 # the basis stays valid until the next solve rebuilds it, which
701 # discards this adjoint first
702 self.adjoint = LSPGAdjoint(self.pod.basis, S, comm, parameter_vector, result.solution)
703 elif S is not None:
704 S.destroy()
705 return result
707 def destroy(self):
708 self.discard_adjoint()
709 if self.pod is not None:
710 self.pod.destroy()
711 self.pod = None
712 if self.W is not None:
713 self.W.destroy()
714 self.W = None