Python API¶

enable_x64 explicitly enables JAX float64 arrays process-wide. Call it before constructing meshes or tracing functions; importing LMhdX does not change precision. enable_x64, a case dtype, ChannelProblem and Q2DProblem also set JAX’s jax_default_matmul_precision to 'highest' when it is unset. Float32 matrix products on Ampere GPUs then run in float32 rather than TensorFloat-32, which was 3e-4 from float64 on an RTX A4000. A value you set yourself, in jax.config or through JAX_DEFAULT_MATMUL_PRECISION, is kept. Case factories and TOML [case] accept dtype="float32" or dtype="float64" (default). During the transition, constructing a float64 case with x64 disabled enables it with a DeprecationWarning; explicit activation avoids that warning. The case dtype controls fully developed mesh, material, initial/restart and output fields, including compiled derivatives. Float32 is qualified only for the tested low-Hartmann cases, not high-Hartmann validation. Supplied meshes are cast to the case dtype; casting cannot recover precision lost during their construction. See JAX’s dtype and x64 contract.

import lmhdx

lmhdx.enable_x64()
case = lmhdx.make_hartmann_case(ha=2, ny=32, nz=32, dtype="float32")
u, phi, jy, jz, lorentz = lmhdx.solve_fully_developed_fields(case)

The package root is the small, stable convenience surface. Advanced workflows live in the module that owns their concepts.

Root API¶

Area	Names
Staggered core	`ChannelProblem`, `duct_problem`, `solve_steady_state`, `advance`
Cases and solves	`make_hartmann_case`, `make_shercliff_case`, `make_hunt_case`, `make_q2d_case`, `solve_fully_developed_fields`, `evolve_q2d`, `Q2DProblem`, `solve`
Meshes	`generate_rect_duct_mesh_from_faces` (the mesh a case’s solution reports)
Wall models	`WallLayer`, `wall_conductance_ratio`, `effective_pinhole_conductance_ratio`, `tangential_stack_conductance_ratio`, `normal_stack_leakage_ratio`, `equivalent_single_layer`, `nested_wall_layer_resolution_summary`
Units	`dynamic_to_kinematic_viscosity`, `kinematic_to_dynamic_viscosity`, `hartmann_number`, `reynolds_number`, `interaction_parameter`, `magnetic_reynolds_number`, `magnetic_field_from_hartmann`
Evidence	The energy budget of `lmhdx.core3d` and the analytical, conservation, and packaged benchmark tools in `lmhdx.io`
Runtime	`enable_compilation_cache`

solve(model) accepts a ChannelProblem, CaseSpec or Q2DProblem. A steady fully developed CaseSpec runs on the staggered core through lmhdx.fully_developed, as does solve_fully_developed_fields, and a transient CaseSpec runs implicit Euler steps there too. A duct with an inlet and an outlet is solved by lmhdx.axial.solve_open_duct. A pipe with an inlet and an outlet on the polar grid (lmhdx.axial.fringe_pipe) is solved in the Stokes limit by lmhdx.axial.solve_open_pipe.

duct_problem(hartmann=..., cells=..., wall_conductance=...) builds a square insulating or Hunt duct with meshes that resolve the layers that exist – a/Ha against the walls normal to the field, a/sqrt(Ha) against the others – and solve_steady_state finds its steady state as a differentiable root. advance runs the same physics as one compiled trajectory when the transient is what is wanted. The remaining result types expose converged, status, steps, residual, fields, and diagnostics; specialized solve functions provide restart, progress, logging, and timing hooks in their owning modules.

Cases: schema, builders, meshes, units and walls¶

Fully developed duct cases: their schema, materials, meshes and the solve entry point.

A CaseSpec (geometry, regions, field, conditions, solver and output choices) is solved on the staggered core by lmhdx.fully_developed: steady cases by solve_fully_developed(), transient ones by solve_fully_developed_transient(). This module also holds the material and nondimensional relations (Hartmann, Reynolds and interaction numbers, wall conductance ratios and wall stacks), the structured mesh a solution reports, tabulated imposed fields, TOML run configurations and the streaming solver log.

class lmhdx.cases.BoundaryCondition(name: 'str', kind: 'BoundaryKind', value: 'float | tuple[float, float, float] | None' = None, region: 'str | None' = None, axis: 'str | None' = None, side: 'str | None' = None)¶: Bases: object

class lmhdx.cases.CaseSpec(name: 'str', geometry: 'GeometrySpec', regions: 'tuple[RegionSpec, ...]', magnetic_field: 'MagneticFieldSpec', boundary_conditions: 'tuple[BoundaryCondition, ...]', time_stepper: 'TimeStepperConfig', solver: 'SolverConfig' = <factory>, output: 'OutputSpec' = <factory>, forcing: 'float' = 1.0, initial_velocity: 'float' = 0.0, notes: 'str' = '', dtype: 'str' = 'float64')¶: Bases: object

class lmhdx.cases.Diagnostics(residual_history: 'jnp.ndarray', courant_like: 'jnp.ndarray', ohmic_power: 'jnp.ndarray', time_history: 'jnp.ndarray' = <factory>, u_max_history: 'jnp.ndarray' = <factory>, mean_velocity_history: 'jnp.ndarray' = <factory>, applied_forcing_history: 'jnp.ndarray' = <factory>, current_max_history: 'jnp.ndarray' = <factory>, face_current_max_history: 'jnp.ndarray' = <factory>, emf_max_history: 'jnp.ndarray' = <factory>, lorentz_max_history: 'jnp.ndarray' = <factory>, face_lorentz_max_history: 'jnp.ndarray' = <factory>, potential_residual_history: 'jnp.ndarray' = <factory>, potential_iterations_history: 'jnp.ndarray' = <factory>, linear_residual_history: 'jnp.ndarray' = <factory>, linear_iterations_history: 'jnp.ndarray' = <factory>, volumetric_flow_rate_history: 'jnp.ndarray' = <factory>, mean_current_magnitude_history: 'jnp.ndarray' = <factory>, lorentz_power_history: 'jnp.ndarray' = <factory>, div_current_max_history: 'jnp.ndarray' = <factory>, charge_balance_residual_history: 'jnp.ndarray' = <factory>, gauge_residual_history: 'jnp.ndarray' = <factory>, interface_current_residual_history: 'jnp.ndarray' = <factory>)¶: Bases: object

class lmhdx.cases.GeometrySpec(kind: 'GeometryKind', width: 'float', height: 'float', length: 'float' = 1.0, nx: 'int' = 1, ny: 'int' = 64, nz: 'int' = 64, wall_thickness: 'tuple[float, float, float, float]' = (0.0, 0.0, 0.0, 0.0), wall_cells: 'tuple[int, int, int, int]' = (0, 0, 0, 0), hartmann_layer_cells: 'int | None' = None, wall_model: 'WallModel' = 'auto')¶: Bases: object

class lmhdx.cases.LoggingSpec(enabled: 'bool' = True, banner: 'bool' = True, print_footer: 'bool' = True, flush: 'bool' = True, step_stride: 'int' = 1)¶: Bases: object

class lmhdx.cases.MHDState(u: 'jnp.ndarray', phi: 'jnp.ndarray', jy: 'jnp.ndarray', jz: 'jnp.ndarray', lorentz_x: 'jnp.ndarray', time: 'float', residual: 'float')¶: Bases: object

class lmhdx.cases.MagneticFieldSpec(kind: 'MagneticFieldKind', value: 'tuple[float, float, float] | None' = None, fn: 'Callable[[jnp.ndarray, jnp.ndarray], jnp.ndarray] | None' = None, table_path: 'str | None' = None, ramp_start: 'float' = 0.0, ramp_duration: 'float' = 0.0)¶: Bases: object

