dragonfly_sim.core.Euler_model

Classes

CompressibleEulerModel

The p = 0 DG (cell-centred finite-volume) compressible Euler model, and

Module Contents

class dragonfly_sim.core.Euler_model.CompressibleEulerModel(mesh, poly_order=0, gamma=1.4, flux_function=hll_flux, extra_variables=())

The p = 0 DG (cell-centred finite-volume) compressible Euler model, and the host every other flow model in the package is built on.

State: the conservative mean flow (rho, rho u, rho E), n_mean = 2 + dim components per cell, optionally followed by extra_variables transported scalars (n_state components in all). The scalars are convected with the mean flow (NS_utils.extended_flux) and are not positivity-limited. With no extra variables every form below is the plain Euler one.

Plug-in points. CompressibleRANSModel fills these; for the plain Euler model they are all empty or None, and the residual, Jacobian and derivative forms are the same expressions as without them.

face_traces () -> (U+, U-) interior-facet traces of the

inviscid flux (MUSCL reconstruction)

residual_terms callables () -> UFL form added to the residual boundary_terms callables (U_b, tag, kind) -> UFL form or None,

added to every boundary condition

boundary_scalars (kind, U, inflow) -> list of the extra variables’

boundary values (kind: slip, farfield, wall); required when n_extra > 0

freestream_scalars () -> list of the extra variables’ free-stream values wall_bc (tag) -> None, replaces the slip wall in define_wall_bc wall_traction (U, n) -> traction added to p n in wall_force_density time_scale_terms callables (U) -> added to the PTC inverse time scale derived_fields DerivedField objects (utils/derived_fields.py):

fields the residual reads that are not unknowns (a reconstructed gradient, cell centres, the wall distance). Every derivative assembler adds their chain-rule terms.

problem_factory (F, J) -> nonlinear problem, replacing

NonlinearProblem_mod in set_up_solver

extra_scales () -> list of the extra variables’ typical

magnitudes (numbers or UFL expressions of Constants); required by state_scales when n_extra > 0

Every boundary condition records (U_b, tag, kind) in boundary_states.

Time integration (core/time_integration.py): use_ptc makes solve_system a pseudo-transient continuation solve; enable_time_stepping() adds a BDF physical time term for unsteady runs. Both live in newton_forms() only, so self.F stays the steady residual that the adjoint differentiates.

classmethod _apply_ksp_options(ksp, ksp_options)

Override KSP settings through the PETSc options DB, after the model’s own: a dict of option names without the prefix, e.g. {“ksp_type”: “dgmres”, “ksp_gmres_restart”: 400}. The overrides last until the next KSP is configured here (the preconditioner reads its options at the first solve), then are deleted. An option the PETSc build does not know (a KSP type it was not built with) raises a PETSc error.

_assemble_physical_residual()

Re-assemble F_physical into self.residual_vec_physical and return its norm. Derived fields (e.g. a reconstructed gradient) are brought up to date with the current state first.

_boundary_state(U_b_mean, kind, inflow=None)

The full boundary state: the mean-flow state and, with extra variables, their values from the boundary_scalars plug-in.

static _destroy_condensed(condensed)

_ensure_residual_vector()

static _form_key(form)

Cache key of a form: its signature AND the identity of every coefficient and constant in it. The signature alone is the same for two forms that differ only in which Constant they hold (e.g. the lift and drag directions), and a compiled form is bound to its objects.

_freestream_scalars()

_inflow_direction_numpy(alpha)

_inflow_direction_ufl()

_lambda_max_expr(U)

The largest wave speed |u| + c of the state U (mean flow), with rho and p floored as in Euler_utils.signed_wave_speeds.

_make_time_integrator()

_matrix_free_adjoint_solver(prefix, rtol, atol, max_it, gmres_restart, ksp_options=None)

_release_condensed_jacobian()

Free assemble_dRdu_mat’s cached condensed dR/dU; rebuilt when next needed.

_report_max_residual_location()

Print the magnitude, solution component and physical coordinate of the largest entry in the current residual_vec_physical.

Collective: every rank participates in the allgather, and the lookup that decides which rank owns the max dof is derived from a collective Vec.max(), so all ranks agree on the answer by construction.

_report_positivity_limiting_node(x, dx, theta, initial_theta, min_theta=1e-06)

Print which node bounded the Newton step length.

_return_freed_memory()

Collect garbage, flush petsc4py’s queue of dropped parallel objects and hand freed heap back to the OS (glibc malloc_trim). Rebuilding the multi-GB 3D RANS matrices on every solve otherwise fragments the heap, and the resident size creeps up ~2 GB per optimization iteration.

