Coverage for dragonfly_sim/core/windtunnel_model.py: 84%

1import numpy as np 

2from mpi4py import MPI 

3from petsc4py import PETSc 

4import dolfinx 

5 

6import datetime 

7import resource 

8import os 

9import gc 

10 

11import csdl_alpha as csdl 

12 

13from time import perf_counter 

14 

15from dragonfly_sim.core.mesh_manager import Mesh 

16from dragonfly_sim.core.Euler_model import CompressibleEulerModel 

17from dragonfly_sim.core.meshwarping import IDWarp_jax 

18from dragonfly_sim.utils.data_handler import DataStore 

19 

20from dragonfly_sim.utils.petsc_utils import set_petsc_vec_array 

21from dragonfly_sim.utils.derived_fields import destroy_chained_operator 

22from dragonfly_sim.utils.solver_utils import ksp_converged_reason_name 

23from dragonfly_sim.core.postprocessor import DG_postprocessor 

24 

25 

26# TODO: Look into using entropy variables instead of conservative variables; 

27# these would potentially allow us to do simulations from low Mach all 

28# the way to supersonic 

29class DG_windtunnel_model(csdl.experimental.CustomImplicitOperation): 

30 # The solve/adjoint callbacks contain MPI collectives (PETSc, mpi4py) and deform 

31 # the shared mesh, so they must never overlap -- also on a single process. 

32 ordered_callbacks = True 

33 

34 def __init__(self, mesh, boundary_dict, shape_parameterization=None, moments=None, 

35 mesh_inner_bdry_function=None, wall_tags=None, symmetry_tags=None, 

36 poly_order=1, gamma=1.4, 

37 filename_suffix="test", 

38 asm_overlap=None, ilu_levels=None, 

39 model_class=CompressibleEulerModel, model_kwargs=None, 

40 reduced_order_model=None): 

41 """ 

42 boundary_dict: the free-stream state {'rho', 'p', 'M', 'alpha'} (density, 

43 pressure, Mach number, angle of attack in rad). The whole outer boundary 

44 is a characteristic far field built from it, so the angle of attack can 

45 change by any amount during an optimization; there is no inflow/outflow 

46 division to choose or to keep valid. 

47 

48 The shape of the object being optimized (the wall) is defined in one of 

49 three ways, in this order of precedence: 

50 1. `mesh_inner_bdry_function`, a function of the boundary facet midpoints 

51 (shape (dim, n)) returning a boolean wall mask. It can encode 

52 exclusions the FFD block cannot. If the mesh also carries boundary 

53 facet tags (utils/mesh_io_utils.py:load_dolfinx_mesh stores them on 

54 the mesh), a warning lists them: they could be used instead. 

55 2. `wall_tags`: the mesh's boundary facet tags (values or names, one or 

56 several) that are the wall, and `symmetry_tags` those that are 

57 symmetry planes (slip walls). Every other boundary is far field. 

58 In 3D without tags, the symmetry plane is y = 0. 

59 3. The shape parameterization's FFD block corners 

60 (shape_parameterization.inner_boundary_function()), which bound the 

61 box the inner boundary facets are contained in. 

62 At least one is required. 

63 

64 shape_parameterization: the wall-shape parameterization, e.g. a 

65 shape_design.FFDShapeParameterization -- the FFD block and the layers 

66 that define the shape design variables and their bounds. The model does 

67 not need to know what those variables are: set_up_sim builds the IDWarp 

68 mesh warper (self.mesh_warper) and hands the wall to 

69 shape_parameterization.setup(); deform_mesh() turns its wall 

70 displacement into mesh node motions; evaluate() takes those and the 

71 flat shape parameter vector: shape optimization. Without it (None) no 

72 warper is built, the mesh stays fixed, and evaluate() takes alpha only: 

73 a forward flow analysis, with derivatives with respect to alpha. 

74 solve_forward() runs one flow solve without a CSDL graph at all. 

75 

76 moments: the moments of the wall force to compute, print and differentiate, 

77 a sequence of (name, point, direction[, ref_length]) tuples (see 

78 DG_postprocessor); default none. E.g. the pitching moment about the quarter 

79 chord of a wing, ("pitch", (0.25, 0., 0.), "pitch"), and its root 

80 bending moment, ("root_bending", (0.25, 0., 0.), (1., 0., 0.)). 

81 

82 model_class / model_kwargs select the flow model: CompressibleEulerModel 

83 (the default) or e.g. core.RANS_model.CompressibleRANSModel with 

84 model_kwargs={'Re': 6e6}. The body's wall is the model's 

85 define_wall_bc (slip for Euler, adiabatic no-slip for RANS). 

86 

87 reduced_order_model: a core.reduced_order_model.ReducedOrderModel to 

88 try before every flow solve (see _solve), or None for full-order 

89 solves only. Where its solution is accepted, the derivatives use the 

90 adjoint of its projected equations (see apply_inverse_jacobian). 

91 """ 

92 super().__init__() 

93 

94 self.asm_overlap = asm_overlap 

95 self.ilu_levels = ilu_levels 

96 

97 # --- Linear-solve convergence criteria ---------------------------- 

98 # We use one set of convergence criteria for all linear solves. 

99 # This does not cover the Newton solver for the Euler model itself, 

100 # but all other things (adjoint, Helmholtz, elasticity-based  