exception lmhdx.cases.NumericalFailure¶

Bases: RuntimeError

Raised when a solver produces nonfinite numerical state.

class lmhdx.cases.OutputSpec(directory: 'str | None' = None, write_paraview: 'bool' = True, write_csv_profiles: 'bool' = True, write_npz: 'bool' = True, write_json_summary: 'bool' = True, write_plots: 'bool' = False, copy_input_file: 'bool' = True, history_stride: 'int' = 0)¶: Bases: object

class lmhdx.cases.RegionSpec(name: str, kind: Literal['fluid', 'solid'], conductivity: float, density: float | None = None, viscosity: float | None = None)¶

Bases: object

Material-region specification.

viscosity is kinematic viscosity nu in m^2/s. Dynamic viscosity mu should be converted before constructing a case.

class lmhdx.cases.RestartLogInfo(enabled: 'bool' = False, path: 'str | None' = None, start_time: 'float' = 0.0, reset_histories: 'bool' = True)¶: Bases: object

class lmhdx.cases.RestartSpec(enabled: 'bool' = False, path: 'Path | None' = None, reset_histories: 'bool' = True, write_restart: 'bool' = False, restart_filename: 'str | None' = None)¶: Bases: object

class lmhdx.cases.RunConfig(case: 'CaseSpec', logging: 'LoggingSpec' = <factory>, restart: 'RestartSpec' = <factory>, input_path: 'Path | None' = None)¶: Bases: object

class lmhdx.cases.Solution(mesh: 'StructuredMesh', state: 'MHDState', diagnostics: 'Diagnostics', case_name: 'str', converged: 'bool | None' = None, status: 'str' = 'not_recorded', steps: 'int' = 0)¶

Bases: object

property fields: MHDState¶: Final fully developed MHD state.

property residual: float¶: Terminal normalized solver residual.

class lmhdx.cases.SolverConfig(kind: 'SolverKind' = 'fully_developed_inductionless', mode: 'SolveMode' = 'steady', time_scheme: 'TimeSchemeKind' = 'implicit_euler')¶: Bases: object

class lmhdx.cases.SolverStepRecord(step_index: 'int', time: 'float', u_max: 'float', mean_velocity: 'float', current_max: 'float', lorentz_max: 'float', residual: 'float', potential_residual: 'float', potential_iterations: 'float', linear_residual: 'float', linear_iterations: 'float', applied_forcing: 'float', courant_like: 'float', ohmic_power: 'float', volumetric_flow_rate: 'float', div_current_max: 'float', gauge_residual: 'float', interface_current_residual: 'float', charge_balance_residual: 'float' = 0.0, potential_initial_residual: 'float' = 0.0, linear_initial_residual: 'float' = 0.0)¶: Bases: object

class lmhdx.cases.StructuredMesh(x_faces: 'jnp.ndarray', y_faces: 'jnp.ndarray', z_faces: 'jnp.ndarray', geometry: 'str' = 'rect_duct', point_coordinates: 'jnp.ndarray | None' = None, fluid_mask: 'jnp.ndarray | None' = None, sigma: 'jnp.ndarray | None' = None, region_ids: 'jnp.ndarray | None' = None, region_names: 'tuple[str, ...]' = ())¶: Bases: object

class lmhdx.cases.TimeStepperConfig(dt: 'float', t_final: 'float', max_steps: 'int')¶: Bases: object

class lmhdx.cases.WallLayer(name: str, conductivity: float, thickness: float, cells: int = 1)¶

Bases: object

One solid layer in a fluid-facing electrical wall stack.

lmhdx.cases.dynamic_to_kinematic_viscosity(dynamic_viscosity: float, density: float) → float¶: Return kinematic viscosity nu = mu / rho in m^2/s.

lmhdx.cases.effective_pinhole_conductance_ratio(*, intact_conductance_ratio: float, metal_conductance_ratio: float, pinhole_fraction: float) → float¶: Return the area-weighted conductance ratio for a pinholed coating.

lmhdx.cases.equivalent_single_layer(layers: Sequence[WallLayer], *, name: str = 'equivalent_wall') → WallLayer¶: Return one layer with the same tangential surface conductance.

lmhdx.cases.generate_rect_duct_mesh_from_faces(*, y_faces: Array, z_faces: Array, length: float = 1.0, nx: int = 1) → StructuredMesh¶: Build a rectangular duct mesh from explicit cross-section faces.

lmhdx.cases.hartmann_number(*, magnetic_field: float, length_scale: float, conductivity: float, density: float, kinematic_viscosity: float) → float¶: Return Ha = B a sqrt(sigma / (rho nu)).

lmhdx.cases.interaction_parameter(*, magnetic_field: float, length_scale: float, conductivity: float, density: float, velocity: float) → float¶: Return N = sigma B^2 a / (rho U).

lmhdx.cases.kinematic_to_dynamic_viscosity(kinematic_viscosity: float, density: float) → float¶: Return dynamic viscosity mu = rho nu in Pa s.

lmhdx.cases.magnetic_field_from_hartmann(*, hartmann: float, length_scale: float, conductivity: float, density: float, kinematic_viscosity: float) → float¶: Return B from a target Hartmann number using kinematic viscosity.

lmhdx.cases.magnetic_reynolds_number(*, velocity: float, length_scale: float, conductivity: float, magnetic_permeability: float = 1.2566370614359173e-06) → float¶: Return Rm = mu0 sigma U a.

lmhdx.cases.make_divergence_free_cross_section_field(*, width: float, height: float, base_bz: float, perturbation: float = 0.15)¶: Return an analytic cross-sectional field with dBy/dy + dBz/dz = 0.

lmhdx.cases.make_hartmann_case(ha: float = 20.0, width: float = 2.0, height: float = 2.0, ny: int = 96, nz: int = 96, conductivity: float = 1.0, density: float = 1.0, viscosity: float = 1.0, output_dir: str | None = None, dtype: str = 'float64') → CaseSpec¶: Build an insulating rectangular Hartmann-duct reference case.

lmhdx.cases.make_hunt_case(ha: float = 20.0, width: float = 2.0, height: float = 2.0, ny: int = 72, nz: int = 72, wall_cells: int = 8, wall_thickness: float = 0.1, insulator_cells: int | None = None, insulator_thickness: float | None = None, fluid_conductivity: float = 1.0, wall_conductance_ratio: float = 0.05, wall_conductivity: float | None = None, insulator_conductivity: float | None = None, insulator_conductivity_ratio: float = 1e-12, density: float = 1.0, viscosity: float = 1.0, output_dir: str | None = None, dtype: str = 'float64') → CaseSpec¶: Build a Hunt duct with conducting Hartmann and insulating side walls.

lmhdx.cases.make_shercliff_case(ha: float = 20.0, width: float = 2.0, height: float = 2.0, ny: int = 96, nz: int = 96, conductivity: float = 1.0, density: float = 1.0, viscosity: float = 1.0, output_dir: str | None = None, dtype: str = 'float64') → CaseSpec¶: Build an all-insulating rectangular Shercliff-duct case.

lmhdx.cases.nested_wall_layer_resolution_summary(layers: Sequence[WallLayer], *, minimum_cells_per_layer: int = 3) → dict[str, object]¶: Return mesh-resolution metrics for a wall-layer stack.

lmhdx.cases.normal_leakage_ratio(*, coating_conductivity: float, coating_thickness: float, fluid_conductivity: float, length_scale: float) → float¶: Return normal shunt ratio g_perp.

lmhdx.cases.normal_stack_leakage_ratio(layers: Sequence[WallLayer], *, fluid_conductivity: float, length_scale: float) → float¶: Return the normal leakage ratio for layers in series.

