Coverage for dragonfly_sim/core/postprocessor.py: 90%
305 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1import numpy as np
2from mpi4py import MPI
3from petsc4py import PETSc
4import ufl
5import dolfinx
7import csdl_alpha as csdl
9from dragonfly_sim.utils.postprocessor_utils import moment_about
10from dragonfly_sim.utils.Euler_utils import pressure
13class MomentSpec:
14 """One moment of the wall force: about `point`, along unit `direction`
15 (a 3-vector in 2D too, where only the z entry is read), normalised in the
16 coefficient by q*S_ref*ref_length. The values are dolfinx Constants, so the
17 forms built on them stay valid when they change."""
19 def __init__(self, mesh, name, point, direction, ref_length=None):
20 self.name = name
21 self.point = dolfinx.fem.Constant(mesh, point)
22 self.direction = dolfinx.fem.Constant(mesh, direction)
23 # None: the postprocessor's reference chord
24 self.ref_length = None if ref_length is None else dolfinx.fem.Constant(mesh, float(ref_length))
25 self.point_array = point
26 self.direction_array = direction
29class DG_postprocessor(csdl.CustomExplicitOperation):
30 """
31 Forces, and moments about user-given points and directions, of the wall
32 force (pressure, plus the viscous traction for a viscous model), with their
33 derivatives for CSDL.
35 moments: sequence of (name, point, direction[, ref_length]) tuples, each
36 giving the moment M = d . (r x F), r = X - point, of the wall force F along
37 the direction d (right-hand rule); CSDL outputs M_<name> and c_m_<name>.
38 point has one entry per dimension; direction is a nonzero vector,
39 normalised to unit length, or the string "pitch": the nose-up-positive
40 pitching moment, (0, 1, 0) in 3D (y spanwise) and (0, 0, -1) in 2D. In 2D
41 only the out-of-plane direction (0, 0, +-1) exists. c_m_<name> is
42 M / (q S_ref ref_length), the reference chord by default. The directions
43 are fixed in the geometry axes and the points fixed in space, also while
44 a shape design variable moves the wall. Only M_<name> is differentiated:
45 S_ref and the reference chord change with the 3D geometry, so constrain
46 M_<name>, not c_m_<name>. For a half-wing on a symmetry plane y = 0, the
47 root bending moment is ("root_bending", (x0, 0, 0), (1, 0, 0)).
49 Each moment costs its own value, state-derivative and mesh-derivative
50 assembly (a chain-rule evaluation per form for derived-field models).
51 Since M_k = d_k . int(X x f) - (d_k x x0_k) . int(f), all moments are
52 linear combinations of 6 base functionals (3 in 2D), so only 6 derivative
53 evaluations would suffice regardless of the number of moments; this is not
54 yet exploited.
55 """
56 ordered_callbacks = True # as DG_windtunnel_model: collectives and the shared mesh
58 def __init__(self, Mesh_obj, sim_model, wall_meshtag_val, moments=None,
59 ref_area_dim=None, ref_chord_dim=None, testing_mode=False,
60 p_inf_dim=None, L_ref_dim=1.0):
61 super().__init__()
63 self.testing_mode = testing_mode
65 self.mesh_obj = Mesh_obj
66 self.sim_model = sim_model
67 self.wall_meshtag_val = wall_meshtag_val
68 self.moment_specs = self._parse_moments(moments)
70 # Reference quantities
71 self.ref_area_dim = ref_area_dim # reference area (if None, set to 1 in 2D and computed planform area in 3D)
72 self.ref_chord_dim = ref_chord_dim # reference chord length (if None, set to 1)
73 self.p_inf_dim = p_inf_dim # dimensionalized freestream pressure, used for computing physical forces
74 self.L_ref_dim = L_ref_dim # dimensionalized length scale, used for computing physical forces
76 # Compute quantities that are constant over the course of an optimization
77 self.compute_initial_and_constant_quantities()
80 def _parse_moments(self, moments):
81 # Validate the user's (name, point, direction[, ref_length]) tuples and
82 # turn them into MomentSpecs
83 dim = self.sim_model.dimensions
84 specs = []
85 for entry in (moments if moments is not None else ()):
86 if not isinstance(entry, (tuple, list)) or len(entry) not in (3, 4):
87 raise ValueError("a moment is a (name, point, direction[, ref_length]) tuple, got {!r}".format(entry))
88 name, point, direction = entry[:3]
89 ref_length = entry[3] if len(entry) == 4 else None
90 if not isinstance(name, str) or not name.isidentifier():
91 raise ValueError("moment name {!r} is not a valid identifier".format(name))
92 if name in [spec.name for spec in specs]:
93 raise ValueError("duplicate moment name {!r}".format(name))
94 point = np.asarray(point, dtype=np.double)
95 if point.shape != (dim,):
96 raise ValueError("moment {!r}: point needs {} entries, got {}".format(name, dim, point.shape))
97 if isinstance(direction, str):
98 if direction != "pitch":
99 raise ValueError("moment {!r}: unknown direction preset {!r}".format(name, direction))
100 direction = [0., 1., 0.] if dim == 3 else [0., 0., -1.]
101 elif np.ndim(direction) == 0 and dim == 2:
102 direction = [0., 0., direction] # scalar shorthand for the out-of-plane axis
103 direction = np.asarray(direction, dtype=np.double)
104 if direction.shape != (3,):
105 raise ValueError("moment {!r}: direction needs 3 entries, got {}".format(name, direction.shape))
106 norm = np.linalg.norm(direction)
107 if not norm > 0.0:
108 raise ValueError("moment {!r}: direction is zero".format(name))
109 direction = direction / norm
110 if dim == 2 and (direction[0] != 0.0 or direction[1] != 0.0):
111 raise ValueError("moment {!r}: in 2D only the out-of-plane direction (0, 0, +-1) exists".format(name))
112 if ref_length is not None and not float(ref_length) > 0.0:
113 raise ValueError("moment {!r}: ref_length must be positive".format(name))
114 specs.append(MomentSpec(self.mesh_obj.mesh, name, point, direction, ref_length))
115 return specs
117 def compute_initial_and_constant_quantities(self):
118 # This function computes quantities that are unchanged during an optimization,
119 # as well as some quantities that need to be initialized here and are modified later
120 n = self.sim_model.dimensions
121 if self.p_inf_dim is not None and self.L_ref_dim is not None:
122 self.force_scale = float(self.p_inf_dim) * float(self.L_ref_dim) ** (n - 1)
123 self.moment_scale = float(self.p_inf_dim) * float(self.L_ref_dim) ** n
124 else:
125 self.force_scale = 1.0
126 self.moment_scale = 1.0
128 rho_in = dolfinx.fem.Constant(self.mesh_obj.mesh, self.sim_model.farfield_conditions['rho'])
129 u_ref = self.sim_model.farfield_conditions['u_ref']
131 self._S_ref_const = dolfinx.fem.Constant(self.mesh_obj.mesh, self._compute_S_ref())
132 self._chord_ref_const = dolfinx.fem.Constant(self.mesh_obj.mesh, self._compute_chord_ref(self._S_ref_const.value))
134 self.C_infty = 0.5*rho_in*(u_ref**2)*self._S_ref_const
136 # Drag and lift vector directions
137 self.psi_drag = self.sim_model.farfield_conditions['flow_dir']
138 self.psi_lift = dolfinx.fem.Constant(
139 self.mesh_obj.mesh, np.zeros((self.sim_model.dimensions,), dtype=np.double))
141 self.set_angle_of_attack(self.sim_model.farfield_conditions['alpha'])
143 def _compute_S_ref(self):
144 # Compute the reference area or chord length for force, moment coefficient calculation
145 if self.ref_area_dim is not None:
146 # Skip the calculation, use the user-defined value instead
147 return self.ref_area_dim
148 if self.sim_model.dimensions == 2:
149 return 1.0
150 elif self.sim_model.dimensions == 3:
151 # S_ref is equal to 0.5 times the wing planform area, calculated as
152 # the wing surface projected onto the z=0 plane
153 S_ref = 0.5 * self.mesh_obj.mesh.comm.allreduce(
154 dolfinx.fem.assemble_scalar(dolfinx.fem.form(
155 abs(self.mesh_obj.n[2]) * self.mesh_obj.form_manager.ds(self.wall_meshtag_val))), MPI.SUM)
156 PETSc.Sys.Print("Reference (planform) area: {}".format(S_ref))
157 return S_ref
159 def _compute_chord_ref(self, S_ref_value):
160 # Compute the reference chord length for moment coefficient calculation
161 if self.ref_chord_dim is not None:
162 # Skip the calculation, use the user-defined value instead
163 return self.ref_chord_dim
164 if self.sim_model.dimensions == 2:
165 return 1.0
166 elif self.sim_model.dimensions == 3:
167 # Calculate the mean geometric chord of the current geometry
168 bbox = self.mesh_obj.compute_wall_bounding_box(self.wall_meshtag_val)
169 span = bbox[1][1] - bbox[1][0]
170 return S_ref_value / span
172 def refresh_geometry_dependent_reference_quantities(self):
173 # Recompute the reference area and chord length after a geometry update
174 self._S_ref_const.value = self._compute_S_ref()
175 self._chord_ref_const.value = self._compute_chord_ref(self._S_ref_const.value)
177 def set_angle_of_attack(self, alpha):
178 # update the angle of attack, and propagate the update to the lift vector direction
179 # NOTE: The drag vector direction is handled differently, and is propagated automatically
180 self.sim_model.set_angle_of_attack(alpha)
182 inflow_vec = self.sim_model.farfield_conditions['inflow_angle']
183 if self.sim_model.dimensions == 2:
184 lift_unit_vec = np.array([-inflow_vec[1], inflow_vec[0]], dtype=np.double)
185 elif self.sim_model.dimensions == 3:
186 # wing pitches about the spanwise (y) axis, so the lift direction
187 # is the cross product of the inflow and y-direction
188 zero_direction_vec = np.zeros((3,))
189 zero_direction_vec[1] = 1.
190 lift_unit_vec = np.cross(inflow_vec, zero_direction_vec)
192 self.psi_lift.value = lift_unit_vec
194 def _pressure_force_vector(self, u_dof_vec):
195 """The force density on the wall: pressure, plus the viscous traction
196 for a viscous model (sim_model.wall_force_density)."""
197 self.sim_model.set_u_vec(u_dof_vec)
198 if self.sim_model.derived_fields:
199 # the viscous traction reads the reconstructed gradient
200 self.sim_model.update_derived_fields()
201 return self.sim_model.wall_force_density(self.sim_model.u_vec, self.mesh_obj.n)
202 p = pressure(self.sim_model.u_vec, self.sim_model.gamma)
203 return p*self.mesh_obj.n
205 def _assemble_value_and_coefficient(self, form, normalizer):
206 raw_val = self.mesh_obj.mesh.comm.allreduce(
207 dolfinx.fem.assemble_scalar(dolfinx.fem.form(form)), MPI.SUM)
208 coeff_val = self.mesh_obj.mesh.comm.allreduce(
209 dolfinx.fem.assemble_scalar(dolfinx.fem.form((1.0/normalizer)*form)), MPI.SUM)
210 return raw_val, coeff_val
212 def _assemble_state_derivative(self, form, label):
213 if self.sim_model.derived_fields:
214 # total derivative through the reconstructed gradient
215 vec = self.sim_model.state_derivative(form)
216 PETSc.Sys.Print("{} norm: {}".format(label, vec.norm()))
217 out = vec.getArray().copy()
218 vec.destroy()
219 return out
220 d_duvec = ufl.derivative(form, self.sim_model.u_vec)
221 d_duvec_vec = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(d_duvec))
222 d_duvec_vec.assemble()
224 PETSc.Sys.Print("{} norm: {}".format(label, d_duvec_vec.norm()))
226 output_vec = d_duvec_vec.getArray().copy()
227 d_duvec_vec.destroy()
228 return output_vec
230 def _assemble_mesh_derivative(self, form, label):
231 if self.sim_model.derived_fields:
232 # total derivative through the derived fields (gradient, cell
233 # centres, wall distance)
234 dK = self.sim_model.mesh_derivative(form)
235 PETSc.Sys.Print("{} norm: {}".format(label, dK.norm()))
236 d_dmesh = dK.getArray().copy()
237 dK.destroy()
238 return (self.mesh_obj.dof_to_node_map.T @ d_dmesh.reshape(-1, self.mesh_obj.mesh.geometry.dim)).flatten()
240 spatial_coordinate = ufl.SpatialCoordinate(self.mesh_obj.mesh)
242 args = form.arguments()
243 n = max(a.number() for a in args) if args else -1
244 dX = ufl.Argument(self.mesh_obj.coordinate_space, n + 1)
246 a_deriv_form = dolfinx.fem.form(ufl.derivative(form, spatial_coordinate, dX))
248 dK = dolfinx.fem.petsc.assemble_vector(a_deriv_form)
249 dK.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
251 PETSc.Sys.Print("{} norm: {}".format(label, dK.norm()))
253 d_dmesh = dK.getArray().copy()
254 dK.destroy()
256 d_dmesh_shape_1_dim = self.mesh_obj.dof_to_node_map.T@d_dmesh.reshape(-1, self.mesh_obj.mesh.geometry.dim)
257 return d_dmesh_shape_1_dim.flatten()
259 def _assemble_alpha_derivative(self, form):
260 """d(form)/dalpha at fixed U through the derived fields, e.g. a
261 reconstructed gradient whose boundary facets carry the alpha-dependent
262 far-field state; zero without derived fields. The force directions are
263 Constants (set_angle_of_attack), so their rotation is not part of it."""
264 if not self.sim_model.derived_fields:
265 return 0.0
266 return self.sim_model.alpha_derivative(form)
268 def compute_force_in_direction(self, u_dof_vec, direction_vec):
269 force_density = self._pressure_force_vector(u_dof_vec)
270 return ufl.dot(direction_vec, force_density)*self.mesh_obj.form_manager.ds(self.wall_meshtag_val)
272 def compute_cl(self, u_dof_vec):
273 lift = self.compute_force_in_direction(u_dof_vec, self.psi_lift)
274 return self._assemble_value_and_coefficient(lift, self.C_infty)
276 def compute_L_deriv(self, u_dof_vec):
277 lift = self.compute_force_in_direction(u_dof_vec, self.psi_lift)
278 return self._assemble_state_derivative(lift, "dL_duvec")
280 def compute_L_meshderiv(self, u_dof_vec):
281 lift = self.compute_force_in_direction(u_dof_vec, self.psi_lift)
282 return self._assemble_mesh_derivative(lift, "dL_dmesh")
284 def compute_cd(self, u_dof_vec):
285 drag = self.compute_force_in_direction(u_dof_vec, self.psi_drag)
286 return self._assemble_value_and_coefficient(drag, self.C_infty)
288 def compute_D_deriv(self, u_dof_vec):
289 drag = self.compute_force_in_direction(u_dof_vec, self.psi_drag)
290 return self._assemble_state_derivative(drag, "dD_duvec")
292 def compute_D_meshderiv(self, u_dof_vec):
293 drag = self.compute_force_in_direction(u_dof_vec, self.psi_drag)
294 return self._assemble_mesh_derivative(drag, "dD_dmesh")
296 def _moment_normalizer(self, spec):
297 ref_length = self._chord_ref_const if spec.ref_length is None else spec.ref_length
298 return self.C_infty * ref_length
300 def compute_moment_form(self, u_dof_vec, spec):
301 """The moment M = d . (r x F) of the wall force about spec.point along
302 spec.direction, as a scalar form."""
303 force_density = self._pressure_force_vector(u_dof_vec)
304 coords = ufl.SpatialCoordinate(self.mesh_obj.mesh)
305 return moment_about(force_density, coords - spec.point, spec.direction,
306 self.sim_model.dimensions) * self.mesh_obj.form_manager.ds(self.wall_meshtag_val)
308 def compute_moment(self, u_dof_vec, spec):
309 """(M, c_m) of one moment."""
310 moment = self.compute_moment_form(u_dof_vec, spec)
311 return self._assemble_value_and_coefficient(moment, self._moment_normalizer(spec))
313 def compute_moment_deriv(self, u_dof_vec, spec):
314 moment = self.compute_moment_form(u_dof_vec, spec)
315 return self._assemble_state_derivative(moment, "dM_{}_duvec".format(spec.name))
317 def compute_moment_meshderiv(self, u_dof_vec, spec):
318 moment = self.compute_moment_form(u_dof_vec, spec)
319 return self._assemble_mesh_derivative(moment, "dM_{}_dmesh".format(spec.name))
321 def compute_moment_alphaderiv(self, u_dof_vec, spec):
322 # explicit alpha dependence at fixed state: the direction is fixed in the geometry axes,
323 # so only the derived fields carry any
324 return self._assemble_alpha_derivative(self.compute_moment_form(u_dof_vec, spec))
326 def forces_and_coefficients(self, u_dof_vec=None, alpha=None):
327 """
328 Lift and drag and their coefficients of a flow state, and the moments
329 and their coefficients, outside any CSDL graph: no FFD block, mesh
330 motion or recorder needed.
332 u_dof_vec: the owned state entries (default: the model's current
333 u_vec); alpha: angle of attack in radians (default: the model's). For a
334 viscous model, D_friction / c_d_friction is the part of the drag from
335 the viscous wall traction. Collective. Returns a dict; "moments" maps
336 each moment name to dict(M=..., c_m=...).
337 """
338 if u_dof_vec is None:
339 u_dof_vec = self.sim_model.u_vec.x.petsc_vec.getArray().copy()
340 if alpha is not None:
341 self.set_angle_of_attack(alpha)
342 self.refresh_geometry_dependent_reference_quantities()
343 L, c_l = self.compute_cl(u_dof_vec)
344 D, c_d = self.compute_cd(u_dof_vec)
345 out = dict(L=L, D=D, c_l=c_l, c_d=c_d)
346 out["moments"] = {}
347 for spec in self.moment_specs:
348 M, c_m = self.compute_moment(u_dof_vec, spec)
349 out["moments"][spec.name] = dict(M=M, c_m=c_m)
350 if getattr(self.sim_model, "wall_traction", None) is not None:
351 # drag of the viscous traction alone
352 self.sim_model.set_u_vec(u_dof_vec)
353 self.sim_model.update_derived_fields()
354 t = self.sim_model.wall_traction(self.sim_model.u_vec, self.mesh_obj.n)
355 form = ufl.dot(self.psi_drag, t) * self.mesh_obj.form_manager.ds(self.wall_meshtag_val)
356 out["D_friction"], out["c_d_friction"] = self._assemble_value_and_coefficient(form, self.C_infty)
357 return out
359 def evaluate(self, u_vec, mesh_deformation=None, alpha=None, shape_param_inputs=None):
360 """Forces and coefficients of the state u_vec. mesh_deformation is the
361 volume mesh motion of a shape optimization; leave it out (None) for a
362 forward analysis on the fixed mesh. alpha is required."""
363 if alpha is None:
364 raise ValueError("evaluate needs alpha")
365 # declare input quantities
366 self.declare_input('u_vec', u_vec)
367 if mesh_deformation is not None:
368 self.declare_input('mesh_deformation', mesh_deformation)
369 self.declare_input('alpha', alpha)
370 if shape_param_inputs is not None:
371 self.declare_input('shape_param_inputs', shape_param_inputs)
373 # declare output variables
374 c_l = self.create_output('c_l', (1,))
375 c_d = self.create_output('c_d', (1,))
376 L = self.create_output('L', (1,))
377 D = self.create_output('D', (1,))
378 M = {spec.name: self.create_output('M_' + spec.name, (1,)) for spec in self.moment_specs}
379 c_m = {spec.name: self.create_output('c_m_' + spec.name, (1,)) for spec in self.moment_specs}
381 # Force and moment coefficients are only used as monitoring outputs;
382 # to avoid CSDL from computing its derivatives we set these to zero
383 input_names = ('u_vec', 'alpha') + (('mesh_deformation',) if mesh_deformation is not None else ())
384 coefficient_names = ('c_l', 'c_d') + tuple('c_m_' + spec.name for spec in self.moment_specs)
385 for _input_name in input_names:
386 for _output_name in coefficient_names:
387 self.declare_derivative_parameters(of=_output_name, wrt=_input_name, dependent=False)
388 if shape_param_inputs is not None:
389 for _output_name in coefficient_names + ('L', 'D') + tuple('M_' + spec.name for spec in self.moment_specs):
390 self.declare_derivative_parameters(of=_output_name, wrt='shape_param_inputs',
391 dependent=False)
393 # construct output of the model
394 output = csdl.VariableGroup()
395 output.c_l = c_l
396 output.c_d = c_d
397 output.L = L
398 output.D = D
399 for spec in self.moment_specs:
400 setattr(output, 'M_' + spec.name, M[spec.name])
401 setattr(output, 'c_m_' + spec.name, c_m[spec.name])
403 return output
405 def compute(self, input_vals, output_vals):
406 PETSc.Sys.Print("starting compute in DG_postprocessor...")
408 u_vec = input_vals['u_vec']
409 alpha = input_vals['alpha'][0]
410 deform_array = input_vals['mesh_deformation'] if 'mesh_deformation' in input_vals else None
411 shape_parameters = input_vals['shape_param_inputs'] if 'shape_param_inputs' in input_vals else None
413 # Update the angle of attack and downstream quantities
414 self.set_angle_of_attack(alpha)
416 # Apply the mesh deformations to the mesh (none in a forward analysis)
417 if deform_array is not None:
418 self.mesh_obj.apply_node_motions(deform_array, design=shape_parameters)
420 # check whether the mesh and solution on it are valid
421 if (self.mesh_obj.is_valid and self.sim_model.last_solve_converged and not np.isnan(np.linalg.norm(u_vec))) or self.testing_mode:
422 # update geometry-dependent quantities, under the influence of the deformed mesh
423 self.refresh_geometry_dependent_reference_quantities()
425 # compute all output quantities and print them
426 lift_val, c_l = self.compute_cl(u_vec)
427 output_vals['L'] = lift_val
428 output_vals['c_l'] = c_l
429 drag_val, c_d_initial = self.compute_cd(u_vec)
430 output_vals['D'] = drag_val
431 output_vals['c_d'] = c_d_initial #+ max(0., 10.*(0.285 - c_l))
433 PETSc.Sys.Print("c_d: {}, c_l: {}".format(np.round(c_d_initial, 6), np.round(output_vals['c_l'], 6)))
434 PETSc.Sys.Print("D: {}, L: {} (nondimensional)".format(np.round(output_vals['D'], 6), np.round(output_vals['L'], 6)))
436 # print dimensional forces, just for monitoring purposes
437 PETSc.Sys.Print("D: {}, L: {} (dimensional)".format(
438 np.round(output_vals['D'] * self.force_scale, 6),
439 np.round(output_vals['L'] * self.force_scale, 6)))
441 for spec in self.moment_specs:
442 moment_val, output_vals['c_m_' + spec.name] = self.compute_moment(u_vec, spec)
443 output_vals['M_' + spec.name] = moment_val
444 PETSc.Sys.Print("moment {} (point {}, direction {}): M: {} (nondimensional), c_m: {}, M: {} (dimensional)".format(
445 spec.name, np.round(spec.point_array, 6), np.round(spec.direction_array, 6),
446 np.round(moment_val, 6), np.round(output_vals['c_m_' + spec.name], 6),
447 np.round(moment_val * self.moment_scale, 6)))
448 else:
449 PETSc.Sys.Print("Mesh or solution is not valid, setting outputs to high drag.")
450 output_vals['L'], output_vals['c_l'] = 0.0, 0.0
451 output_vals['D'], output_vals['c_d'] = 10., 1.0
452 for spec in self.moment_specs:
453 output_vals['M_' + spec.name], output_vals['c_m_' + spec.name] = 0.0, 0.0
455 def compute_derivatives(self, input_vals, outputs_vals, derivatives):
456 # Compute the derivatives of the output quantities w.r.t. the inputs
457 PETSc.Sys.Print("starting compute_derivatives in DG_postprocessor...")
459 u_vec = input_vals['u_vec']
460 alpha = input_vals['alpha'][0]
461 deform_array = input_vals['mesh_deformation'] if 'mesh_deformation' in input_vals else None
462 shape_parameters = input_vals['shape_param_inputs'] if 'shape_param_inputs' in input_vals else None
464 # Update the angle of attack and downstream quantities
465 self.set_angle_of_attack(alpha)
467 # Apply the mesh deformations to the mesh and update geometry-dependent quantities
468 if deform_array is not None:
469 self.mesh_obj.apply_node_motions(deform_array, design=shape_parameters)
470 self.refresh_geometry_dependent_reference_quantities()
472 # Compute all force and moment derivatives w.r.t. the nondimensional state and mesh vertex movement
473 derivatives['L', 'u_vec'] = self.compute_L_deriv(u_vec)[None, :]
474 derivatives['D', 'u_vec'] = self.compute_D_deriv(u_vec)[None, :]
475 for spec in self.moment_specs:
476 derivatives['M_' + spec.name, 'u_vec'] = self.compute_moment_deriv(u_vec, spec)[None, :]
478 if deform_array is not None:
479 derivatives['L', 'mesh_deformation'] = self.compute_L_meshderiv(u_vec)[None, :]
480 derivatives['D', 'mesh_deformation'] = self.compute_D_meshderiv(u_vec)[None, :]
481 for spec in self.moment_specs:
482 derivatives['M_' + spec.name, 'mesh_deformation'] = self.compute_moment_meshderiv(u_vec, spec)[None, :]
484 # Compute the force derivatives w.r.t. angle of attack changes at fixed u_vec: the
485 # rotation of the lift and drag directions (dL = -D dalpha, dD = L dalpha), plus
486 # the forces' dependence on alpha through the derived fields
487 L_val, _ = self.compute_cl(u_vec)
488 D_val, _ = self.compute_cd(u_vec)
489 dL = self._assemble_alpha_derivative(self.compute_force_in_direction(u_vec, self.psi_lift))
490 dD = self._assemble_alpha_derivative(self.compute_force_in_direction(u_vec, self.psi_drag))
491 derivatives['L', 'alpha'] = np.array([[-D_val + dL]], dtype=np.float64)
492 derivatives['D', 'alpha'] = np.array([[ L_val + dD]], dtype=np.float64)
493 for spec in self.moment_specs:
494 derivatives['M_' + spec.name, 'alpha'] = np.array(
495 [[self.compute_moment_alphaderiv(u_vec, spec)]], dtype=np.float64)