User guide

A DRAGonFLy case is a CSDL model built from three operations: the mesh warper, the flow model (DG_windtunnel_model) and the force postprocessor (DG_postprocessor). This page walks through setting one up, in the order of the scripts in examples/. Run the script with MPI (OMP_NUM_THREADS=1 mpirun -n <ranks> python script.py).

1. Mesh

Read a DOLFINx mesh with load_dolfinx_mesh, which handles XDMF, gmsh .msh and structured CGNS files and keeps the boundary tags the file carries:

from dragonfly_sim.utils.mesh_io_utils import load_dolfinx_mesh

mesh, info = load_dolfinx_mesh("naca0012_euler_mesh_quad_v2.xdmf")

A mesh read with plain DOLFINx (xdmf.read_mesh()) works too, with the boundary function.

The solver needs to know which boundary facets are the wall (the body) and which are symmetry planes; all others are far field. Either pass a function of the facet midpoints (mesh_inner_bdry_function, with airfoil_inner_bdry_function and wing_inner_bdry_function in dragonfly_sim.utils.mesh_manager_utils for the example meshes), or name the mesh’s own boundary tags (wall_tags, symmetry_tags). The function takes precedence and, if the mesh also has tags, a warning lists them. Giving the shape parameterization a block_corner_list is a third way to find the wall. The formats, both options and the checks are described in Meshes.

Requirements:

Orientation: the freestream flows in \(+x\). In 2D, y is vertical; in 3D, y is spanwise and z vertical.
3D: model a half-wing whose root lies on the symmetry plane \(y = 0\).
Far field: the outer boundary should be far from the body. All of it that is neither wall nor symmetry plane gets the characteristic far-field condition.

2. Flow conditions

boundary_dict = {'rho': 1.0, 'M': 0.8, 'p': 1.0, 'alpha': np.radians(2)}

This free-stream state is all that is specified: density, pressure, Mach number and angle of attack (radians). The whole outer boundary is a characteristic (Riemann-invariant) far field built from it, which takes each facet’s inflow or outflow character from the flow itself, so the angle of attack can be changed by any amount, for example in an optimization over a wide range of angles. The far field is subsonic: \(M < 1\).

The former form {'inlet': {...}, 'outlet': {...}} is no longer accepted; pass the contents of 'inlet' and drop 'outlet'.

3. Flow model

from dragonfly_sim.core.windtunnel_model import DG_windtunnel_model
from dragonfly_sim.core.shape_design import FFDShapeParameterization, ControlPointMotions
from dragonfly_sim.utils.mesh_manager_utils import airfoil_inner_bdry_function

# 5 x 3 FFD control points; they move vertically, within +/-0.01
shape = FFDShapeParameterization([5, 3], [2, 2], layers=[ControlPointMotions({1: 0.01})])
model = DG_windtunnel_model(mesh, boundary_dict, shape, moments=[("pitch", (0.25, 0.), "pitch")],
                            mesh_inner_bdry_function=airfoil_inner_bdry_function,
                            poly_order=0, filename_suffix="my_case", asm_overlap=1, ilu_levels=1)
model.set_up_sim()

poly_order=0 is required in this release. The shape parameterization holds the FFD block and the design variables with their bounds (see Shape parameterization); the model does not depend on which variables they are. set_up_sim() tags the boundaries, builds the weak form and the mesh warper (model.mesh_warper), and sets up the shape parameterization’s FFD block (model.ffd).

Solver settings are attributes; change them before set_up_sim():

Attribute	Default	Meaning
`euler_residual_conv_limit`	`1e-9`	residual norm at which the flow solve has converged
`max_newton_iterations`	`100`	total Newton step budget per flow solve
`linear_solver_rtol`, `linear_solver_atol`	`1e-20`, `1e-11`	adjoint and tangent linear solves
`linear_solver_max_it`	`1000`	iteration limit of those solves

Constructor options select the flow model:

Argument	Default	Meaning
`model_class`, `model_kwargs`	`CompressibleEulerModel`, `{}`	e.g. `CompressibleRANSModel`, `{'Re': 6e6}` for SA-neg RANS (see RANS, time integration and far fields); the body becomes an adiabatic no-slip wall
`mesh_inner_bdry_function`, `wall_tags`, `symmetry_tags`	`None`	how the wall and symmetry planes are found (see Mesh above)
`reduced_order_model`	`None`	a `ReducedOrderModel` to try before every flow solve (see Reduced-order modeling)

RANS solves use pseudo-transient continuation, so max_newton_iterations counts pseudo steps (~60 from the free stream on the coarse NACA0012 C-grid at \(M = 0.7\)); raise it to a few hundred.

4. Postprocessor

from dragonfly_sim.core.postprocessor import DG_postprocessor

post = DG_postprocessor(model.mesh, model.sim_model, model.WALL_TAG, np.array([0.25, 0.]), p_inf_dim=101325.)

The aerodynamic centre argument may be left out (quarter chord). ref_area_dim and ref_chord_dim override the reference area and chord of the coefficients. p_inf_dim and L_ref_dim only scale the dimensional forces printed to the console.

5. Design variables and model graph

Build the CSDL graph between recorder.start() and recorder.stop():

recorder = csdl.Recorder(inline=False)
recorder.start()
# ... construct model and post (steps 3 and 4) ...

params = shape.declare_design_variables()   # every shape variable, with its bounds

alpha = csdl.Variable(name='alpha', value=np.array([np.radians(2)]))
alpha.set_as_design_variable(lower=np.radians(1.75), upper=np.radians(2.25), scaler=1/np.radians(1))

mesh_motion = model.deform_mesh()
u = model.evaluate(mesh_motion, params, alpha)
out = post.evaluate(u, mesh_motion, alpha, params)

out.D.set_as_objective()
out.L.set_as_constraint(lower=0.225)
recorder.stop()

post.evaluate returns L, D and, for every moment (see Moments), M_<name>, which have derivatives and may be used as objectives and constraints, and the coefficients c_l, c_d and c_m_<name>, which are for monitoring only. The flattened shape parameter vector params is passed along for bookkeeping; it carries no derivatives.

Moments

The moments of the wall force to compute, print after every evaluation and differentiate are given to DG_windtunnel_model (or DG_postprocessor) as moments, a list of (name, point, direction[, ref_length]) tuples; there are none by default. Each is \(M = \mathbf{d}\cdot(\mathbf{r}\times\mathbf{F})\) with \(\mathbf{r} = \mathbf{x} - \) point, integrated over the wall (right-hand rule). direction is any nonzero vector (normalised), or "pitch": the nose-up-positive pitching moment, \((0, 1, 0)\) in 3D (spanwise \(y\)) and \((0, 0, -1)\) in 2D, where only the out-of-plane direction \((0, 0, \pm 1)\) exists. The root bending moment of a half-wing on the \(y = 0\) symmetry plane, and its pitching moment, both about the root quarter chord:

moments = [("root_bending", (0.25, 0., 0.), (1., 0., 0.)),
           ("pitch",        (0.25, 0., 0.), "pitch")]
model = DG_windtunnel_model(mesh, boundary_dict, shape, moments=moments, ...)
post = DG_postprocessor(model.mesh, model.sim_model, model.WALL_TAG, moments=moments, p_inf_dim=101325.)
out = post.evaluate(u, mesh_motion, alpha, params)
out.M_root_bending.set_as_constraint(upper=...)

The outputs are M_<name> (differentiated with respect to the state, the mesh motion and \(\alpha\)) and c_m_<name> (monitoring only). The coefficient divides by \(\tfrac{1}{2}\rho_\infty|\mathbf{u}_\infty|^2 S_{ref}\,L_k\), where \(L_k\) is the reference chord unless the fourth tuple entry gives another length, e.g. the semi-span for a bending coefficient. The reference area and chord change with the 3D geometry, so constrain M_<name>, not c_m_<name>. The directions are fixed in the geometry axes and the points fixed in space, also while a shape variable moves the wall: a variable that moves the root station does not move the point. The names must be unique identifiers. Each moment costs one assembly of its value and of each derivative.

6. Running

sim = csdl.experimental.JaxSimulator(recorder, gpu=False)
sim.run()                          # one analysis
sim.check_totals(step_size=1e-5)   # compare derivatives with finite differences