lmhdx.cases.require_finite(stage: str, **values) → None¶: Raise with field names when numerical output is nonfinite.

lmhdx.cases.reynolds_number(*, velocity: float, length_scale: float, kinematic_viscosity: float) → float¶: Return Re = U a / nu.

lmhdx.cases.solve(model: ChannelProblem | CaseSpec | Q2DProblem) → SteadySolution | Solution | Q2DResult¶

Solve a duct, a fully developed case, or a Q2D problem.

A lmhdx.core3d.ChannelProblem goes to the staggered core’s steady solve, lmhdx.steady.solve_steady_state(), and so does a steady CaseSpec, through lmhdx.fully_developed.solve_fully_developed(), which reports it on the case’s cross-section; each is compiled once per problem. A transient CaseSpec runs implicit Euler steps on the core, lmhdx.fully_developed.solve_fully_developed_transient(). A duct with an inlet and an outlet is solved by lmhdx.axial.solve_open_duct().

lmhdx.cases.tangential_stack_conductance_ratio(layers: Sequence[WallLayer], *, fluid_conductivity: float, length_scale: float) → float¶: Return the thin-wall tangential ratio for layers in parallel.

lmhdx.cases.wall_conductance_ratio(*, wall_conductivity: float, wall_thickness: float, fluid_conductivity: float, length_scale: float) → float¶: Return thin-wall tangential conductance ratio c.

Inertialess core flow¶

Inertialess core-flow model of a thin-walled rectangular duct in a varying field.

At large Hartmann number and interaction parameter, inertia and viscosity are confined to layers and the core obeys grad p = j x B and Ohm’s law. For a field B_y(x) across a duct of half-width 1 (side walls at z = -1, 1) and half-height a (Hartmann walls at y = -a, a), the pressure is constant along field lines and the three-dimensional core reduces to three functions of two variables (Hua, Walker, Picologlou & Reed, ANL/FPP/TM-228, 1988, eqs. 4a-4c): the core pressure p(x, z), the Hartmann-wall potential phi_t(x, z) and the side-wall potential phi_s(x, y). With beta = 1/B and K = beta^2 + a^2 beta'^2 / 3, on the quadrant 0 <= y <= a, -1 <= z <= 0:

d/dx(beta^2 dp/dx) + K d2p/dz2 = beta' dphi_t/dz          (4a, no normal flow at y = a)
c_t (d2/dx2 + d2/dz2) phi_t    = a beta' dp/dz             (4b, charge in the Hartmann wall)
c_s (d2/dx2 + d2/dy2) phi_s    = -beta dp/dx (x, -1)       (4c, charge in the side wall)

with phi_t = 0 and dp/dz = 0 at z = 0, dphi_s/dy = 0 at y = 0, the corner conditions phi_t(x, -1) = phi_s(x, a) and c_t dphi_t/dz = c_s dphi_s/dy (7d, e), and the side-layer flux closure (14) K dp/dz(x, -1) = beta' phi_t(x, -1) - (beta int_0^a phi_s dy)' / a, which says that the core and side-layer flux Q = -a beta^2 int dp/dx dz - beta int phi_s dy does not change along the duct. Walker’s thin-wall limits follow in a uniform field: the fully developed gradient is -dp/dx = B^2 / (1 + a/c_t + a^2/(3 c_s)) at unit mean velocity, and c_t/(a + c_t) with perfectly conducting side walls.

The equations are the stationarity conditions of one quadratic functional, maximal in p and minimal in the potentials:

L = int int [-a/2 (beta^2 p_x^2 + K p_z^2) + c_t/2 |grad phi_t|^2 - a beta' p dphi_t/dz] dx dz
    + int int c_s/2 |grad phi_s|^2 dx dy - int p(x, -1) q dx,   q = a K dp/dz(x, -1)

which is why the coupled operator is symmetric and indefinite. It is discretized directly, so the matrix is symmetric by construction. As in TM-228 section 3 the grid is staggered in z: phi_t sits on nodes from the corner to z = 0, p at the cell centres between them; phi_s sits on nodes in y and the corner node is shared, so its molecule is split between the Hartmann wall (half a cell in z) and the side wall (half a cell in y) and carries (7d, e) without further equations. Every term is a finite-volume molecule. The pressure at the side wall is extrapolated from the first centre with the closure, p(-1) = p_1 - dz/2 dp/dz(-1), the higher-order expansion TM-228 section 3.2 asks for; in the functional this is the term dz/(4a) int q^2/K dx. The ends are fully developed (5a-f): p is given, dphi/dx = 0 is natural. The solution is then rescaled once to the imposed flow rate (6).

The pressure and potentials are solved together, never segregated: TM-228 found a segregated scheme divergent for small wall conductance. The system is small and two-dimensional, so it is solved directly (SuperLU on the host) inside jax.lax.custom_linear_solve(); derivatives with respect to the wall conductances, the field scale and the drive are exact, and the adjoint is the same factorization. beta = 1/B is floored at beta_max (TM-228 caps it at 1000), and the number of floored nodes is returned.

The model neglects inertia: its error scales as N^(-1/3) (Mistrangelo et al., Fusion Eng. Des. 173, 2021), which at ALEX B2 (N = 540) is as large as the three-dimensional excess it computes.

class lmhdx.coreflow.CoreFlow(x, *, aspect: float = 1.0, nz: int = 16, ny: int = 16)¶

Mesh and assembled structure of the TM-228 core-flow problem on a uniform x grid.

x holds the nx + 1 uniformly spaced stations (the ends are fully developed), aspect is a (the Hartmann-wall half-height over the side-wall half-width), and nz and ny are the cells across the Hartmann and side walls. Everything here is host data fixed by the mesh; the field, conductances and drive enter solve() as traced values.

operator(field, *, c_t, c_s, field_scale=1.0, beta_max=1000.0) → csr_matrix¶: Return the assembled coupled operator on the free unknowns as a host sparse matrix.

solve(field, *, c_t, c_s, field_scale=1.0, drive=1.0, mean_velocity=1.0, beta_max=1000.0) → CoreFlowResult¶

Solve for p, phi_t and phi_s together; differentiable in every traced argument.

c_t and c_s are the Hartmann- and side-wall conductance ratios: a number, an array with one value per station of x, or a callable of x returning one (see layer_conductances()). A face of the x grid conducts at the mean of its two stations, so the functional, and the operator, stay symmetric, and the derivatives in every conductance value are exact. field is B_y at the stations of x (for example from lmhdx.core3d.fringe_field() through midplane_field(), or any 1-D array), multiplied by field_scale. The inlet pressure is drive and the outlet zero; unless mean_velocity is None the solution is then rescaled once so that the mean axial velocity is mean_velocity.

class lmhdx.coreflow.CoreFlowResult(pressure_drop: jax.Array, axial_flux: jax.Array, pressure: jax.Array, phi_top: jax.Array, phi_side: jax.Array, floored_nodes: jax.Array)¶

Core-flow solution on one quadrant, scaled to the requested mean velocity.

pressure is (nx + 1, nz) at the z cell centres, phi_top is (nx + 1, nz + 1) on the z nodes (the last column is z = 0) and phi_side is (nx + 1, ny + 1) on the y nodes (the last column is the corner). axial_flux is the quadrant flux on each x face and pressure_drop is p(x_first) - p(x_last).

axial_flux: Array¶: Alias for field number 1

floored_nodes: Array¶: Alias for field number 5

phi_side: Array¶: Alias for field number 4

phi_top: Array¶: Alias for field number 3

pressure: Array¶: Alias for field number 2

pressure_drop: Array¶: Alias for field number 0

lmhdx.coreflow.fully_developed_gradient(c_t, c_s, aspect: float = 1.0)¶: Return Walker’s fully developed -dp/dx at unit field and unit mean velocity.

lmhdx.coreflow.layer_conductances(c_t, c_s, hartmann: float, field, k, *, floor: float = 1.0)¶

Return the wall conductances with the layers’ finite-Ha conductance at the local field.

At a finite Hartmann number the Hartmann layers conduct like 1/Ha of extra Hartmann-wall conductance and the side layers like k/sqrt(Ha) of extra side-wall conductance; at a station of field B the local Hartmann number is Ha |B|, so c_t + 1/(Ha |B|) and c_s + k/sqrt(Ha |B|) (#204: with these the 3-D core and the core-flow model agree within 1 % at c 0.1, Ha 2e4). k is a number, a callable of the local Hartmann number, or a table (hartmann_numbers, k_values) interpolated in log Ha, for example from side_layer_coefficient(). Where Ha |B| is below floor the correction is frozen at floor. c_t, c_s and field are numbers or one value per station; the result is two arrays for CoreFlow.solve().

lmhdx.coreflow.midplane_field(field) → tuple[ndarray, ndarray]¶: Return the cell-centre x and B_y nearest y = 0 of an lmhdx.core3d.ImposedField.

lmhdx.coreflow.side_layer_coefficient(c: float, hartmann: float, *, cells: int = 48) → float¶

Return k of the side layers in a square duct with both walls of conductance ratio c.

The fully developed gradient of LMhdX’s own 2-D solve on the staggered core (lmhdx.core3d.duct_problem(), thin walls on both axes) is set equal to Walker’s 1/(1 + 1/c_t + 1/(3 c_s)) with c_t = c + 1/Ha, which is inverted for c_s = c + k/sqrt(Ha). One steady solve; cells must resolve the layers (48 to 64 cells from Ha 400 to 2e4, #204).

Ducts with an inlet and an outlet¶

Ducts with an inlet and an outlet: the non-periodic axial direction (plan 1.9b, D26).

The axial axis is the first. Its conditions follow the published practice of HIMAG, FreeMHD, GridapMHD and the 2025 six-code benchmark, chosen so that every solve stays direct and every derivative implicit:

Inlet. The velocity is LMhdX’s own fully developed profile at the inlet field, solved on the same cross-section and scaled to the imposed flow rate. It is array-valued Dirichlet data on the inlet face (lmhdx.grid.BoundaryCondition). The flow rate is exact and the pressure drop is an output; there is no extra unknown.
Outlet. Zero axial gradient of every velocity component and p = 0. The pressure operator is then non-singular, and the axial axis stays diagonalizable (Neumann at the inlet, Dirichlet at the outlet).
No normal current at either end, dphi/dn = (u x B).n, with the potential’s gauge fixed by removing its mean.

The inlet enters as a lift. The fully developed profile carried unchanged along the whole duct is discretely divergence free, so the solution is that lift plus a correction with a zero inlet face. The correction lives in a linear space on which the Stokes-limit operator is symmetric in the face-volume inner product, the two end faces owning half a cell each, so the solve is the same preconditioned conjugate-gradient solve as a periodic duct’s (lmhdx.steady.solve_steady_state()) and differentiates the same way. Far upstream of a field change the lift is the discrete solution, which is why the uniform region carries its fully developed gradient.

With advection on (plan 1.9d) the steady problem is nonlinear. It is solved by Newton’s method from the Stokes-limit solution, each Newton update by flexible GMRES with recycling (solvax.gcrot()) right-preconditioned by the Stokes-limit solve itself, a conjugate-gradient solve run to a loose tolerance: the Oseen operator differs from the Stokes one by the transport, which is small against the Lorentz force when N = Ha^2 / Re is large, so a few outer iterations suffice. Continuation in the flow rate steps toward the target Reynolds number. The root is differentiated by the implicit function theorem (solvax.root_solve()), its tangent and transposed solves being the same preconditioned Krylov solves.

class lmhdx.axial.OpenDuctSolution(velocity: tuple[Field, Field, Field], pressure: Field, potential: Field, currents: tuple[Field, Field, Field], magnetic_field: tuple[Field, Field, Field], residual_norm: jnp.ndarray, initial_residual_norm: jnp.ndarray, iterations: jnp.ndarray)¶

Fields of an inflow-outflow duct, a pytree; pressure is the physical one, zero at the outlet.

currents: tuple[Field, Field, Field]¶: Alias for field number 3

initial_residual_norm: Array¶: Alias for field number 6

iterations: Array¶: Alias for field number 7

magnetic_field: tuple[Field, Field, Field]¶: Alias for field number 4

potential: Field¶: Alias for field number 2

pressure: Field¶: Alias for field number 1

residual_norm: Array¶: Alias for field number 5

velocity: tuple[Field, Field, Field]¶: Alias for field number 0

lmhdx.axial.axial_faces(lower: float, upper: float, core: tuple[float, float], spacing: float, growth: float = 1.1) → ndarray¶

Faces of spacing spacing over core, growing by growth per cell into the buffers.

Buffer cells stop growing at 8 * spacing; the last cell on each side absorbs the remainder, so the ends land exactly on lower and upper.

lmhdx.axial.charge_balance(solution: OpenDuctSolution, problem: ChannelProblem) → Array¶

Largest net current out of a cell, relative to the largest motional current through a cell.

The motional current sigma (u x B).n is what the potential balances; the net current is a small difference of it and the potential gradient, so it is the scale the charge equation is solved to.

lmhdx.axial.fringe_duct(*, hartmann: float, wall_conductance: float = 0.0, half_length: float = 3.0, upstream: float = 15.0, downstream: float = 10.0, spacing: float = 0.25, cells: int = 24, cells_in_layer: int = 6, flow_rate: float = 4.0, solenoidal: bool = False, advection: str = 'off', **controls) → ChannelProblem¶

The ANL fringe (TM-228) in a square duct with an inlet and an outlet.

The field falls as B_y = Ha (1 - sin(pi x / 2 x0)) / 2 over |x| <= half_length and is uniform outside; upstream and downstream half-widths of buffer (D26: 15 and 10) separate the ramp from the ends. solenoidal=False is TM-228’s field, B_y alone, which is divergence free. The conductance applies to all four walls, and the default flow rate is a unit mean velocity; half-width, density, viscosity and conductivity being one, the Reynolds number is the mean velocity, flow_rate / 4, and N = Ha^2 / Re.

lmhdx.axial.fully_developed_inlet(problem: ChannelProblem, flow_rate: float, *, tolerance: float = 1e-09, **controls) → tuple[ndarray, float]¶

Solve the fully developed flow of problem’s first cross-section at its first cell’s field.

The inlet field must be uniform over the cross-section, as it is upstream of a magnet. Returns the axial velocity on that cross-section scaled to flow_rate, and the axial pressure gradient that drives it (negative for a positive flow). The cross-section, conductances and field are the duct’s own, so upstream of any field change the three-dimensional solution reproduces it to round-off.

lmhdx.axial.mass_balance(velocity: tuple[Field, Field, Field]) → Array¶: Largest net volume flux out of a cell, relative to the largest flux through a cell.

lmhdx.axial.open_duct(problem: ChannelProblem, flow_rate: float, **controls) → ChannelProblem¶

Turn the first axis of problem into an inlet and an outlet at the imposed flow_rate.

problem supplies the mesh, walls, field and properties; its first condition and forcing are replaced. controls go to fully_developed_inlet().

lmhdx.axial.pressure_drop(pressure: Field, start: float, end: float) → Array¶: Mean pressure at start minus that at end, interpolated linearly between centres.

lmhdx.axial.solve_open_duct(problem: ChannelProblem, *, field_scale: float | Array = 1.0, tolerance: float = 1e-09, max_iterations: int = 36000, continuation: tuple[float, ...] = (1.0,), max_newton_steps: int = 12, inner_tolerance: float = 0.01, inner_iterations: int = 4000) → OpenDuctSolution¶

Solve an inflow-outflow duct; differentiable in field_scale.

In the Stokes limit, one preconditioned conjugate-gradient solve for the correction to the lift, certified on its residual (tolerance relative to the lift’s residual); a rejected solve raises eagerly and gives nonfinite fields under tracing. max_iterations follows the tolerance rule of #150 (600 restarts of 60).

With advection, Newton’s method from that solution (module docstring): continuation lists the fractions of the flow rate solved in turn, ending at 1; each Newton update is a solvax.gcrot() solve preconditioned by the Stokes-limit CG run to inner_tolerance (at most inner_iterations), and iterations then counts the outer Krylov iterations of all Newton steps. The root is certified at tolerance.

lmhdx.axial.station_flow_rates(velocity: tuple[Field, Field, Field]) → Array¶: The flow rate through every axial face.

lmhdx.axial.station_pressure(pressure: Field) → tuple[ndarray, Array]¶: The axial cell centres and the area-mean pressure over each cross-section.

Fast-diagonal Poisson and Helmholtz solvers¶

Direct Poisson solves by fast diagonalization on a tensor-product grid.

The finite-volume Laplacian of lmhdx.ops separates: on a tensor-product grid it is the Kronecker sum of three one-dimensional operators. Each of those is symmetric once the cell widths are folded in, so diagonalizing them on the host turns a Poisson solve into three tensor contractions and one elementwise divide. The cost is a handful of matrix multiplies rather than an iteration whose count depends on the Hartmann number, the result is exact to round-off rather than to a tolerance, and the solve is a linear map, so it differentiates for free.

The one-dimensional operators are read out of lmhdx.ops itself, by applying the assembled Laplacian to unit vectors on a grid that is one cell wide across the other two axes. The factorization therefore cannot drift away from the stencil the rest of the code uses.

A pure Neumann or fully periodic problem determines its solution only up to a constant. The compatible component is removed from the right-hand side and the returned field has zero volume-weighted mean.

The factorization represents the homogeneous operator. Inhomogeneous boundary data is affine, not linear, so it belongs in the right-hand side; passing a condition that carries a value is refused rather than silently linearized.

precision="mixed" runs the contractions of a float64 solve in float32 and recovers float64 accuracy by defect correction: the residual against the assembled operator is formed in float64 and solved again in float32, refinements times (solvax.iterative_refinement()). The orthogonal transforms keep their relative accuracy mode by mode, so the shifted viscous operator reaches the float64 floor after one correction; the singular Laplacian on a layer-resolving mesh contracts more slowly and takes two, which are the factory defaults. The contraction assumes true float32 matmuls; on Ampere GPUs pin jax_default_matmul_precision to "float32", because TensorFloat-32 stalls the correction near 1e-4. Float32 states are solved in float32 as before.

The fully developed pipe built on the polar solver lives in lmhdx.pipe; its public names stay importable from here.

class lmhdx.poisson.FastDiagonalHelmholtz(grid: Grid, offset: tuple[float, float, float], conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], shift: float, coefficient: float, vectors: tuple[ndarray, ndarray, ndarray], values: tuple[ndarray, ndarray, ndarray], scales: tuple[ndarray, ndarray, ndarray], slices: tuple[slice, slice, slice], precision: str = 'state', refinements: int = 1)¶

A factorized shift * I - coefficient * laplacian for one staggered position.

This is what an implicit viscous solve needs. The operator separates exactly as the pressure Laplacian does, so the same host-side eigendecomposition turns the solve into three contractions and a divide, and the step is no longer bounded by the mesh.

solve(rhs: Field) → Field¶

Return the field this operator maps to rhs.

Prescribed entries are returned as zero: they are boundary data the caller owns, not unknowns this solve may set.

class lmhdx.poisson.FastDiagonalPoisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], vectors: tuple[ndarray, ndarray, ndarray], values: tuple[ndarray, ndarray, ndarray], scales: tuple[ndarray, ndarray, ndarray], singular: bool, precision: str = 'state', refinements: int = 2, corrections: int = 0)¶