101 # mesh warping, etc.) when used. 

102 self.linear_solver_rtol = 1e-20 

103 self.linear_solver_atol = 1e-11 

104 self.linear_solver_max_it = 1000 

105 

106 # If linear solves converge above this tolerance, we raise an error 

107 self.linear_residual_accept_limit = 1e-10 

108 

109 # --- Nonlinear (Euler) convergence criteria ----------------------- 

110 self.euler_residual_conv_limit = 1e-9 

111 # Total Newton-step budget for the FOM solve 

112 self.max_newton_iterations = 100 

113 # KDTree coincidence radius for the dof <-> mesh-node permutation. 

114 # These two orderings are 1:1 by construction and differ only in 

115 # index order, so this tolerance accounts for round-off errors 

116 self.mesh_node_match_tol = 1e-10 

117 

118 if shape_parameterization is not None and not hasattr(shape_parameterization, 'setup'): 

119 raise TypeError("shape_parameterization must be a shape parameterization object " 

120 "(e.g. shape_design.FFDShapeParameterization), got {!r}. The FFD " 

121 "block arguments (ffd_shape, ffd_degree, ffd_block_corner_list, " 

122 "cp_coord_opt_idxs) moved there.".format(shape_parameterization)) 

123 if mesh_inner_bdry_function is None and wall_tags is None and shape_parameterization is not None: 

124 mesh_inner_bdry_function = shape_parameterization.inner_boundary_function() 

125 if mesh_inner_bdry_function is None and wall_tags is None: 

126 raise ValueError("Provide either mesh_inner_bdry_function, wall_tags (the mesh's own " 

127 "boundary tags) or a shape parameterization with FFD block corners " 

128 "to locate the inner boundary facets.") 

129 

130 # Initial instantiation of mesh and simulation model 

131 self.mesh = Mesh(mesh, mesh_inner_bdry_function, 

132 node_match_tol=self.mesh_node_match_tol, 

133 wall_tags=wall_tags, symmetry_tags=symmetry_tags) 

134 if model_class is CompressibleEulerModel: 

135 _euler_kwargs = dict(poly_order=poly_order, gamma=gamma) 

136 self.sim_model = CompressibleEulerModel(self.mesh, **_euler_kwargs) 

137 else: 

138 self.sim_model = model_class(self.mesh, gamma=gamma, **(model_kwargs or {})) 

139 # TODO: Perhaps split up CompressibleEulerModel into an Euler model class and a solver class that facilitates solver interactions;  

140 # the solver has become so complicated recently that spinning it off into its own class would likely help with code readability 

141 

142 if asm_overlap is not None: 

143 self.sim_model.asm_overlap = asm_overlap 

144 if ilu_levels is not None: 

145 self.sim_model.ilu_levels = ilu_levels 

146 

147 # Set by solve_residual_equations, once, as the LAST write to this 

148 # name for the evaluation: whether the solve that produced the 

149 # accepted u_vec actually converged.  

150 # Mainly used for testing purposes. 

151 self.sim_model.last_solve_converged = False 

152 

153 # Mesh warping only with a shape parameterization (see __init__) 

154 self.shape_parameterization = shape_parameterization 

155 self.shape_parameterized = shape_parameterization is not None 

156 

157 # The mesh warper, built in set_up_sim() 

158 self.mesh_warper = None 

159 

160 # Define data storage object for external post-processing 

161 self.data_store = DataStore() 

162 

163 self.previous_solution = None 

164 

165 # POD reduced-order model (None: full-order solves only) 

166 self.reduced_order_model = reduced_order_model 

167 # (parameter vector, converged) of the last evaluation if the FOM 

168 # solved it, else None: a re-evaluation of that design skips the ROM 

169 # (see _solve) 

170 self._fom_design = None 

171 

172 # Naming suffix for exported files 

173 self.filename_suffix = filename_suffix 

174 

175 # Free the flow solver's matrices after every forward solve and the 

176 # condensed dR/dU after every condensed adjoint 

177 # (sim_model.release_solver), so the two never coexist; they are 

178 # rebuilt when next needed. For memory-bound runs (3D RANS). 

179 self.release_solver_memory = False 

180 

181 # Adjoint reuse. The transposed dR/dU and the preconditioned KSP of the 

182 # state last linearized (condensed Jacobian), kept for the adjoint 

183 # solves of every output of that design: {'key': state key, 'ksp': KSP}; 

184 # freed when the next flow solve starts or a different state is 

185 # linearized. 

186 self._adjoint_cache = None 

187 

188 # Options of the adjoint of a model with a condensed Jacobian (default off; 

189 # set_adjoint_options sets them together with those of the reconstruction): 

190 # adjoint_matrix_free the adjoint operator A_UU^T x + B^T (A_UG^T x), without 

191 # forming A or the product A_UG B 

192 # cache_jacvec_operators the dR/dx operator and dR/dalpha of the 

193 # compute_jacvec_product calls of all outputs of one 

194 # design are built once (freed with the adjoint cache) 

195 self.adjoint_matrix_free = False 

196 self.cache_jacvec_operators = False 

197 # PETSc options (names without the "adj_" prefix) overriding the adjoint KSP's settings of the 

198 # model, e.g. {"ksp_type": "hpddm", "ksp_hpddm_type": "gcrodr", "ksp_hpddm_recycle": 20}; the KSP 

