Coverage for dragonfly_sim/utils/ffd_dv_utils.py: 80%
143 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""Construction of the FFD control-point-motion design variable and its bounds.
3The design variable ``cp_motions`` has shape ``ffd_shape + (n_opt,)``, where the
4trailing axis is a *compressed* index over the optimized spatial directions --
5``cp_motions[..., m]`` moves spatial component ``coord_idxs[m]`` of every control
6point, and directions absent from ``coord_idxs`` stay at their baseline value.
7See shape_parameterization.WallFFD.apply_cp_motions for the application of
8that mapping.
10Everything here is expressed through a single ``dv_spec`` dict, keyed by spatial
11direction index, so the choice of optimized directions and the bounds on each of
12them cannot drift apart:
14 dv_spec = {
15 0: 0.02, # x: symmetric +/- 0.02 everywhere
16 2: (-0.05, 0.10), # z: asymmetric, uniform over control points
17 }
19A value may also be a callable ``f(coords) -> (lower, upper)`` for bounds that
20vary over the block; see spanwise_linear_bounds.
22Sectional shape variables (shape_parameterization.SectionalShape) get their
23arrays from build_sectional_dv, scalar planform variables (sweep, aspect ratio,
24...) from build_planform_dv, and ShapeDVSet collects any mix of these into the
25ordered set of shape design variables (shape_design.FFDShapeParameterization
26declares them through its layers).
28Numpy-only apart from ShapeDVSet's variable creation, which imports csdl
29lazily -- so the bounds logic can be exercised without an MPI environment.
30Callers pass their own ``printer`` (typically ``PETSc.Sys.Print``) to
31print_cp_motion_bounds.
32"""
34import numpy as np
35from typing import NamedTuple, Optional
38class CPMotionDV(NamedTuple):
39 """Everything needed to declare ``cp_motions`` as a design variable.
41 coord_idxs : (n_opt,) int array of optimized spatial components, ascending
42 shape : ffd_shape + (n_opt,), the shape of the design variable
43 value : initial value (zeros -- the motions are deltas from baseline)
44 lower/upper: bounds, in physical units, same shape as value
45 scaler : per-entry scaler for set_as_design_variable, or None
46 """
47 coord_idxs: np.ndarray
48 shape: tuple
49 value: np.ndarray
50 lower: np.ndarray
51 upper: np.ndarray
52 scaler: Optional[np.ndarray]
55def cp_dv_directions(dv_spec):
56 """The optimized spatial components named by dv_spec, sorted ascending.
58 Call this before constructing DG_windtunnel_model, which needs
59 cp_coord_opt_idxs up front; build_cp_motion_dv (which needs the baseline
60 control points, and so can only run after set_up_sim) derives the
61 same array again and returns it for cross-checking.
63 Sorting is load-bearing: it makes the design variable's column order
64 independent of the order the dict literal happens to be written in, and
65 shape_parameterization.validate_cp_coord_opt_idxs enforces the same
66 strictly-increasing invariant.
67 """
68 if not isinstance(dv_spec, dict):
69 raise TypeError("dv_spec must be a dict keyed by spatial direction index, "
70 "got {}".format(type(dv_spec).__name__))
71 if len(dv_spec) == 0:
72 raise ValueError("dv_spec is empty: at least one coordinate direction must "
73 "be a design variable")
75 for direction in dv_spec:
76 if not isinstance(direction, (int, np.integer)) or isinstance(direction, bool):
77 raise TypeError("dv_spec keys must be integer spatial direction indices, "
78 "got key {!r} of type {}".format(direction, type(direction).__name__))
79 if direction < 0:
80 raise ValueError("dv_spec keys must be non-negative, got {}".format(direction))
82 return np.array(sorted(int(d) for d in dv_spec), dtype=int)
85def _resolve_direction_bounds(bound_spec, coords, direction):
86 """One direction's bounds as a pair of (num_control_points,) arrays."""
87 num_cps = coords.shape[0]
89 if callable(bound_spec):
90 lower, upper = bound_spec(coords)
91 elif isinstance(bound_spec, (tuple, list)):
92 if len(bound_spec) != 2:
93 raise ValueError("bounds for direction {} must be a (lower, upper) pair, "
94 "got a sequence of length {}".format(direction, len(bound_spec)))
95 lower, upper = bound_spec
96 elif isinstance(bound_spec, (float, int, np.floating, np.integer)):
97 if bound_spec < 0.:
98 raise ValueError("a scalar bound is the symmetric half-range and must be "
99 "non-negative; direction {} got {}".format(direction, bound_spec))
100 lower, upper = -bound_spec, bound_spec
101 else:
102 raise TypeError("bounds for direction {} must be a scalar, a (lower, upper) pair, "
103 "or a callable(coords) -> (lower, upper), got {}"
104 .format(direction, type(bound_spec).__name__))
106 lower = np.broadcast_to(np.asarray(lower, dtype=np.double), (num_cps,))
107 upper = np.broadcast_to(np.asarray(upper, dtype=np.double), (num_cps,))
109 # cp_motions starts at zero, so the initial point has to lie inside the box.
110 # A pinned direction, (0., 0.), satisfies this.
111 if not np.all(lower <= 0.):
112 raise ValueError("lower bounds for direction {} must be non-positive "
113 "(cp_motions starts at zero), got max {}".format(direction, lower.max()))
114 if not np.all(upper >= 0.):
115 raise ValueError("upper bounds for direction {} must be non-negative "
116 "(cp_motions starts at zero), got min {}".format(direction, upper.min()))
118 return lower, upper
121def build_cp_motion_dv(baseline_controlpoints, dv_spec, scaler_mode='bounds'):
122 """Build the cp_motions design variable arrays from a dv_spec.
124 baseline_controlpoints : (*ffd_shape, ndim) control point coordinates, as
125 held by shape_parameterization.WallFFD.baseline_coefficients
126 dv_spec : {direction: bound_spec}, see the module docstring
127 scaler_mode : 'bounds' derives a scaler of 1/half-range so every direction's
128 box becomes +/-1 in the optimizer's space; None emits no scaler.
130 Returns a CPMotionDV. CSDL flattens design variables and their bounds in C
131 order (csdl_alpha/backends/simulator.py, build_opt_metadata), and this
132 function reshapes in C order too, so no manual index bookkeeping is needed:
133 entry [i, j, k, m] bounds the motion of
134 baseline_controlpoints[i, j, k, coord_idxs[m]].
135 """
136 baseline_controlpoints = np.asarray(baseline_controlpoints, dtype=np.double)
137 if baseline_controlpoints.ndim < 2:
138 raise ValueError("baseline_controlpoints must be (*ffd_shape, ndim), got shape {}"
139 .format(baseline_controlpoints.shape))
141 ffd_shape = baseline_controlpoints.shape[:-1]
142 ndim = baseline_controlpoints.shape[-1]
144 coord_idxs = cp_dv_directions(dv_spec)
145 if np.any(coord_idxs >= ndim):
146 raise ValueError("dv_spec names direction(s) {} but the control points only have "
147 "{} spatial components".format(coord_idxs[coord_idxs >= ndim].tolist(), ndim))
149 coords = baseline_controlpoints.reshape(-1, ndim)
150 num_cps = coords.shape[0]
151 n_opt = coord_idxs.shape[0]
153 lower = np.zeros((num_cps, n_opt), dtype=np.double)
154 upper = np.zeros((num_cps, n_opt), dtype=np.double)
155 for m, direction in enumerate(coord_idxs):
156 lower[:, m], upper[:, m] = _resolve_direction_bounds(dv_spec[int(direction)],
157 coords, int(direction))
159 dv_shape = tuple(ffd_shape) + (n_opt,)
160 lower = lower.reshape(dv_shape)
161 upper = upper.reshape(dv_shape)
163 if scaler_mode is None:
164 scaler = None
165 elif scaler_mode == 'bounds':
166 # +/-1 per direction in the optimizer's space. Without this, directions
167 # whose bounds differ by orders of magnitude give SLSQP's QP badly scaled
168 # columns and the tightly bounded ones barely move -- the same failure the
169 # alpha scaler in the drivers exists to prevent.
170 half_range = np.maximum(np.abs(lower), np.abs(upper))
171 # a pinned direction has zero half-range; leave its scaler at 1
172 scaler = np.divide(1.0, half_range, out=np.ones_like(half_range), where=half_range > 0.)
173 else:
174 raise ValueError("scaler_mode must be 'bounds' or None, got {!r}".format(scaler_mode))
176 return CPMotionDV(coord_idxs=coord_idxs,
177 shape=dv_shape,
178 value=np.zeros(dv_shape, dtype=np.double),
179 lower=lower,
180 upper=upper,
181 scaler=scaler)
184def spanwise_linear_bounds(base_bound, scale_at_root, scale_at_ref, y_ref, y_axis=1):
185 """A bound callable whose magnitude varies linearly with the spanwise coordinate.
187 The bound is base_bound*scale_at_root at y = 0 and base_bound*scale_at_ref at
188 y = y_ref, linearly interpolated in between and clipped to that range outside
189 [0, y_ref] so the magnitude never changes sign. Bounds are symmetric:
190 lower = -bound, upper = +bound.
192 Motivating case: a wing FFD block that tapers from root chord to tip chord, so
193 a single scalar bound is a much larger fraction of the local chord at the tip
194 than at the root. Shrinking the bound along the span keeps the allowed
195 deformation a roughly constant fraction of the local chord.
196 """
197 def bound_function(coords):
198 t = np.clip(coords[:, y_axis]/y_ref, 0., 1.)
199 scale = (1. - t)*scale_at_root + t*scale_at_ref
200 bound = base_bound*scale
201 return -bound, bound
203 return bound_function
206def print_cp_motion_bounds(cp_dv, printer=print, axis_names=('x', 'y', 'z')):
207 """Summarise a CPMotionDV, one line per optimized direction.
209 printer is called once per line; pass PETSc.Sys.Print from an MPI driver to
210 keep the output on rank 0.
211 """
212 printer("cp_motions design variable: shape {}, {} optimized direction(s)"
213 .format(cp_dv.shape, cp_dv.coord_idxs.shape[0]))
214 for m, direction in enumerate(cp_dv.coord_idxs):
215 direction = int(direction)
216 name = axis_names[direction] if direction < len(axis_names) else str(direction)
217 lower_m, upper_m = cp_dv.lower[..., m], cp_dv.upper[..., m]
218 line = (" column {} -> direction {} ({}): lower [{:.4g}, {:.4g}], upper [{:.4g}, {:.4g}]"
219 .format(m, direction, name,
220 lower_m.min(), lower_m.max(), upper_m.min(), upper_m.max()))
221 if cp_dv.scaler is not None:
222 scaler_m = cp_dv.scaler[..., m]
223 line += ", scaler [{:.4g}, {:.4g}]".format(scaler_m.min(), scaler_m.max())
224 printer(line)
227def build_sectional_dv(num_values, bounds, scaler_mode='bounds'):
228 """Arrays for one sectional shape variable (twist, camber, ...).
230 num_values : entries of the variable -- one per section, or fewer for a
231 B-spline profile over the sections (see SectionalShape)
232 bounds : a scalar symmetric half-range, a (lower, upper) pair, or a pair of
233 (num_values,) arrays. As for cp_motions the variable starts at zero (a
234 delta from the baseline), so zero must lie inside the box; pin an entry
235 with (0., 0.).
236 scaler_mode : as in build_cp_motion_dv
238 Returns a CPMotionDV with coord_idxs=None.
239 """
240 shape = (int(num_values),)
241 coords = np.zeros((shape[0], 1))
242 lower, upper = _resolve_direction_bounds(bounds, coords, 'sectional')
243 lower, upper = np.array(lower), np.array(upper)
245 if scaler_mode is None:
246 scaler = None
247 elif scaler_mode == 'bounds':
248 half_range = np.maximum(np.abs(lower), np.abs(upper))
249 scaler = np.divide(1.0, half_range, out=np.ones_like(half_range), where=half_range > 0.)
250 else:
251 raise ValueError("scaler_mode must be 'bounds' or None, got {!r}".format(scaler_mode))
253 return CPMotionDV(coord_idxs=None, shape=shape, value=np.zeros(shape, dtype=np.double),
254 lower=lower, upper=upper, scaler=scaler)
257def build_planform_dv(baseline, bounds, scaler_mode='bounds', name='planform'):
258 """Arrays for one scalar planform variable in absolute units.
260 Unlike cp_motions and the sectional variables, a planform variable holds
261 the quantity itself (an aspect ratio of 8.6, a sweep of 0.44 rad), so it
262 starts at its baseline value and bounds is an absolute (lower, upper)
263 pair that must bracket it.
265 scaler_mode : 'bounds' scales by 1/half-range about the baseline (the
266 larger side), so the box is about +/-1 around the start in the
267 optimizer's space; None emits no scaler.
269 Returns a CPMotionDV of shape (1,) with coord_idxs=None.
270 """
271 baseline = float(baseline)
272 if not isinstance(bounds, (tuple, list, np.ndarray)) or len(bounds) != 2:
273 raise ValueError("bounds for planform variable {!r} must be an absolute "
274 "(lower, upper) pair, got {!r}".format(name, bounds))
275 lower, upper = (float(b) for b in bounds)
276 if not lower <= baseline <= upper:
277 raise ValueError("bounds ({}, {}) for planform variable {!r} must bracket its "
278 "baseline value {} (the variable starts there)"
279 .format(lower, upper, name, baseline))
280 lower, upper = np.array([lower]), np.array([upper])
282 if scaler_mode is None:
283 scaler = None
284 elif scaler_mode == 'bounds':
285 half_range = max(upper[0] - baseline, baseline - lower[0])
286 scaler = np.array([1.0 / half_range if half_range > 0. else 1.0])
287 else:
288 raise ValueError("scaler_mode must be 'bounds' or None, got {!r}".format(scaler_mode))
290 return CPMotionDV(coord_idxs=None, shape=(1,), value=np.array([baseline]),
291 lower=lower, upper=upper, scaler=scaler)
294class ShapeDVSet:
295 """The ordered set of shape design variables of a driver.
297 ``flat_variable`` concatenates them, in declaration order and each variable
298 flattened in C order, into the one vector the model takes as its
299 ``shape_parameters``.
300 """
302 def __init__(self):
303 self.names = []
304 self.dvs = []
305 self.variables = []
307 def add(self, name, dv):
308 """Create the csdl design variable for `dv` (a CPMotionDV) and return it."""
309 import csdl_alpha as csdl
310 if name in self.names:
311 raise ValueError("shape design variable {!r} added twice".format(name))
312 var = csdl.Variable(name=name, value=np.array(dv.value, dtype=np.double))
313 var.set_as_design_variable(lower=dv.lower, upper=dv.upper, scaler=dv.scaler)
314 self.names.append(name)
315 self.dvs.append(dv)
316 self.variables.append(var)
317 return var
319 def __getitem__(self, name):
320 return self.variables[self.names.index(name)]
322 def dv(self, name):
323 return self.dvs[self.names.index(name)]
325 @property
326 def size(self):
327 return int(sum(np.prod(dv.shape) for dv in self.dvs))
329 @property
330 def shape(self):
331 return (self.size,)
333 def flat_variable(self):
334 """All shape variables as one csdl vector, in declaration order.
336 This is what the model takes as its ``shape_parameters``. A single
337 variable is returned as is.
338 """
339 import csdl_alpha as csdl
340 if not self.variables:
341 raise ValueError("no shape design variables added")
342 if len(self.variables) == 1:
343 return self.variables[0]
344 return csdl.concatenate([v.reshape((int(np.prod(v.shape)),)) for v in self.variables],
345 axis=0)