A factorized Laplacian that solves laplacian(u) = rhs in three contractions.

residual_norm(solution: Field, rhs: Field) → Array¶: Return the maximum absolute residual of a candidate solution.

solve(rhs: Field) → Field¶: Return the field whose Laplacian is rhs.

class lmhdx.poisson.FastDiagonalPolarPoisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], radial_vectors: ndarray, radial_values: ndarray, radial_scale: ndarray, axial_vectors: ndarray, axial_values: ndarray, axial_scale: ndarray, singular: bool, shift: float = 0.0, coefficient: float = -1.0, precision: str = 'state', refinements: int = 2, wall_conductance: float = 0.0)¶

A factorized polar Laplacian: one radial eigendecomposition per azimuthal mode.

solve(rhs: Field) → Field¶

Return the field this operator maps to rhs.

The operator is shift*I - coefficient*laplacian; the Poisson factory passes (0, -1), which leaves the Laplacian itself.

solve_with_wall(rhs: Field) → tuple[Field, Field]¶

Return the cell potential and, on the radial faces, the outer thin wall’s potential.

As lmhdx.ops.thin_wall_current() reads it: the outer entry is the sheet, the inner entry repeats the adjacent cell so it carries no current.

class lmhdx.poisson.FastDiagonalThinWallPoisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], vectors: tuple[ndarray, ndarray, ndarray], values: tuple[ndarray, ndarray, ndarray], scales: tuple[ndarray, ndarray, ndarray], singular: bool, precision: str = 'state', refinements: int = 2, corrections: int = 0, conductance: tuple[float, float, float] = (0.0, 0.0, 0.0), operators: tuple[ndarray, ...] = (), weights: tuple[ndarray, ...] = (), corner_gain: ndarray | None = None, pads: tuple[tuple[int, int], ...] | None = None, fractions: tuple[tuple[float, float], ...] | None = None, corner_fix: tuple | None = None)¶

The potential Laplacian closed by thin conducting walls, factorized exactly.

Each conducting axis carries a sheet node at both walls (assemble_thin_wall_operator()), so the solve is still three contractions and a divide, symmetric in the cell volumes extended by c times the wall area. Where two conducting walls meet, the corner node joins the two sheets in series, so the charge one delivers is what the other receives (Hua et al. 1988 split the corner the same way). The Kronecker sum would also let that node conduct along the edge, which has no sheet area; a rank-four Woodbury correction per mode of the third axis removes it, and vanishes when that axis does not vary.

solve(rhs: Field) → Field¶: Return the cell potential whose charge balance is rhs.

solve_with_walls(rhs: Field) → tuple[Field, tuple[Field | None, Field | None, Field | None]]¶

Return the cell potential, with zero volume mean, and each conducting axis’s sheet potentials.

Sheet potentials come back on the faces normal to their axis, as lmhdx.ops.thin_wall_current() reads them: wall entries set, interior zero.

lmhdx.poisson.assemble_axis_laplacian(grid: Grid, axis: int, condition: BoundaryCondition) → ndarray¶

Return the dense one-dimensional Laplacian lmhdx.ops applies along axis.

The other two axes are collapsed to a single cell with a homogeneous Neumann condition, which contributes nothing, so the result is exactly the stencil the three-dimensional operator uses along axis.

lmhdx.poisson.assemble_radial_laplacian(grid: Grid, condition: BoundaryCondition) → ndarray¶

Return the dense radial Laplacian the polar flux form applies.

