dragonfly_sim.core.Euler_model

Classes

CompressibleEulerModel

Module Contents

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

_assemble_physical_residual()

Re-assemble F_physical into self.residual_vec_physical and return its norm.

_inflow_direction_numpy(alpha)

_inflow_direction_ufl()

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

_write_derived_fields(solution_writer, pressure_writer, write_counter)

Interpolate the derived fields off the current u_vec and write them.

u_vec must already hold the state being exported.

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_inlet_outlet_conditions for why alpha_var must be the single stored instance.

Sparse in practice: alpha only enters the residual through u_in in the subsonic inflow BC, so the result is supported on inlet-adjacent dofs alone. The interior terms, the outflow BC (which prescribes p) 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.

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

compute_interior_integral_terms()

compute_weakform()

define_elements_functionspaces()

define_inlet_outlet_conditions(bc_dict)

define_slipwall_bc(meshtag)

define_sol_export_files(func_space, file_name_addendum=None)

define_subsonic_inflow_bc(meshtag)

define_subsonic_outflow_bc(meshtag)

define_trial_testfunctions()

interpolate_solution_vector(interpolant=None)

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.

set_angle_of_attack(alpha)

set_mach(M)

set_u_vec(u_vec_entries)

set_up_solver()

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. Returns the physical residual vector (residual_vec_physical).

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

alpha_tagging = None

alpha_tagging_tol

asm_overlap = 2

bcs = []

boundary_conditions = None

dimensions

export_solutions = False

flux_function

fom_newton_inner_max_it = 5

fom_pressure_writer = None

forms

functionspaces = None

gamma = 1.4

ilu_levels = 2

initial_conditions = None

mesh_obj

poly_order = 0

pressure_func = None

report_positivity_limiter = True

solver_converged = False

solver_is_set_up = False