199 # is kept for all outputs of a design, so a recycling method carries its subspace between them 

200 self.adjoint_ksp_options = None 

201 self._dRdx_cache = None 

202 self._dRdalpha_cache = None 

203 

204 # Start each condensed adjoint solve from a first-order solution: the 

205 # solution of the same transposed system with the compact (nearest 

206 # neighbour, gradients frozen) part A_UU of the Jacobian alone -- the 

207 # first-order Jacobian the preconditioner is built from -- solved to a 

208 # relative tolerance first_order_adjoint_rtol. It reuses the 

209 # preconditioner's factors, so costs only Krylov iterations. Applies to 

210 # models with a condensed Jacobian (MUSCL / RANS); the guess is always 

211 # used, and its effect is logged (self._last_first_order_guess). 

212 self.adjoint_first_order_guess = False 

213 self.first_order_adjoint_rtol = 1e-2 

214 self._last_first_order_guess = None 

215 

216 # Write mesh_deformation_<suffix>.bp (the node motions of every design 

217 # evaluation). The solution files are written either way. 

218 self.write_mesh_deformation = True 

219 

220 # Write the output frame (solution and mesh deformation) only for the 

221 # first evaluation of each design: an optimizer that evaluates the 

222 # same design again (CSDL re-runs the model for the derivatives) then 

223 # does not duplicate the frame. Frames keep the eval_idx of that first 

224 # evaluation, so their labels have gaps. 

225 self.write_once_per_design = False 

226 self._last_written_design = None 

227 

228 

229 # Discretization and physics-related parameters 

230 self.poly_order:int = poly_order 

231 self.gamma:float = gamma 

232 self.boundary_condition_dict = boundary_dict 

233 self.moments = moments 

234 

235 # Dolfinx mesh tags for the various boundaries 

236 self.WALL_TAG = 2 # walls (including wing/airfoil boundaries) 

237 self.SYMMETRY_TAG = 5 # symmetry planes (3D: y = 0, or the mesh's symmetry_tags) 

238 self.FARFIELD_TAG = 6 # characteristic far field: the whole outer boundary 

239 

240 # global evaluation index of solve_residual_equations 

241 self.eval_idx = 0 

242 

243 def set_up_sim(self): 

244 # pass the various tolerances and limits to the Euler sim model 

245 self.sim_model.residual_conv_limit = self.euler_residual_conv_limit 

246 self.sim_model.max_newton_iterations = self.max_newton_iterations 

247 

248 # use the free-stream dictionary to define the far-field conditions 

249 self.sim_model.define_farfield_conditions(self.boundary_condition_dict) 

250 

251 # compute boundary facet boolean arrays 

252 inner_facets = self.mesh.inner_bdry_facet_mask 

253 

254 # Find the symmetry face facets (only relevant in 3D at the moment) 

255 # NOTE: The symmetry plane is currently hardcoded at y=0 

256 has_symmetry = self.mesh.symmetry_facet_mask is not None or self.sim_model.dimensions == 3 

257 if self.mesh.symmetry_facet_mask is not None: 

258 symmetry_facets = self.mesh.symmetry_facet_mask # from the mesh's boundary tags 

259 elif self.sim_model.dimensions == 3: 

260 symmetry_facets = self.mesh.bdry_midpoints[:, 1] <= 1e-6 

261 else: 

262 symmetry_facets = np.array([False]*inner_facets.shape[0]) 

263 

264 # Everything that is neither wall nor symmetry plane is far field. The 

265 # characteristic far-field state decides facet by facet whether the flow 

266 # enters or leaves, so no inflow/outflow division (and nothing to 

267 # re-tag when the angle of attack changes) is needed. 

268 farfield_facets = ~np.asarray(inner_facets, dtype=bool) & ~np.asarray(symmetry_facets, dtype=bool) 

269 

270 # Print the local number of facets of each type for each MPI rank 

271 PETSc.Sys.syncPrint("Rank {}, number of wing boundary facets: {}".format(self.mesh.mesh.comm.Get_rank(), np.count_nonzero(inner_facets))) 

272 PETSc.Sys.syncPrint("Rank {}, number of symmetry facets: {}".format(self.mesh.mesh.comm.Get_rank(), np.count_nonzero(symmetry_facets))) 

273 PETSc.Sys.syncPrint("Rank {}, number of far-field facets: {}".format(self.mesh.mesh.comm.Get_rank(), np.count_nonzero(farfield_facets))) 

274 PETSc.Sys.syncFlush() 

275 

276 # Pair the boundary facet sets with the various mesh tags 

277 bdry_meshtags_and_ds_list = [(farfield_facets, self.FARFIELD_TAG), (inner_facets, self.WALL_TAG)] 

278 if has_symmetry: 

279 bdry_meshtags_and_ds_list += [(symmetry_facets, self.SYMMETRY_TAG)] 

280 

281 # Build the boundary facet tags and define the corresponding quadrature measure 

282 self.mesh.tag_boundary_facets(bdry_meshtags_and_ds_list) 

283 self.mesh.form_manager.build_boundary_measure( 

284 measure_metadata={'quadrature_degree': 3*self.poly_order + 1}) 

285 

286 self.mesh.configure_output(self.filename_suffix, 

287 deformation=self.shape_parameterized and self.write_mesh_deformation) 

288 

289 if self.shape_parameterized: 