The azimuth is collapsed to a single periodic cell, which contributes nothing because its two faces carry the same value, and the axial direction to a single cell with a homogeneous Neumann condition. The metric factors of the azimuth and the axis cancel between the face areas and the cell volume, so what is left is exactly (1/r) d/dr (r d/dr) as the production stencil discretizes it, including the zero-area face on the axis.

lmhdx.poisson.assemble_staggered_axis_operator(grid: Grid, axis: int, offset: tuple[float, float, float], condition: BoundaryCondition) → ndarray¶

Return the dense one-dimensional operator lmhdx.ops applies along axis.

As with assemble_axis_laplacian(), the other axes are collapsed to a single cell under a homogeneous Neumann condition so they contribute nothing, and the operator is read out of the production stencil rather than written a second time.

lmhdx.poisson.assemble_thick_wall_operator(grid: Grid, axis: int, layers: tuple) → tuple[ndarray, ndarray, tuple[float, float]]¶

Return the one-dimensional potential operator with a wall resolved in cells at each end, and its weights.

layers is (lower, upper): None for an insulating wall, else (ratios, widths), each cell’s conductivity over the fluid’s (one number for all) and its width, from the fluid outwards, insulated outside; adjacent cells join through their half cells in series. A wall cell is weighted by the ratio times its width, so the tangential axes conduct through it at the wall’s conductivity and the operator stays a Kronecker sum, as for assemble_thin_wall_operator(); the wall must span the fluid’s tangential extent. The fluid reaches the first wall cell through the two half cells in series. Also returned, per end, is the fraction of the potential difference across the fluid’s half cell, which FastDiagonalThinWallPoisson.solve_with_walls() uses to report the interface potential lmhdx.ops.thin_wall_current() reads.

lmhdx.poisson.assemble_thin_wall_operator(grid: Grid, axis: int, conductance: tuple[float, float]) → tuple[ndarray, ndarray]¶

Return the one-dimensional potential operator with a sheet node on each conducting wall, and its weights.

A thin wall of conductance ratio c has a potential of its own: the fluid reaches it across the half cell, j_n = (phi_P - phi_w)/(h_P/2), and the sheet carries that current along itself, j_n = -c lap_t phi_w (Walker’s condition). The sheet is one more node on the wall-normal axis, weighted by c times the wall’s face measure where a cell is weighted by its own measure, so the tangential axes act on it as on a cell and the operator stays a Kronecker sum. The rows are read out of the production stencil: the coupling is the difference between the prescribed-value and the insulating wall cell. An end with zero conductance gets no node.

lmhdx.poisson.azimuthal_eigenvalues(grid: Grid) → ndarray¶

Return the eigenvalue of the azimuthal second difference for every mode.

The discrete Fourier basis diagonalizes a uniform periodic second difference, so mode m contributes -4 sin^2(pi m / N) / dtheta^2 divided by r^2. That last division is what stops the polar Laplacian separating into a sum of one-dimensional operators – and what makes it separate again once the azimuth is transformed, one radial operator per mode.

lmhdx.poisson.fast_diagonal_helmholtz(grid: Grid, offset: tuple[float, float, float], conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], *, shift: float = 1.0, coefficient: float = 1.0, precision: str = 'state', refinements: int = 1) → FastDiagonalHelmholtz¶

Factorize shift * I - coefficient * laplacian at one staggered position.

precision="mixed" solves float64 right-hand sides in float32 with refinements float64 corrections (module docstring).

lmhdx.poisson.fast_diagonal_poisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], *, precision: str = 'state', refinements: int = 2) → FastDiagonalPoisson¶

Factorize the separable Laplacian for grid under conditions.

precision="mixed" solves float64 right-hand sides in float32 with refinements float64 corrections (module docstring).

lmhdx.poisson.fast_diagonal_polar_poisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], *, shift: float = 0.0, coefficient: float = -1.0, precision: str = 'state', refinements: int = 2, wall_conductance: float = 0.0) → FastDiagonalPolarPoisson¶

Factorize shift*I - coefficient*laplacian, one eigendecomposition per azimuthal mode.

The default leaves the Laplacian itself. A positive shift with a positive coefficient is the damped operator that preconditions a pipe at large Hartmann number. precision="mixed" solves float64 right-hand sides in float32 with refinements float64 corrections (module docstring).

wall_conductance closes the outer radius of the potential Laplacian with a thin conducting wall: a sheet node on the radial operator (assemble_thin_wall_operator()) that sees the azimuthal eigenvalue at the wall radius, so each mode is still one eigendecomposition.

lmhdx.poisson.fast_diagonal_thin_wall_poisson(grid: Grid, conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition], conductance: tuple[float, float, float], *, precision: str = 'state', refinements: int = 2, layers: tuple | None = None) → FastDiagonalThinWallPoisson¶

Factorize the potential Laplacian with a thin wall of ratio conductance[axis] on both walls of an axis.

Zero keeps the insulating closure of conditions; a conducting axis must have insulating (Neumann) walls to replace, and one axis at least must not conduct. layers gives axes walls resolved in cells instead (assemble_thick_wall_operator()), (lower, upper) per axis or None, with no thin wall. Two axes may have them when the third has one cell: the Kronecker sum would give a corner cell the product of its two walls’ ratios, which _corner_fix() corrects exactly.

lmhdx.poisson.free_slice(grid: Grid, axis: int, offset_value: float, condition: BoundaryCondition) → slice¶

Return the entries of a staggered axis that a solve may change.

A cell-centred axis is free everywhere. A face-centred axis with a wall has its two boundary faces prescribed, and a periodic one carries a duplicate of its first face at the end. Solving on anything else would either invent a value for a prescribed face or treat one face as two unknowns.

Fully developed pipe flow¶

The pipe names remain importable from lmhdx.poisson, where they lived in 1.9.

Fully developed flow along a circular pipe in a transverse magnetic field.

A pipe is the other half of the duct validation, and the geometry the ALEX B1 experiment used. Nothing about the physics changes – the same inductionless system, the same insulating or thin conducting wall – but the coordinates do, and with them what separates and what does not.

Fully developed means the velocity is axial and depends only on the cross section, \(\mathbf u = u(r,\theta)\hat z\). With \(\mathbf B = B\hat x\) the Lorentz force is then purely axial as well, so the transverse momentum equations are satisfied by rest and the whole problem is one scalar equation coupled to the potential:

\[\nabla^2 u + B\,\partial_y\varphi - B^2 u + f = 0, \qquad \nabla\cdot\mathbf J = 0,\]

with \(\mathbf J = -\nabla\varphi + \mathbf u\times\mathbf B\). Both unknowns are cell centred, which is why this needs none of the staggered vector machinery of lmhdx.core3d: there is no cross flow to project.

The electromotive force is where the coordinates show. A uniform Cartesian field is not uniform in polar components, \(B_r = B\cos\theta\) and \(B_\theta = -B\sin\theta\), so \(\mathbf u\times\mathbf B\) has \(B u\sin\theta\) through a radial face and \(B u\cos\theta\) through an azimuthal one. The same face currents then carry the potential equation and the axial force, which is the consistency the whole package is built on.

The system is linear in the velocity, so it is one preconditioned Krylov solve, not a Newton iteration. The preconditioner is the damped operator \((B^2 - \nabla^2)^{-1}\), factorized exactly by lmhdx.poisson.fast_diagonal_polar_poisson(), which is what keeps the iteration count from growing with the Hartmann number.

class lmhdx.pipe.PipeProblem(grid: Grid, hartmann: float, wall_conductance: float = 0.0, forcing: float = 1.0)¶

A circular pipe of unit radius in a transverse field of magnitude hartmann.