_write_output_fields(write_counter)

Interpolate the output fields off the current u_vec and write one frame. u_vec must already hold the state being exported. Collective.

add_boundary_condition(meshtag, U_b, kind, inviscid_flux='llf')

Weakly impose the boundary state U_b on ds(meshtag) and register it.

inviscid_flux=”llf” is the one-sided Lax-Friedrichs boundary flux against U_b (Euler_utils.boundary_flux); “exact” is the physical flux of U_b itself, for states that already carry the upwinding (the characteristic far field) or that must have no mass flux (no-slip wall). The boundary_terms plug-ins add their terms (e.g. viscous fluxes) to the same integral.

alpha_derivative(form)

dJ/dalpha of a scalar functional J at fixed U: its explicit alpha dependence plus that of the derived fields (e.g. a reconstructed gradient reading the alpha-dependent far-field states). A number. Collective.

apply_krylov_solver_settings(ksp, max_it=100, gmres_restart=None, monitor_convergence=True, rtol=0.001, atol=None)

Configure ksp through the PETSc options DB.

rtol/atol are parameters rather than literals because this method is shared by three very different solves. The DEFAULTS are the forward Euler Newton inner KSP’s settings – a loose rtol=1e-3 carrying the Eisenstat-Walker forcing (see set_up_solver’s atol comment), and no atol at all, i.e. PETSc’s 1e-50 floor. The adjoint and tangent solves pass the model’s standardized linear-solve pair in instead.

atol=None means “leave the key unset” and NOT “set zero”: writing 0 would be a real change from PETSc’s 1e-50 default.

Both of those call sites previously called ksp.setTolerances() AFTER this method, silently overriding the rtol written here – so this method’s rtol only ever governed the forward solve, which was readable from neither end. Passing tolerances in keeps one source of truth per KSP.

assemble_dRdalpha_vec(u_vec_entries=None)

The derivative dR/d(alpha) of the converged residual, assembled as a Vec.

A VECTOR, not a matrix, because alpha is a single global scalar: the derivative of the residual with respect to it is one column, and ufl.diff against the stored alpha_var yields a rank-1 form that assembles straight into it. Nothing has to be contracted afterwards.

Note this is ufl.diff, not the ufl.derivative used by every other assembler here. ufl.derivative only accepts Coefficients and SpatialCoordinate and rejects a Constant outright; ufl.diff against a ufl.variable-wrapped Constant is the route that works. See define_farfield_conditions for why alpha_var must be the single stored instance.

Sparse in practice: alpha only enters the residual through the free-stream velocity in the far-field boundary state, so the result is supported on far-field-adjacent dofs alone. The interior terms and the slip wall (which uses the facet normal) carry no alpha dependence.

Caches its own form and Vec.

assemble_dRdu_mat(u_vec_entries=None)

The forward Newton Jacobian, assembled.

Caches its own form/matrix, mirroring compute_dRdxmesh_mat’s convention. With a reconstructed gradient this is the exact condensed dR/dU of the steady residual (an aij Mat).

compute_L2_norm(input_vec, normalization_vec=None, wall_tag=None)

L2 norm of the state-space function with the owned entries input_vec, over the domain, or over the boundary facets tagged wall_tag; every component of the state is included. With normalization_vec, also the norm relative to that function’s: returns (norm, relative norm). Collective.

compute_dRdu_mat(u_vec_entries=None)

The forward Newton Jacobian, dR/du, as a UFL expression.

compute_dRdxmesh_mat(u_vec_entries, V)

compute_initial_conditions_from_farfield()

compute_interior_integral_terms()

compute_weakform()

define_elements_functionspaces()

define_farfield_bc(meshtag, transverse=None)

Characteristic (Riemann-invariant) far field on ds(meshtag), for the whole outer boundary: no inflow/outflow division is made, each facet takes its character from the sign of the characteristic speeds, so the angle of attack can change freely. Outgoing acoustic waves leave the domain (exactly at normal incidence) instead of reflecting off a fixed state. The inviscid flux is the physical flux of the characteristic state, which already carries the upwinding; extra variables take their free-stream values where the flow enters and the interior values where it leaves.

transverse (core/farfield.py:TransverseFarfield) adds the transverse correction to the incoming invariant: the ‘riemann2’ far field of unsteady runs.

The state depends on alpha through the free-stream velocity, so dR/dalpha follows automatically.

define_farfield_conditions(farfield)

