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, typically from XDMF:

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "naca0012_euler_mesh_quad_v2.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh()

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; it is split automatically into inflow and outflow.

  • Wall: identified by a function mesh_inner_bdry_function(x) that takes facet midpoints of shape (dim, n) and returns a boolean mask. dragonfly_sim.utils.mesh_manager_utils provides airfoil_inner_bdry_function and wing_inner_bdry_function for the meshes used in the examples. Alternatively, pass ffd_block_corner_list, and the wall facets inside that FFD block are used.

2. Flow conditions

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

Only subsonic inflow and outflow are supported. The inflow/outflow split of the far field is frozen at the initial alpha; changing the angle of attack by more than 1 degree from it raises an error, so keep the design variable bounds within that range.

3. Flow model

from dragonfly_sim.core.windtunnel_model import DG_windtunnel_model
from dragonfly_sim.utils.mesh_manager_utils import airfoil_inner_bdry_function
from dragonfly_sim.utils.ffd_dv_utils import cp_dv_directions

dv_spec = {1: 0.01}   # control points move vertically, within +/-0.01
model = DG_windtunnel_model(mesh, boundary_dict, ffd_shape=[5, 3], aero_center=np.array([0.25, 0.]),
                            mesh_inner_bdry_function=airfoil_inner_bdry_function,
                            cp_coord_opt_idxs=cp_dv_directions(dv_spec),
                            poly_order=0, ffd_degree=[2, 2],
                            filename_suffix="my_case", asm_overlap=1, ilu_levels=1)
model.set_up_sim()

poly_order=0 is required in this release. cp_coord_opt_idxs lists the spatial directions in which control points may move (see Shape parameterization). set_up_sim() tags the boundaries, builds the weak form, the FFD block (model.ffd) and the mesh warper (model.mesh_warper).

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

log_cm_alpha

True

print \(dc_m/d\alpha\) after each solve (one extra linear solve)

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.)

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():

from dragonfly_sim.utils.ffd_dv_utils import ShapeDVSet, build_cp_motion_dv

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

shape_dvs = ShapeDVSet()
cp_motions = shape_dvs.add('cp_motions', build_cp_motion_dv(model.ffd.baseline_coefficients, dv_spec))
coefficients = model.ffd.apply_cp_motions(model.ffd.baseline_variable(), cp_motions,
                                          cp_dv_directions(dv_spec))

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))

params = shape_dvs.flat_variable()
mesh_motion = model.mesh_warper.evaluate(model.ffd.wall_displacement(coefficients), params)
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 M, which have derivatives and may be used as objectives and constraints, and the coefficients c_l, c_d, c_m and c_m_alpha, which are for monitoring only. The flattened shape parameter vector params is passed along for bookkeeping; it carries no derivatives.

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.

Output files

Written to the working directory, named with filename_suffix:

  • FOM_solution_*.bp, FOM_pressure_*.bp, mesh_deformation_*.bp: solution, pressure and mesh deformation per evaluation, in ADIOS2/VTX format (open in ParaView);

  • output_meshtags.xdmf: the boundary tags, to check the inflow/outflow/wall/symmetry split;

  • memory_logs/: peak memory per MPI rank.

model.data_store collects per-evaluation results; model.data_store.write_store_to_numpy_file(save_folder, save_filename) saves them to an existing folder.