from modopt import CSDLAlphaProblem, SLSQP
prob = CSDLAlphaProblem(problem_name='shape_opt', simulator=sim)
SLSQP(prob, solver_options={'ftol': 1e-8, 'maxiter': 100}).solve()

Each flow solve starts from the previous converged solution.

Forward analysis (no shape variables)

Leave out the shape parameterization: set_up_sim() then builds no FFD block and no mesh warper (model.ffd and model.mesh_warper stay None), the mesh stays fixed, and the angle of attack is the only input.

model = DG_windtunnel_model(mesh, boundary_dict, moments=[("pitch", (0.25, 0.), "pitch")],
                            mesh_inner_bdry_function=airfoil_inner_bdry_function, poly_order=0)
model.set_up_sim()

# one flow solve, no CSDL graph: the state, a convergence flag, and the forces
u, converged, forces = model.solve_forward(np.radians(2))
print(forces["c_l"], forces["c_d"], forces["moments"]["pitch"]["c_m"])

# or a CSDL graph with alpha as its only input
u = model.evaluate(alpha=alpha)
out = model.postprocessor.evaluate(u, alpha=alpha)

DG_postprocessor needs no FFD definitions either: post.forces_and_coefficients(u, alpha) returns L, D, c_l, c_d, the moments (a dict of {"M": ..., "c_m": ...} per name) and, for a viscous model, the friction drag D_friction, c_d_friction of any state, outside CSDL, and post.evaluate(u, alpha=alpha) leaves out the mesh deformation. See examples/airfoil_analysis.py. Time-dependent runs are set up directly from the model classes, as in examples/flow_analysis.py.

Reduced-order model

Pass a ReducedOrderModel to try a POD reduced-order model before every flow solve, in solve_forward and in the CSDL graph alike:

from dragonfly_sim.core.reduced_order_model import ReducedOrderModel

rom = ReducedOrderModel(rb_size=20, eta_threshold=1e-10)
model = DG_windtunnel_model(mesh, boundary_dict, ..., reduced_order_model=rom)

Once two converged solutions exist, each evaluation first solves the least-squares Petrov-Galerkin problem in their POD basis. Its solution is used if \(\eta = \|R\|/\|u\| <\) eta_threshold; otherwise the full-order solve runs as usual and its solution is added to the basis. run_fom_every_evaluation=True runs the full-order solve every time, to compare the two. Basis, inner product and stopping rules are described in Reduced-order modeling.

Output files

Written to the working directory, named with filename_suffix:

FOM_solution_*.bp, mesh_deformation_*.bp: solution and mesh deformation per evaluation, in ADIOS2/VTX format (open in ParaView). The solution file holds separate fields: density, momentum (vector), energy (\(\rho E\)), each transported turbulence variable (e.g. rho_nu_tilde), velocity (vector) and pressure. With a shape parameterization, the solution is written on the deformed mesh it was solved on, and the mesh deformation on the baseline mesh (apply it with Warp By Vector); frame \(k\) of each file belongs to the same evaluation. write_mesh_deformation = False leaves out the mesh deformation file, and write_once_per_design = True writes one frame per design when the optimizer evaluates a design twice;
output_meshtags.xdmf: the boundary tags, to check the far-field/wall/symmetry tags;
memory_logs/: peak memory per MPI rank.

model.data_store collects per-evaluation results (with a reduced-order model also its wall times, \(\eta\), coefficients, acceptance, whether the LSPG solver converged, errors against the full-order solve, and the snapshot singular values: the ROM_* and rom_* entries); model.data_store.write_store_to_numpy_file(save_folder, save_filename) saves them to an existing folder.

To write other fields, use dragonfly_sim.utils.filewriter.FileWriter(filename, comm, function_space, mesh_deformation=False, time_dependent=False) and call interpolate_and_write(f, write_counter=...), or interpolate_and_write(f, time=t) for a time-dependent writer. With mesh_deformation=True each frame is written on the mesh geometry at the time of the write. Piecewise-constant and continuous fields are gathered to rank 0 and written as one block, so ParaView shows no cracks at the MPI partition boundaries.