Set the free-stream state the far-field boundary condition relaxes to.

farfield is {‘rho’, ‘p’, ‘M’, ‘alpha’}: density, pressure, Mach number and angle of attack (rad). This is all the user states; the boundary treatment (inflow, outflow) follows from it per facet in define_farfield_bc, through the characteristic state.

NOTE: the far field is the subsonic characteristic one, so M < 1. A supersonic free stream needs a different boundary state in define_farfield_bc (all variables from the free stream at inflow, all from the interior at outflow); this is the one place that would change.

define_slipwall_bc(meshtag)

define_sol_export_files(func_space=None, file_name_addendum=None, mesh_deformation=False)

Open FOM_solution_<suffix>.bp, one file holding every output field (output_field_expressions) as a separate, named field; the mesh is stored once per frame for all of them. Collective.

It is the ALIGNED series: exactly one frame per design evaluation, written by write_solution_output() at eval_idx, so step k here is step k in mesh_deformation_<suffix>.bp. It carries the accepted solution – see DG_windtunnel_model.solve_residual_equations. mesh_deformation=True writes each frame on the deformed mesh it was solved on. func_space is unused (the fields define their own spaces); it is kept for callers.

define_trial_testfunctions()

define_wall_bc(meshtag)

The wall of the body: a slip wall, unless a wall_bc plug-in (e.g. the RANS model’s no-slip wall) replaces it.

static destroy_adjoint_operator(ksp)

Destroy a KSP of linearized_solver and, for a matrix-free one, what its operator’s context owns.

enable_time_stepping(dual_time=False)

Add the BDF physical time term to the Newton forms (and, with dual_time, the pseudo-time term as well). Call before set_up_solver. Returns the model’s TimeIntegrator.

incoming_state_converged()

True, with solver_converged set, if the current u_vec already satisfies residual_conv_limit. Needs only the residual form, not the solver: the caller can skip set_up_solver and solve_system. The incoming state routinely is converged, since every evaluation that re-solves at an unchanged design arrives here converged, and a full Newton step would walk it backwards. Leaves the residual in residual_vec_physical either way. Collective.

interpolate_solution_vector(interpolant=None)

linearized_solver(transpose=False, prefix='lin_', rtol=1e-10, atol=1e-11, max_it=1000, gmres_restart=200, release_operands=False, matrix_free=False, ksp_options=None)

The KSP of dR/dU x = rhs (or its transpose, for the adjoint) at the current state, for models whose dR/dU is a condensed product (has_condensed_jacobian). FGMRES on the exact operator, right preconditioned with ASM/ILU of the compact A_UU part – the forward solve’s compact_pc pairing – or MUMPS with direct_linear_solver. The transposes are formed explicitly. With release_operands and transpose, the condensed dR/dU (and its blocks) is freed once the transposes exist, before the preconditioner is set up: only the transposes are used from then on (release_solver). The KSP holds its operators and, after its first solve, the preconditioner’s factors, so it solves any number of right-hand sides at this state; the caller destroys it. Collective.

matrix_free (with transpose): the operator is the python Mat of TransposedCondensedOperator, A_UU^T x + B^T (A_UG^T x), which never forms A or the product A_UG B; its context owns the matrices and is destroyed with destroy_adjoint_operator(ksp).

mean_flow(U)

The mean-flow slice (rho, rho u, rho E) of a state expression.

mesh_derivative(form)

dJ/dx of a scalar functional J through the derived fields: a Vec on the mesh’s coordinate_space.

newton_forms()

The residual/Jacobian pair Newton drives: model.F with the time terms the model uses (core/time_integration.py). Without PTC or time stepping this is (self.F, self.J), the plain steady solve.

output_field_expressions()

The fields of the solution file, as (name, value shape, UFL expression of u_vec): the conservative state as separate fields – density rho, momentum rho*u (a vector), energy rho*E, and each transported scalar under its own name (e.g. rho_nu_tilde) – then the velocity u (a vector) and the pressure p.

positivity_step_length(x, dx, tau=1.0)

Step length for the update x - theta dx: (theta, cap), with the step-length cap = min((tau / (||dx|| / max(||x||, 1)))^0.5, 1) and theta <= cap the largest step keeping density and pressure positive in every cell (transported extra variables are not limited). Used by the Newton step limiter and by the reduced-order solve (core/reduced_order_model.py). Collective.

relative_residual_norm(full_norm)

eta = ||r|| / ||u||, the dimensionless residual measure.