property conditions: tuple[BoundaryCondition, BoundaryCondition, BoundaryCondition]¶: Conditions on the potential: insulating at the wall, periodic elsewhere.

factorization()¶: Factorize the potential Laplacian, with the thin wall’s sheet node when the wall conducts.

preconditioner()¶: Factorize the damped operator that preconditions the velocity solve.

lmhdx.pipe.flow_rate(velocity: Field) → float¶: Return the mean axial velocity over the cross-section.

lmhdx.pipe.pipe_grid(radial: int, azimuthal: int, hartmann: float, *, cells_in_layer: int = 6) → Grid¶: Return a polar grid whose radial cells resolve the 1/Ha wall layer.

lmhdx.pipe.pipe_problem(*, hartmann: float, radial: int = 48, azimuthal: int = 64, wall_conductance: float = 0.0) → PipeProblem¶: Build a pipe whose mesh resolves the layer its Hartmann number implies.

lmhdx.pipe.pipe_residual(velocity: Field, problem: PipeProblem, factorization) → Field¶: Return the steady axial momentum residual of a candidate velocity.

lmhdx.pipe.solve_pipe(problem: PipeProblem, *, tolerance: float = 1e-11, max_restarts: int = 40)¶

Return the axial velocity and the potential of a fully developed pipe.

The system is linear in the velocity, so this is one preconditioned Krylov solve. The preconditioner inverts the damped operator exactly, which is the stiff part of the problem and the reason the iteration does not lengthen with the Hartmann number.

Quasi-two-dimensional flow¶

Quasi-two-dimensional inductionless MHD on periodic planes.

class lmhdx.q2d.Q2DDiagnostics(kinetic_energy_initial: float, kinetic_energy_final: float, enstrophy_final: float, energy_budget_residual: float, max_divergence: float, max_courant: float)¶: Compact physical and numerical gates for a Q2D solve.

class lmhdx.q2d.Q2DProblem(initial_vorticity: Array, forcing: Array | None = None, length: tuple[float, float] = (6.283185307179586, 6.283185307179586), viscosity: float = 0.01, hartmann_friction: float = 0.1, dt: float = 0.01, steps: int = 100, history_stride: int = 0, adjoint_checkpoint_size: int | None = None, energy_budget_tolerance: float = 0.001)¶

Periodic Sommeria–Moreau vorticity problem.

forcing is a vorticity source. The linear Hartmann-layer closure is -hartmann_friction * vorticity; lengths, viscosity, time, and forcing must use one consistent unit system. All inputs determine a common real working dtype through JAX promotion, with at least float32 precision. energy_budget_tolerance bounds the normalized trapezoidal energy defect for host acceptance; it does not change the numerical trajectory.

class lmhdx.q2d.Q2DResult(problem: Q2DProblem, x: Array, y: Array, vorticity: Array, velocity_x: Array, velocity_y: Array, frame_times: Array, vorticity_history: Array, diagnostics: Q2DDiagnostics, status: str)¶

Final Q2D fields, optional frames, and acceptance diagnostics.

property converged: bool¶: Whether the finite trajectory passed acceptance, not steady convergence.

property fields: Q2DResult¶: Return the field-bearing result, matching the common solve contract.

lmhdx.q2d.evolve_q2d(initial_vorticity: Array, *, forcing: Array | None = None, length: tuple[float, float] = (6.283185307179586, 6.283185307179586), viscosity: float | Array = 0.01, hartmann_friction: float | Array = 0.1, dt: float | Array = 0.01, steps: int = 100, adjoint_checkpoint_size: int | None = None) → tuple[Array, Array, Array]¶

Return final Q2D fields through a JIT- and autodiff-safe numerical core.

Array state, forcing, domain lengths, viscosity, Hartmann friction, and timestep are differentiable. steps and adjoint_checkpoint_size are static controls; the default checkpoint schedule retains O(sqrt(steps)) trajectory states in reverse mode. State, forcing and coefficients use JAX’s common real dtype, with at least float32 precision; float64 requires enabling JAX x64 before constructing inputs.

lmhdx.q2d.make_q2d_case(*, shape: tuple[int, int] = (64, 64), length: tuple[float, float] = (6.283185307179586, 6.283185307179586), mode: tuple[int, int] = (1, 1), amplitude: float = 1.0, viscosity: float = 0.01, hartmann_friction: float = 0.1, dt: float = 0.01, steps: int = 100, history_stride: int = 0, energy_budget_tolerance: float = 0.001) → Q2DProblem¶: Build a Taylor–Green decay case with an analytical exponential rate.

lmhdx.q2d.solve_q2d(problem: Q2DProblem) → Q2DResult¶: Evolve IFRK4 fields and enforce Courant and normalized energy-budget gates.

Fully developed cases on the staggered core¶

Fully developed duct cases solved on the staggered core.

A CaseSpec for a rectangular duct becomes a ChannelProblem with one periodic axial cell, solved by the conjugate-gradient steady solve of lmhdx.steady and reported on the case’s cross-section. That is the route lmhdx.solve() and lmhdx.solve_fully_developed_fields() take.

Mesh. ny and nz are the fluid cells along y and z. The faces follow lmhdx.core3d.duct_problem(): the walls normal to the field are clustered to the Hartmann layer delta = sqrt(rho nu / sigma) / |B|, the others to the side layer sqrt(a delta), with a the half-width along the field, six cells in each layer (hartmann_layer_cells overrides) and the gentlest stretching that spans the duct; an odd count adds one centre cell. A 2 x 2 duct with unit properties gets exactly the faces of duct_problem(hartmann=Ha, cells=n).

Walls. An insulating wall is the homogeneous Neumann closure. The conducting walls of a layered_duct follow geometry.wall_model. "thin" makes each a sheet of conductance sigma_w t_w / sigma (lmhdx.poisson) with no cells of its own, which needs equal walls on one axis. "resolved" gives each wall wall_cells uniform cells of its own conductivity, insulated outside: one wall, two different walls, or a layer that no boundary names (ChannelProblem.wall_layers); a corner cell takes the nearer wall’s material, as the retired cell-centred solver assigned it, and the reported fields cover the fluid. The default "auto" is thin where that holds and resolved otherwise. A wall stack of several materials is a ChannelProblem with per-cell ratios.

Field. A constant field is three numbers; an analytic or tabulated one is sampled at the cell centres as an ImposedField, and the layers follow its peak transverse strength. An axial component is kept: it adds no electromotive force to the axial flow, and on a varying field it can drive a secondary flow, which the retired cell-centred solver dropped.

Drive. forcing is the axial force density. With zero forcing and an inlet_flow_rate boundary, the flow rate is met by scaling the unit-drive solution, which is exact because the problem is linear in the drive.

The case’s time stepper does not enter a steady solve: the steady state is one preconditioned CG solve to a relative residual of 1e-9. Without lmhdx.enable_x64() it runs in float32 to 1e-5, which the solve reaches at Ha 20 on 32 cells and not at Ha 100 on 48, where it raises; a float32 case with float64 enabled is solved in float64 and returned in float32.

Fully developed duct design: throughput, pumping power and their derivatives.

At fixed field, materials and geometry, flow is linear in force density: Q = G f. One unit-drive solve gives G, eliminating the drive from fixed-flow optimization as f = Q_target / G with df/dQ = 1/G. For a uniform pressure-gradient drive over length L, pressure drop is f L and hydraulic power is f L Q. These are isothermal segment quantities, excluding entry/exit losses, manifolds and thermal effects.

A CaseSpec is solved on the staggered core through lmhdx.fully_developed, so both families share one solver. The channel_* functions are the ChannelProblem-native counterparts, sharing the same linear response, reusing DuctResponse, pressure_drop() and hydraulic_power(), and the same segment-quantity disclaimer above. They apply in the Stokes limit only: Q = G f holds because the steady residual is affine in the drive when advection is "off", so channel_flow_response() rejects any other value.