290 # The JAX-based IDWarp implementation as mesh warper 

291 self.mesh_warper = IDWarp_jax(self.mesh, wall_meshtag_val=self.WALL_TAG, 

292 anchor_meshtag_vals=(self.FARFIELD_TAG,)) 

293 

294 PETSc.Sys.Print("Initialized mesh warping object") 

295 

296 # The shape parameterization gets the wall nodes (replicated, in the 

297 # numbering of the mesh's wall node set) and facets; its 

298 # wall_displacement() is the (M_global, dim) array mesh_warper.evaluate 

299 # takes. 

300 wall_nodes = self.mesh.boundary_node_set((self.WALL_TAG,)) 

301 self.shape_parameterization.setup(self.mesh.mesh.comm, wall_nodes.coords, 

302 wall_nodes.global_idx, 

303 self.mesh.wall_facets(self.WALL_TAG)) 

304 else: 

305 PETSc.Sys.Print("No shape parameterization given: fixed mesh, forward flow analysis") 

306 

307 # define the integration measure for the mesh cells and non-boundary facets 

308 self.mesh.form_manager.build_cell_and_interior_facet_measures( 

309 measure_metadata={'quadrature_degree': 3*self.poly_order + 1}) 

310 

311 # Define the required Dolfinx objects for the simulation, including the initial conditions 

312 self.sim_model.define_elements_functionspaces() 

313 self.sim_model.define_trial_testfunctions() 

314 self.sim_model.compute_initial_conditions_from_farfield() 

315 self.sim_model.interpolate_solution_vector() 

316 

317 # We output the facet mesh tags to an xdmf file, so we can inspect the BCs in paraview 

318 self.mesh.mesh.topology.create_connectivity(self.sim_model.dimensions-1, self.sim_model.dimensions) 

319 with dolfinx.io.XDMFFile(self.mesh.mesh.comm, "output_meshtags.xdmf", "w") as xdmf: 

320 # Write the Mesh object to the file 

321 xdmf.write_mesh(self.mesh.mesh) 

322 # Write the MeshTags object to the file 

323 xdmf.write_meshtags(self.mesh.meshtags, self.mesh.mesh.geometry) 

324 

325 # Define the boundary condition terms of the Euler weak form 

326 self.sim_model.define_farfield_bc(self.FARFIELD_TAG) 

327 self.sim_model.define_wall_bc(self.WALL_TAG) 

328 if has_symmetry: 

329 self.sim_model.define_slipwall_bc(self.SYMMETRY_TAG) 

330 

331 # Construct the Euler weak form for forward evaluation 

332 self.sim_model.compute_weakform() 

333 

334 # Define the file names and objects for solution writes to external files for Paraview 

335 self.sim_model.define_sol_export_files(self.sim_model.functionspaces["V"], 

336 file_name_addendum=self.filename_suffix, 

337 mesh_deformation=self.shape_parameterized) 

338 

339 # We instantiate a separate DG_postprocessor object here so we can print  

340 # forces and moments and coefficients after each forward solve 

341 self.postprocessor = DG_postprocessor( 

342 self.mesh, self.sim_model, self.WALL_TAG, 

343 moments=self.moments, 

344 p_inf_dim=101325.) 

345 

346 if self.reduced_order_model is not None: 

347 self.reduced_order_model.set_up(self.sim_model) 

348 

349 def _log_memory(self, stage): 

350 # We write a log file for each MPI rank that documents the peak local memory usage 

351 # at various points during an optimization run. 

352 mem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss 

353 rank = MPI.COMM_WORLD.Get_rank() 

354 fname = 'memory_logs/r{}_{}.log'.format(rank, self.filename_suffix) 

355 os.makedirs('memory_logs', exist_ok=True) 

356 with open(fname, 'a') as f: 

357 f.write('{} {} {} {} {}\n'.format( 

358 self.eval_idx, datetime.datetime.now(), os.getpid(), mem, stage)) 

359 

360 @property 

361 def ffd(self): 

362 """The shape parameterization's FFD block (None without one, or before 

363 set_up_sim).""" 

364 if self.shape_parameterization is None: 

365 return None 

366 return getattr(self.shape_parameterization, 'ffd', None) 

367 

368 def wall_geometry(self): 

369 """Geometric constraint quantities of the FFD-deformed wall (area/volume, 

370 thickness, planform); the shape parameterization's 

371 shape_parameterization.WallGeometry. 

372 

373 Built on first use. Collective: call it on every rank. 

374 """ 

375 if self.ffd is None: 

376 raise RuntimeError("wall_geometry needs a shape parameterization (an FFD block), " 

377 "set up by set_up_sim") 

378 return self.shape_parameterization.geometry 

379 

380 def deform_mesh(self): 

381 """Mesh node motions (CSDL) for the shape parameterization's current 

382 wall displacement: the first input of evaluate(). Call after the 

383 shape design variables are declared.""" 

384 if not self.shape_parameterized: 

385 raise RuntimeError("deform_mesh needs a shape parameterization") 

386 shape = self.shape_parameterization 

387 return self.mesh_warper.evaluate(shape.wall_displacement(), shape.parameter_vector()) 

388 

389 def evaluate(self, mesh_node_motions:csdl.Variable = None, shape_parameters:csdl.Variable = None, 

390 alpha:csdl.Variable = None): 