Absolute residual norms are not comparable across problems: they scale with sqrt(n_dofs) and with the magnitude of the conserved variables, so a threshold calibrated on one mesh is meaningless on the next. Dividing by the current ||u|| removes both, which is what makes a single tolerance transferable between the 2D and 3D cases.

Collective (Vec.norm). Returns nan rather than raising when u_vec has blown up, so a diverged solve still prints a number.

release_derived_caches()

Free what the derived fields cache between calls (Jacobians kept while the state, alpha and mesh are unchanged). Collective.

release_solver()

Free the forward solver (SNES, KSP and its preconditioner), the Newton problem’s matrices, and the cached condensed dR/dU of assemble_dRdu_mat. set_up_solver and assemble_dRdu_mat rebuild them when next needed. For memory-bound runs, e.g. a 3D RANS shape optimization: the forward solve’s matrices and the adjoint’s then never coexist. Collective.

set_angle_of_attack(alpha)

set_mach(M)

set_u_vec(u_vec_entries)

set_up_solver()

solve_linearized(rhs, out, **kwargs)

Solve dR/dU x = rhs once (see linearized_solver for the arguments): returns the KSP (iterations, residual norm and reason are read off it).

solve_system()

Newton solve of the Euler system, starting from the current u_vec.

SNESNewtonSolver runs in passes of fom_newton_inner_max_it steps, and the physical residual is tested after each pass, until it drops below residual_conv_limit or the max_newton_iterations budget is spent. With use_ptc it is instead a pseudo-transient continuation solve (TimeIntegrator.solve_steady), max_newton_iterations pseudo steps of one Newton step each. Returns the physical residual vector (residual_vec_physical).

state_derivative(form)

dJ/dU of a scalar functional J (rank-0 UFL form), through the derived fields: a Vec in the state layout.

state_scale_expressions()

Typical magnitudes of the state components: the free-stream [rho, rho U (per direction), p/(gamma - 1) + rho U^2/2] of the mean flow, then the extra_scales plug-in’s values for the extra variables (numbers or UFL expressions of Constants).

state_scales()

state_scale_expressions evaluated to floats, one per state component. Collective (UFL expressions are evaluated by assembly).

steady_residual()

Within this block the Newton problem (self.problem) assembles the steady residual and its Jacobian: the pseudo-time and physical time terms of core/time_integration.py are switched off (their Constants set to zero) and restored on exit. A no-op without time terms.

update_derived_fields()

Bring every derived field up to date with the current state. Collective; a no-op without derived fields.

update_geometry_fields()

Refresh what depends on the mesh coordinates alone (cell centres, wall distance, the reconstruction’s cell volumes). Runs automatically after every Mesh.apply_node_motions / reset_nodes once the model has derived fields (see compute_weakform). A no-op for the plain Euler model.

wall_force_density(U, n)

Force per unit wall area on the body: pressure, plus the wall_traction plug-in’s viscous traction when there is one.

write_solution_output(write_counter, solution_array=None)

Write one aligned output frame at index write_counter.

Called once per design evaluation from DG_windtunnel_model.solve_residual_equations with write_counter = eval_idx, so FOM_solution_<suffix>.bp steps in lockstep with mesh_deformation_<suffix>.bp.

solution_array is the accepted solution of the evaluation, passed explicitly rather than read off self.u_vec. Passing None writes whatever u_vec currently holds.

F = None

J = None

_applied_ksp_option_keys

property _chain

_output_fields = []

_time_stepping = None

asm_overlap = 2

bcs = []

boundary_scalars = None

boundary_states = []

boundary_terms = []

derived_fields = []

dimensions

direct_linear_solver = False

export_solutions = False

extra_scales = None

extra_variables = ()

face_traces = None

farfield_conditions = None

flux_function

fom_newton_inner_max_it = 5

fom_solution_writer = None

forms

freestream_scalars = None

functionspaces = None

gamma = 1.4

property has_condensed_jacobian

True when dR/dU is a condensed product (a derived field depends on U), which the default ASM/ILU-on-dR/dU solves cannot precondition; the linearized solves then go through solve_linearized.

ilu_levels = 2

initial_conditions = None

krylov_diagonal_scale = True

krylov_rtol = 0.001

krylov_type = 'gmres'

last_step_theta = 1.0

mesh_obj

n_extra = 0

n_mean

n_state

poly_order = 0

problem_factory = None

ptc_settings

ptc_warm_start = True

report_positivity_limiter = True

residual_terms = []

solver_converged = False

solver_is_set_up = False

time_integrator = None

time_scale_terms = []

use_ptc = False

wall_bc = None

wall_traction = None