class lmhdx.fully_developed.DuctResponse(flow_per_unit_drive: Array, magnetic_field_scale: Array)¶

The linear throughput response of one duct at one field strength.

drive_for(target_flow_rate: float | Array) → Array¶: Return the drive that delivers target_flow_rate exactly.

lmhdx.fully_developed.case_mesh(case: CaseSpec) → StructuredMesh¶: Return the fluid cross-section the core solves case on, as a StructuredMesh.

lmhdx.fully_developed.channel_cross_section_weights(problem: ChannelProblem) → Array¶

Return the cross-section integration weights of a ChannelProblem.

Axis 0 is the flow axis of every channel this package builds (see lmhdx.core3d.duct_problem()), so a cell’s weight is its transverse (y, z) area alone, independent of the axial spacing – unlike fluid_cell_areas(), a channel carries no fluid mask, so every transverse cell counts.

lmhdx.fully_developed.channel_drive_for_flow_rate(problem: ChannelProblem, target_flow_rate: float | Array, *, magnetic_field_scale: float | Array = 1.0) → Array¶: Return the axial force density whose fully developed flow rate is target_flow_rate.

lmhdx.fully_developed.channel_fixed_flow_hydraulic_power(problem: ChannelProblem, target_flow_rate: float, length: float, *, magnetic_field_scale: float | Array = 1.0) → Array¶

Return the hydraulic power needed to hold a throughput at a given field, on a ChannelProblem.

The drive is eliminated analytically, so this is a function of the design inputs alone and is differentiable through the solve.

lmhdx.fully_developed.channel_flow_rate(problem: ChannelProblem, velocity: Array) → Array¶

Integrate an axial velocity slice over a ChannelProblem cross-section.

velocity is the flow-axis component at one axial station – every station carries the same value by periodicity – shaped like channel_cross_section_weights().

lmhdx.fully_developed.channel_flow_response(problem: ChannelProblem, *, magnetic_field_scale: float | Array = 1.0) → DuctResponse¶

Measure G = Q(f = 1) on a ChannelProblem, the flow a unit axial drive produces.

One solve determines the whole drive-to-flow relation because the Stokes residual is affine in the drive; ChannelProblem.advection must be "off", since otherwise the residual carries -div(uu) and Q = G f breaks down.

lmhdx.fully_developed.channel_problem(case: CaseSpec) → ChannelProblem¶

Return the staggered-core problem a fully developed case solves, at unit drive.

Raise ValueError for what a case cannot mean: fluids of different properties (a CaseSpec gives regions no geometry), an imposed current density, unequal thin walls, a conducting wall without a layered duct. A constant field stays three numbers; an analytic or tabulated one is sampled at the cell centres as an ImposedField.

lmhdx.fully_developed.drive_for_flow_rate(case: CaseSpec, target_flow_rate: float | Array, *, magnetic_field_scale: float | Array = 1.0) → Array¶: Return the force density whose fully developed flow rate is target_flow_rate.

lmhdx.fully_developed.fluid_cell_areas(case: CaseSpec) → Array¶

Return the cross-section weights of the fluid mesh a case is solved on.

That mesh is lmhdx.fully_developed.case_mesh(), the one the velocity of lmhdx.solve_fully_developed_fields() and lmhdx.solve() lives on.

lmhdx.fully_developed.hydraulic_power(drive: float | Array, flow_rate: float | Array, length: float) → Array¶

Return the isothermal hydraulic power of a fully developed segment.

This is the pressure drop times the throughput. It excludes entry and exit losses, manifolds and every thermal effect, so it is a segment quantity and not a blanket pumping budget.

lmhdx.fully_developed.linear_flow_response(case: CaseSpec, *, magnetic_field_scale: float | Array = 1.0) → DuctResponse¶

Measure G = Q(f = 1), the flow rate a unit drive produces.

One solve determines the whole drive-to-flow relation because the problem is linear in the drive.

lmhdx.fully_developed.pressure_drop(drive: float | Array, length: float) → Array¶: Return the pressure drop of a uniform pressure-gradient drive over length.

lmhdx.fully_developed.solve_fully_developed(case: CaseSpec, *, logger=None, start_time: float = 0.0) → Solution¶

Solve a fully developed case to its steady state on the core and report it as a Solution.

status is "converged" with residual the relative steady residual ||R(u)|| / ||R(0)||; a solve that fails raises. steps is zero: the solve is one CG call with no outer iterations. The diagnostics hold one record.

lmhdx.fully_developed.solve_fully_developed_fields(case: CaseSpec, *, forcing: float | Array | None = None, magnetic_field_scale: float | Array = 1.0) → tuple[Array, Array, Array, Array, Array]¶

Return the steady velocity, potential, currents and Lorentz force on case_mesh().

All five are cell-centred (ny, nz) arrays in the case’s dtype: axial velocity, potential (zero volume mean), the y and z currents averaged from their faces, and the axial Lorentz force density. forcing and magnetic_field_scale are continuous design inputs, differentiable through the implicit solve of lmhdx.steady.solve_steady_state(); the case itself is static, so close over it outside jax.jit(). A solve that fails raises NumericalFailure when called with concrete inputs and gives nonfinite fields under tracing.

lmhdx.fully_developed.solve_fully_developed_transient(case: CaseSpec, logger=None, *, initial_state: MHDState | None = None, initial_diagnostics: Diagnostics | None = None, append_diagnostics: bool = False, restart_info=None) → Solution¶

Run a transient fully developed case on the core, by implicit Euler steps.

Each step is (u - u_n)/dt = A u + b for the core’s Stokes-limit operator A and drive b, solved by the steady solve’s CG with a mass term and its preconditioner factorized at the step (see _transient_programs()). The steps run at time_stepper.dt from initial_state (or rest plus initial_velocity) to t_final, at most max_steps, compiled as one scan between kept records. A ramped field scales the Lorentz force step by step. With an inlet_flow_rate and zero forcing each step meets the flow rate exactly. output.history_stride keeps every stride-th step and the last (0, the last alone); residual_history is the step’s largest velocity change, linear_iterations_history its CG iterations, and status is "completed". A step whose CG fails gives nonfinite fields, which raise.

lmhdx.fully_developed.volumetric_flow_rate(case: CaseSpec, velocity: Array) → Array¶: Integrate an axial velocity over the fluid cross-section.

Outputs, validation reports and the command line¶

Outputs, evidence reports and the command line.

A solved case is written as ParaView files, midplane and centreline profiles, NPZ archives, restart bundles and figures; the analytical, conservation and benchmark reports (Hartmann acceptance, profile metrics, solver benchmarks) read the same solutions; and lmhdx (main()) runs a named or TOML case, a validation or a benchmark from the command line.

class lmhdx.io.AcceptanceReport(case_name: 'str', l2_error: 'float', linf_error: 'float', l2_threshold: 'float', linf_threshold: 'float', passed_l2: 'bool', passed_linf: 'bool', passed: 'bool')¶

class lmhdx.io.AnalyticComparison(coordinate: 'jnp.ndarray', simulated: 'jnp.ndarray', reference: 'jnp.ndarray', l2_error: 'float', linf_error: 'float')¶

class lmhdx.io.ProfileSymmetry(axis: 'str', mean_abs_error: 'float', max_abs_error: 'float')¶

class lmhdx.io.RestartBundle(path: 'Path', state: 'MHDState', diagnostics: 'Diagnostics', metadata: 'dict[str, object]', y_faces: 'np.ndarray', z_faces: 'np.ndarray', geometry_kind: 'str')¶