391 """ 

392 The flow state u_vec as a CSDL output. 

393 

394 Shape optimization (a shape parameterization was given): 

395 evaluate(mesh_node_motions, shape_parameters, alpha), with 

396 mesh_node_motions from deform_mesh() and shape_parameters the flat 

397 shape parameter vector. Forward analysis on the fixed mesh: 

398 evaluate(alpha=alpha). alpha is always required. 

399 """ 

400 if alpha is None: 

401 raise ValueError("evaluate needs alpha") 

402 if (mesh_node_motions is None) != (shape_parameters is None): 

403 raise ValueError("pass mesh_node_motions and shape_parameters together, or neither") 

404 if mesh_node_motions is not None and not self.shape_parameterized: 

405 raise ValueError("mesh_node_motions given, but there is no shape parameterization / " 

406 "mesh warper (pass shape_parameterization to DG_windtunnel_model)") 

407 

408 # set inputs using self.declare_input 

409 if mesh_node_motions is not None: 

410 self.declare_input('mesh_node_motions', mesh_node_motions) 

411 self.declare_input('shape_parameters', shape_parameters) 

412 self.declare_input('alpha', alpha) 

413 

414 # The output on this rank is only the local portion of the solution vector 

415 u_vec = self.create_output('u_vec', (self.sim_model.u_vec.x.petsc_vec.local_size,)) 

416 

417 if mesh_node_motions is not None: 

418 # We have declared shape_parameters as an input here, just for logging and computation purposes, 

419 # not because we actually need it; mesh_node_motions carries the effect of the shape_parameters. 

420 self.declare_derivative_parameters(of='u_vec', wrt='shape_parameters', dependent=False) 

421 

422 return u_vec 

423 

424 def solve_forward(self, alpha=None): 

425 """ 

426 One flow solve on the current (undeformed) mesh, without CSDL: set 

427 alpha (radians; default: the boundary dictionary's), solve, print the 

428 coefficients. Warm-starts from the previous solve. Call set_up_sim() 

429 first. Returns (u_vec array, converged, forces), forces being 

430 DG_postprocessor.forces_and_coefficients (L, D, c_l, c_d, the moments, and 

431 the friction drag for a viscous model). 

432 """ 

433 if alpha is not None: 

434 self.sim_model.set_angle_of_attack(alpha) 

435 self.postprocessor.set_angle_of_attack(alpha) 

436 self.sim_model.compute_weakform() 

437 solution, converged = self._solve( 

438 self._parameter_vector(None, self.sim_model.farfield_conditions['alpha'])) 

439 self.sim_model.last_solve_converged = converged 

440 if np.isfinite(np.linalg.norm(solution)): 

441 self.previous_solution = solution 

442 return solution, converged, self.postprocessor.forces_and_coefficients(solution) 

443 

444 def solve_residual_equations(self, input_vals, output_vals): 

445 # Define inputs 

446 mesh_node_motions = input_vals['mesh_node_motions'] if 'mesh_node_motions' in input_vals else None 

447 shape_parameters = input_vals['shape_parameters'] if 'shape_parameters' in input_vals else None 

448 alpha = input_vals['alpha'][0] 

449 

450 # Update the angle of attack to its current value 

451 self.sim_model.set_angle_of_attack(alpha) 

452 

453 PETSc.Sys.Print("angle of attack: {} deg".format(np.degrees(alpha))) 

454 PETSc.Sys.Print("shape parameters: {}".format(shape_parameters)) 

455 PETSc.Sys.Print("solve_residual_equations, stepping into mesh.apply_node_motions") 

456 if mesh_node_motions is not None: 

457 PETSc.Sys.Print("mesh_node_motions 2-norm: {}".format(np.linalg.norm(mesh_node_motions))) 

458 

459 # Apply the mesh node motions to the mesh 

460 self.mesh.apply_node_motions(mesh_node_motions, design=shape_parameters) 

461 

462 # If the mesh is not valid (e.g. has inverted cells), set the output to zero and don't run any simulations 

463 if not self.mesh.is_valid: 

464 PETSc.Sys.Print("Mesh is not valid, skipping simulations.") 

465 output_vals['u_vec'] = np.zeros((self.sim_model.u_vec.x.petsc_vec.local_size,)) 

466 self._fom_design = None 

467 if self.reduced_order_model is not None: 

468 self.reduced_order_model.discard_adjoint() 

469 return 

470 

471 # re-define the current weak form if necessary 

472 self.sim_model.compute_weakform() 

473 

474 # --- flow solve (ROM, else FOM) ----------------------------------- 

475 solution_output, solution_converged = self._solve(self._parameter_vector(shape_parameters, alpha)) 

476 

477 self._log_memory("after flow solve, before setting u_vec output") 

478 

479 self.sim_model.last_solve_converged = solution_converged 

480 

481 PETSc.Sys.Print("setting u_vec output...") 

482 output_vals['u_vec'] = solution_output 

483 

484 # Warm start: keep the accepted solution unless it has blown up 

485 if np.isfinite(np.linalg.norm(solution_output)): 

486 PETSc.Sys.Print("Updating warm start solution") 

487 self.previous_solution = solution_output 

488 else: 

489 PETSc.Sys.Print("Not updating warm start solution, using previous value") 

490 

491 design = np.concatenate([[alpha], np.ravel(shape_parameters) 

492 if shape_parameters is not None else []]) 

493 if (self.write_once_per_design and self._last_written_design is not None 

494 and np.array_equal(design, self._last_written_design)): 

495 PETSc.Sys.Print("Same design as the last written frame: not writing eval_idx {}".format( 

496 self.eval_idx)) 

497 else: 

498 PETSc.Sys.Print("Writing aligned output frame at eval_idx {} (FOM solution)".format( 

499 self.eval_idx)) 

500 # we write the deformation and the solution outputs to the same time step indices 

501 if mesh_node_motions is not None: 

502 self.mesh.write_deformation_output(mesh_node_motions, self.eval_idx) 

503 self.sim_model.write_solution_output(self.eval_idx, 

504 solution_array=solution_output) 

505 self._last_written_design = design 

506 

507 # Explicitly collect garbage and cleanup PETSc objects 

508 gc.collect() 

509 PETSc.garbage_cleanup(self.mesh.mesh.comm) 

510 

511 self._log_memory("after Euler PDE solve") 

512 

513 self.eval_idx += 1 

514 

515 @staticmethod 

516 def _parameter_vector(shape_parameters, alpha): 

517 """The ROM snapshot parameter vector of an evaluation: the flat shape 

518 parameters (if any), then alpha.""" 

519 parts = [] if shape_parameters is None else [np.asarray(shape_parameters, dtype=float).ravel()] 

520 return np.concatenate(parts + [np.array([float(alpha)])]) 

521 

522 def _solve(self, parameter_vector): 

523 """ 

524 The accepted flow solution of one evaluation: (owned solution entries, 

525 converged). 

526 

527 Without a reduced-order model this is the full-order (FOM) solve from 

528 the warm start. With one, the ROM is tried first (once it holds enough 

529 snapshots), from the warm start projected onto the basis, and its 

530 solution is accepted when eta = ||R||/||u|| < eta_threshold; the ROM 

531 then holds the adjoint of that solution (rom.adjoint), which 

532 apply_inverse_jacobian uses. Otherwise the FOM solves from the warm 

533 start exactly as without a ROM, and its solution, if converged and 

534 finite, becomes a snapshot at `parameter_vector`. 

535 

536 A re-evaluation of the design the FOM solved last skips the ROM: the 

537 FOM restarts from its own solution (no Newton steps if that 

538 converged), so the design keeps the FOM solution and the full-order 

539 adjoint. CSDL re-runs the model before every derivative evaluation, 

540 and by then the ROM would reproduce that solution, now a snapshot, and 

541 accept it. 

542 """ 

543 # the forward solve needs the memory of the last design's adjoint operators 

544 self._release_adjoint_cache() 

545 rom = self.reduced_order_model 

546 if rom is None: 

547 return self._run_fom_solve() 

548 rom.discard_adjoint() 

549 

550 repeat = self._fom_design is not None and np.array_equal(parameter_vector, self._fom_design[0]) 

551 rom_result = rom_forces = None 

552 if repeat: 

553 PETSc.Sys.Print("Re-evaluating the design the FOM solved last: skipping the ROM") 

554 else: 

555 self.sim_model.interpolate_solution_vector(self.previous_solution) 

556 rom_result = rom.solve(self.sim_model, parameter_vector) 

557 if rom_result is not None: 

558 rom_forces = self._log_rom_solve(rom_result) 

559 accept_rom = rom_result is not None and rom_result.accepted 

560 if accept_rom and not rom.run_fom_every_evaluation: 

561 PETSc.Sys.Print("Accepting ROM solution") 

562 self._fom_design = None 

563 return rom_result.solution, True 

564 

565 solution_fom, fom_converged = self._run_fom_solve() 

566 if rom_forces is not None and fom_converged: 

567 self._log_rom_error(rom_forces) 

568 if accept_rom: 

569 PETSc.Sys.Print("Accepting ROM solution (FOM run for comparison)") 

570 self._fom_design = None 

571 return rom_result.solution, True 

572 

573 PETSc.Sys.Print("Accepting FOM solution") 

574 finite = self.mesh.mesh.comm.allreduce(bool(np.all(np.isfinite(solution_fom))), op=MPI.LAND) 

575 converged = fom_converged and finite 

576 if converged and repeat and self._fom_design[1]: 

577 PETSc.Sys.Print("The snapshot of this design is already stored") 

578 elif converged: 

579 self.sim_model.set_u_vec(solution_fom) 

580 rom.add_snapshot(self.sim_model, parameter_vector) 

581 sigma = rom.snapshot_singular_values() 

582 if self.mesh.mesh.comm.Get_rank() == 0 and sigma is not None: 

583 self.data_store.global_snapshot_sing_vals_per_iteration += [(self.eval_idx, sigma)] 

584 else: 

585 PETSc.Sys.Print("Not adding a ROM snapshot: the FOM solve did not converge") 

586 self._fom_design = (np.array(parameter_vector, dtype=float), converged) 

587 return solution_fom, fom_converged 

588 

589 def _log_rom_solve(self, result): 

590 # coefficients of the ROM solution, logged like the FOM's 

591 _, cl = self.postprocessor.compute_cl(result.solution) 

592 _, cd = self.postprocessor.compute_cd(result.solution) 

593 cms = self._moment_coefficients(result.solution) 

594 if self.mesh.mesh.comm.Get_rank() == 0: 

595 self.data_store.ROM_walltime += [(self.eval_idx, result.walltime)] 

596 self.data_store.rom_relative_residuals += [(self.eval_idx, result.eta)] 

597 self.data_store.ROM_force_coefficients += [(self.eval_idx, cd, cl) + cms] 

598 self.data_store.ROM_meets_threshold += [(self.eval_idx, result.accepted)] 

599 PETSc.Sys.Print("ROM coefficients: c_d: {}, c_l: {}".format(cd, cl) + self._moment_log_string(cms)) 

600 return (cd, cl) + cms 

601 

602 def _log_rom_error(self, rom_forces): 

603 # ROM-vs-FOM coefficient errors against the FOM solve just run 

604 errors = tuple(float(r - f) for r, f in zip(rom_forces, self._last_fom_coefficients)) 

605 if self.mesh.mesh.comm.Get_rank() == 0: 

606 self.data_store.rom_coefficient_errors += [(self.eval_idx,) + errors] 

607 PETSc.Sys.Print("ROM - FOM coefficient errors: c_d: {:.3e}, c_l: {:.3e}".format(*errors[:2]) 

608 + "".join(", c_m_{}: {:.3e}".format(spec.name, e) 

609 for spec, e in zip(self.postprocessor.moment_specs, errors[2:]))) 

610 

611 def _moment_coefficients(self, solution): 

612 # the coefficient of every moment of the postprocessor, in order 

613 return tuple(self.postprocessor.compute_moment(solution, spec)[1] 

614 for spec in self.postprocessor.moment_specs) 

615 

616 def _moment_log_string(self, coefficients): 

617 return "".join(", c_m_{}: {}".format(spec.name, c) 

618 for spec, c in zip(self.postprocessor.moment_specs, coefficients)) 

619 

620 def _run_fom_solve(self): 

621 # Run the FOM model and output its solution 

622 

623 # store pre-FOM time stamp for performance profiling 

624 fom_pre_time = perf_counter() 

625 

626 # use the most recent converged solution as warm start 

627 initial_condition = self.previous_solution 

628 self.sim_model.interpolate_solution_vector(initial_condition) 

629 

630 # solve the FOM 

631 # an incoming state that already meets the convergence criterion 

632 # needs no solver (SNES, KSP/PC, matrices): skip building it 

633 if self.sim_model.incoming_state_converged(): 

634 residual_norm = self.sim_model.residual_vec_physical.norm() 

635 else: 

636 self.sim_model.set_up_solver() 

637 residual_norm = self.sim_model.solve_system().norm() 

638 fom_converged = self.sim_model.solver_converged 

639 eta = self.sim_model.relative_residual_norm(residual_norm) 

640 solution_fom = self.sim_model.u_vec.x.petsc_vec.getArray().copy() 

641 if self.release_solver_memory: 

642 self.sim_model.release_solver() 

643 

644 self._log_memory("after FOM solve, before doing anything else") 

645 

646 # store post-FOM time stamp for performance profiling 

647 fom_post_time = perf_counter() 

648 PETSc.Sys.Print("FOM wall time: {}".format(fom_post_time - fom_pre_time)) 

649 

650 L_fom, cl_fom = self.postprocessor.compute_cl(solution_fom) 

651 D_fom, cd_fom = self.postprocessor.compute_cd(solution_fom) 

652 cms_fom = self._moment_coefficients(solution_fom) 

653 self._last_fom_coefficients = (cd_fom, cl_fom) + cms_fom 

654 

655 if self.mesh.mesh.comm.Get_rank() == 0: 

656 self.data_store.FOM_walltime += [(self.eval_idx, fom_post_time - fom_pre_time)] 

657 self.data_store.fom_relative_residuals += [(self.eval_idx, eta)] 

658 self.data_store.FOM_force_coefficients += [(self.eval_idx, cd_fom, cl_fom) + cms_fom] 

659 PETSc.Sys.Print("FOM coefficients: c_d: {}, c_l: {}".format(cd_fom, cl_fom) + self._moment_log_string(cms_fom)) 

660 PETSc.Sys.Print("FOM relative residual norm eta: {} (||r||: {})".format( 

661 eta, residual_norm)) 

662 

663 return solution_fom, fom_converged 

664 

665 def apply_inverse_jacobian(self, input_vals, outputs, d_outputs, d_residuals, mode): 

666 pre_apply_inverse_jacobian_time = perf_counter() 

667 # Here we solve the adjoint equation when calculating derivatives for optimization 

668 print("starting apply_inverse_jacobian in DG_windtunnel_model...") 

669 

670 mesh_node_motions = input_vals['mesh_node_motions'] if 'mesh_node_motions' in input_vals else None 

671 shape_parameters = input_vals['shape_parameters'] if 'shape_parameters' in input_vals else None 

672 alpha = input_vals['alpha'][0] 

673 u_vec = outputs['u_vec'] # numpy array 

674 

675 if mode == 'rev': 

676 # compute d_residual = (dr_du^-1)*d_outputs 

677 

678 # apply the input mesh node motions to the mesh of the sim model 

679 # and set the angle of attack before constructing the weak form 

680 # and the corresponding solution that we are linearizing around 

681 if mesh_node_motions is not None: 

682 PETSc.Sys.Print("apply_inverse_jacobian, stepping into mesh.apply_node_motions") 

683 self.mesh.apply_node_motions(mesh_node_motions, design=shape_parameters) 

684 self.sim_model.set_angle_of_attack(alpha) 

685 self.sim_model.compute_weakform() 

686 self.sim_model.set_u_vec(u_vec) 

687 

688 rom = self.reduced_order_model 

689 if (rom is not None and rom.adjoint is not None 

690 and rom.adjoint.matches(self._parameter_vector(shape_parameters, alpha), u_vec)): 

691 # The accepted solution is the ROM's: the adjoint of its 

692 # projected (LSPG) equations, lifted to the full state, which 

693 # compute_jacvec_product then uses like the full-order one 

694 d_residuals['u_vec'] = rom.adjoint.solve(self.sim_model, d_outputs['u_vec']) 

695 self._log_adjoint("rom", pre_apply_inverse_jacobian_time) 

696 return 

697 

698 if self.sim_model.has_condensed_jacobian: 

699 self._condensed_adjoint(d_outputs, d_residuals, 

700 self._adjoint_state_key(u_vec, alpha, mesh_node_motions)) 

701 if self.release_solver_memory: 

702 self.sim_model.release_solver() 

703 PETSc.garbage_cleanup() 

704 self._log_adjoint("fom", pre_apply_inverse_jacobian_time) 

705 return 

706 

707 self._drdu_mat = self.sim_model.assemble_dRdu_mat() 

708 if not hasattr(self, '_d_output_petsc'): 

709 self._d_output_petsc = self._drdu_mat.createVecLeft() 

710 self._d_residual_petsc = self._drdu_mat.createVecLeft() 

711 self._F_form = dolfinx.fem.form(self.sim_model.F) 

712 self._residual_vec_petsc = dolfinx.fem.petsc.create_vector( 

713 dolfinx.fem.extract_function_spaces(self._F_form)) 

714 

715 with self._residual_vec_petsc.localForm() as loc: 

716 loc.set(0.0) 

717 dolfinx.fem.petsc.assemble_vector(self._residual_vec_petsc, self._F_form) 

718 self._residual_vec_petsc.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) 

719 PETSc.Sys.Print("max residual entry: {}".format(self._residual_vec_petsc.max()[1])) 

720 residual_norm = self._residual_vec_petsc.norm() 

721 PETSc.Sys.Print("forward eval residual norm: {}".format(residual_norm)) 

722 

723 # (no assembly here any more: assemble_dRdu_mat above already 

724 # returns the matrix assembled. It hands back the same cached 

725 # Mat object every call, so the KSP's operator reference below 

726 # stays valid.) 

727 if not hasattr(self, '_adj_solver'): 

728 self._adj_solver = PETSc.KSP().create(self.mesh.mesh.comm) 

729 self._adj_solver.setOperators(self._drdu_mat) 

730 self._adj_solver.setOptionsPrefix("adj_") 

731 

732 # define the adjoint solver settings 

733 self.sim_model.apply_krylov_solver_settings( 

734 self._adj_solver, max_it=self.linear_solver_max_it, 

735 gmres_restart=200, 

736 monitor_convergence=False, 

737 rtol=self.linear_solver_rtol, 

738 atol=self.linear_solver_atol) 

739 

740 set_petsc_vec_array(self._d_output_petsc, d_outputs['u_vec']) 

741 

742 # solve the transposed linear system and output the solution 

743 self._adj_solver.solveTranspose(self._d_output_petsc, self._d_residual_petsc) 

744 

745 rnorm = self._adj_solver.getResidualNorm() 

746 PETSc.Sys.Print("adj_ solve: {} iterations, residual norm {} [{}]".format( 

747 self._adj_solver.getIterationNumber(), rnorm, 

748 ksp_converged_reason_name(self._adj_solver.getConvergedReason()))) 

749 

750 if rnorm > self.linear_residual_accept_limit: 

751 raise ValueError( 

752 "apply_inverse_jacobian residual norm is too high: {:.6e} > " 

753 "linear_residual_accept_limit {:.1e}".format( 

754 rnorm, self.linear_residual_accept_limit)) 

755 

756 PETSc.Sys.Print("d_residual norm: {}".format(self._d_residual_petsc.norm())) 

757 

758 if mesh_node_motions is not None: 

759 PETSc.Sys.Print("mesh_node_motions min, max: {}, {}".format(mesh_node_motions.min(), mesh_node_motions.max())) 

760 PETSc.Sys.Print("shape parameters: {}".format(input_vals['shape_parameters'])) 

761 PETSc.Sys.Print("alpha: {}".format(np.degrees(alpha))) 

762 

763 d_residuals['u_vec'] = self._d_residual_petsc.getArray() 

764 

765 PETSc.garbage_cleanup() 

766 

767 self._log_adjoint("fom", pre_apply_inverse_jacobian_time) 

768 else: 

769 raise ValueError("Reverse mode required") 

770 

771 def _log_adjoint(self, kind, start_time): 

772 # which adjoint ('rom' or 'fom') the derivative of the last evaluation 

773 # used, and its wall time 

774 walltime = perf_counter() - start_time 

775 if self.mesh.mesh.comm.Get_rank() == 0: 

776 self.data_store.adjoint_solves += [(self.eval_idx - 1, kind, walltime)] 

777 PETSc.Sys.Print("apply_inverse_jacobian total time ({} adjoint): {}".format(kind, walltime))