API Reference

benchmarks

class benchmarks.LidDrivenCavity3D

Bases: benchmark

Three-dimensional lid-driven cavity benchmark.

This benchmark sets up a Navier–Stokes problem on a unit cube with a smooth lid-driven boundary condition and graded mesh refinement.

generate_mesh(Lx=1.0, Ly=1.0, Lz=1.0, lc=None, out='./geo/mesh/box.msh')

Generate a 3D tetrahedral mesh of a rectangular box using Gmsh.

The mesh is created using OpenCASCADE (OCC) geometry and includes:

  • A single volume representing the domain.

  • Physical groups for all boundary faces (xmin, xmax, ymin, ymax, zmin, zmax).

  • Graded refinement near: * the top face (z = Lz), * the face at x = Lx.

  • Background mesh size field based on a distance-threshold strategy.

  • MFEM-compatible MSH 2.2 ASCII output.

Parameters:

  • Lx (float) – Length of the domain in the x-direction.

  • Ly (float) – Length of the domain in the y-direction.

  • Lz (float) – Length of the domain in the z-direction.

  • lc (float or None) – Base mesh size. If None, defaults to min(Lx, Ly, Lz) / 10.

  • out (str) – Output path for the generated .msh file.

Raises:

RuntimeError – If the volume or required boundary faces cannot be identified.

Returns:

None

Return type:

None

class benchmarks.benchmark(name, SimulationHelper, SimulationDataProcessor, dts, refinements, orders)

Bases: object

Harness for running benchmark tests in a reproducible manner.

This class wires together:

plot_euler()

Pull Euler results and plot them using the local pipeline.

plot_local()

Collect and plot results from local runs.

run_euler()

Run the full parameter sweep using the Euler backend.

run_local()

Run the full parameter sweep locally.

SimIO

class SimIO.ExactManipulationsSpace(coords)

Bases: object

Symbolic differential operators for vector calculus in 3D.

This class provides symbolic implementations of common vector calculus operators using sympy, including gradient, divergence, curl, Laplacian, and convective terms.

Parameters:

coords (list[sympy.Symbol]) – Spatial coordinate symbols (e.g. [x, y, z]).

convective_term(u)

Compute the convective term using the identity

\[(\mathbf{u} \cdot \nabla) \mathbf{u} = -\mathbf{u} \times (\nabla \times \mathbf{u})\]
Parameters:

u (sympy.Matrix) – Velocity vector field.

Returns:

Convective term.

Return type:

sympy.Matrix

curl(u)

Compute the 3D curl of a vector field.

Assumes coordinates are ordered as [x, y, z].

Parameters:

u (sympy.Matrix) – Vector field Matrix([u1, u2, u3]).

Returns:

Curl vector.

Return type:

sympy.Matrix

divergence(u)

Compute the divergence of a vector field.

Parameters:

u (sympy.Matrix) – Vector field.

Returns:

Scalar divergence.

Return type:

sympy.Expr

grad(scalar_field)

Compute the gradient of a scalar field.

Parameters:

scalar_field (sympy.Expr) – Scalar symbolic expression.

Returns:

Gradient vector.

Return type:

sympy.Matrix

laplacian_vector(vec_field)

Compute the component-wise Laplacian of a vector field.

Parameters:

vec_field (sympy.Matrix) – Vector field.

Returns:

Vector Laplacian.

Return type:

sympy.Matrix

class SimIO.ExactManipulationsSpaceTime(coords, t)

Bases: ExactManipulationsSpace

Extension of ExactManipulationsSpace to space–time settings.

This class adds symbolic time differentiation to the spatial vector calculus operators provided by the base class.

time_derivative(u, t)

Compute the time derivative of a symbolic expression.

\[\frac{\partial u}{\partial t}\]
Parameters:

  • u (sympy.Expr or sympy.Matrix) – Symbolic expression or vector field.

  • t (sympy.Symbol) – Time variable.

Returns:

Time derivative of u.

Return type:

sympy.Expr or sympy.Matrix

class SimIO.IBVPNavierStokes(u_init: MutableDenseMatrix, nu: float, coords, t, u_boundary=Matrix([[0], [0], [0]]), force=Matrix([[0], [0], [0]]))

Bases: ExactManipulationsSpaceTime

Symbolic representation of an initial-boundary value problem (IBVP) for the incompressible Navier–Stokes equations in 3D.

The problem has the form

\[\partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} - \nu \Delta \mathbf{u} + \nabla p = \mathbf{f},\]

with prescribed initial condition

\[\mathbf{u}(\cdot, 0) = \mathbf{u}_0,\]

and boundary condition

\[\mathbf{u}|_{\partial \Omega} = \mathbf{u}_\text{boundary}.\]
Parameters:

  • u_init (sympy.Matrix) – Initial velocity field (3×1 vector).

  • nu (float) – Kinematic viscosity.

  • coords (list[sympy.Symbol]) – Spatial coordinate symbols (e.g. [x, y, z]).

  • t (sympy.Symbol) – Time symbol.

  • u_boundary (sympy.Matrix) – Prescribed boundary velocity. Defaults to zero.

  • force (sympy.Matrix) – External forcing term. Defaults to zero.

Raises:

ValueError – If u_init is not a 3×1 vector.

get_force()

Return the external forcing term.

Returns:

Forcing vector.

Return type:

sympy.Matrix

get_u_boundary()

Return the prescribed boundary velocity.

Returns:

Boundary velocity.

Return type:

sympy.Matrix

get_u_init()

Return the initial velocity field.

Returns:

Initial velocity.

Return type:

sympy.Matrix

get_viscosity()

Return the viscosity parameter.

Returns:

Kinematic viscosity.

Return type:

float

get_vorticity_init()

Return the initial vorticity field.

\[\boldsymbol{\omega}_0 = \nabla \times \mathbf{u}_0\]
Returns:

Initial vorticity.

Return type:

sympy.Matrix

class SimIO.IBVPNavierStokesSolution(u: MutableDenseMatrix, p: Expr, nu: float, coords, t, u_boundary=Matrix([[0], [0], [0]]), force=Matrix([[0], [0], [0]]))

Bases: IBVPNavierStokes

Symbolic exact solution of the incompressible Navier–Stokes IBVP.

This class represents a fully specified velocity–pressure pair \((\mathbf{u}, p)\) depending on space and time. The initial condition is automatically extracted as

\[\mathbf{u}_0 = \mathbf{u}(\cdot, 0).\]

The class inherits the Navier–Stokes IBVP structure from IBVPNavierStokes.

Parameters:

  • u (sympy.Matrix) – Time-dependent velocity field.

  • p (sympy.Expr) – Time-dependent pressure field.

  • nu (float) – Kinematic viscosity.

  • coords (list[sympy.Symbol]) – Spatial coordinate symbols (e.g. [x, y, z]).

  • t (sympy.Symbol) – Time symbol.

  • u_boundary (sympy.Matrix) – Prescribed boundary velocity. Defaults to zero.

  • force (sympy.Matrix) – External forcing term. Defaults to zero.

get_p()

Return the pressure field.

Returns:

Pressure.

Return type:

sympy.Expr

get_u()

Return the velocity field.

Returns:

Velocity field.

Return type:

sympy.Matrix

get_vorticity()

Compute the vorticity of the velocity field.

\[\boldsymbol{\omega} = \nabla \times \mathbf{u}\]
Returns:

Vorticity vector.

Return type:

sympy.Matrix

class SimIO.NavierStokesSimulationHelper(exact_equations: IBVPNavierStokes, name: str, mesh='./extern/mfem/data/ref-cube.mesh', visualisation=1, printlevel=0)

Bases: SimulationHelper

Simulation helper specialized for Navier–Stokes IBVP data.

This subclass converts symbolic Navier–Stokes problem data (initial/boundary conditions, forcing, and optional exact solutions) into a JSON configuration dictionary understood by the external solver.

Parameters:

  • exact_equations (IBVPNavierStokes) – Navier–Stokes IBVP data container providing viscosity, forcing, initial conditions, and optionally boundary/exact data.

  • name (str) – Simulation/benchmark name used for config and output folder names.

  • mesh (str) – Path to the MFEM mesh file.

  • visualisation (int) – Visualisation flag forwarded to the solver.

  • printlevel (int) – Verbosity level forwarded to the solver.

base_config_file()

Build the base JSON configuration dictionary for the external solver.

The returned dictionary always includes:

  • mesh path, solver type, visualisation, printlevel

  • viscosity

  • forcing term

  • initial velocity and initial vorticity

Additionally, the following fields are included when available:

  • boundary_data_u if boundary velocity is provided

  • exact_data_u and exact_data_w if exact solution fields are provided

Returns:

Configuration dictionary to be written as JSON.

Return type:

dict

class SimIO.SimulationDataProcessor(name: str)

Bases: object

Collect and post-process simulation output data.

This class scans solver output files in out/data/<name>/ and produces convergence plots (error vs refinement) and optional conservation plots. It can also pull result files from the Euler cluster via SFTP.

Parameters:

name (str) – Simulation/benchmark name used to locate input/output folders.

add_reference_triangles(ax, order, x_anchor, y_anchor)

Add a dotted reference line indicating expected convergence order.

The reference indicator is anchored at (x_anchor, y_anchor) and assumes that one refinement step corresponds to halving the mesh size.

Parameters:

  • ax (matplotlib.axes.Axes) – Matplotlib axes to draw on.

  • order (int) – Expected convergence order (used as exponent in \(O(h^p)\)).

  • x_anchor (float) – Anchor x-location (typically the last refinement index).

  • y_anchor (float) – Anchor y-location (typically the last error value).

Returns:

None

Return type:

None

collect_data()

Collect error data from solver output CSV files.

This scans files matching:

out/data/<name>/<name>_conv_order*_ref*_vars.csv

Each file is parsed and the last row is used (no time matching). Errors are read from columns listed in error_columns.

The returned structure is grouped by polynomial order:

  • data[order]["refs"]: 1D array of refinement levels

  • data[order]["errs"][col]: 1D array of error values for column col

Returns:

Collected data grouped by order.

Return type:

dict[int, dict[str, object]]

Raises:

ValueError – If any required error column is missing in a CSV file.

plot_conserved_variables()

Plot conserved/diagnostic quantities over time from solver CSV output.

This scans the same output files as collect_data() and, for each column listed in cons_columns, produces a time series plot saved to:

out/plots/<name>/conservation/<name>_<column>_order<k>_ref<r>.png

Requires that the CSV files contain a time_full column.

Returns:

None

Return type:

None

Raises:

ValueError – If time_full or any column in cons_columns is missing in a CSV file.

plot_convergence(show_plot=False, reference_order=<function SimulationDataProcessor.<lambda>>)

Plot convergence curves (log-scale error vs refinement index).

Requires collect_data() to have been called first.

One figure is generated per error column in error_columns and saved to:

out/plots/<name>/convergence/<name>_<error_column>
Parameters:

  • show_plot (bool) – If True, display the plot window; otherwise close after saving.

  • reference_order (collections.abc.Callable[[int], int]) – Function mapping polynomial order to reference slope order, used for the \(O(h^p)\) indicator. Signature: reference_order(order: int) -> int.

Returns:

None

Return type:

None

pull_data_from_euler()

Pull solver output CSV files from the Euler cluster via SFTP.

For each config file found in:

data/config/<name>/

this method attempts to download the corresponding output file:

out/data/<name>/<config_basename>_vars.csv

from the remote path under /cluster/home/<user>/dualfieldmfem/ into the local out/data/<name>/ directory.

Missing files are skipped with a warning.

Returns:

None

Return type:

None

class SimIO.SimulationHelper(name: str)

Bases: object

Base helper class for running simulation sweeps and managing configuration I/O.

This class provides utilities to: - convert SymPy expressions to C-style strings for code generation, - generate per-(order, refinement) JSON config files, - run all generated configs locally via an external executable, - submit/run configs on the Euler cluster via SSH.

Subclasses must implement base_config_file().

Parameters:

name (str) – Simulation/benchmark name used for config and output folder names.

base_config_file() → dict

Return a base configuration dictionary used by generate_config_files().

Subclasses must override this method.

Raises:

NotImplementedError – Always, unless overridden in a subclass.

Returns:

Base configuration dictionary.

Return type:

dict

expr_to_c(expr)

Convert a SymPy expression to a compact C expression string.

The conversion uses sympy.printing.ccode.ccode() and then rewrites coordinate symbols x0, x1, x2 to array access x[0], x[1], x[2].

Parameters:

expr (sympy.Expr) – SymPy expression to convert.

Returns:

C expression string with x[i] indexing and no spaces.

Return type:

str

generate_config_files(T: float, dt: Callable[[int, int], float], refinements: Callable[[int], Iterable[int]], orders: List[int], tol=1e-08)

Generate JSON configuration files for a convergence/sweep study.

Config files are written to ./data/config/<name>/. Any existing directory at that location is removed and recreated.

For each order in orders and each refinement in refinements(order), this method writes a file named:

<name>_conv_order<order>_ref<refinement>.json
Parameters:

  • T (float) – Final simulation time.

  • dt (collections.abc.Callable[[int, int], float]) – Timestep selection rule. Signature: dt(order: int, refinement: int) -> float.

  • refinements (collections.abc.Callable[[int], collections.abc.Iterable[int]]) – Refinement selection rule. Signature: refinements(order: int) -> Iterable[int].

  • orders (list[int]) – Polynomial orders to run.

  • tol (float) – Solver tolerance written to the config file.

Returns:

None

Return type:

None

run_all_configs()

Run all generated configuration files locally using an external executable.

Expects configuration files in:

./data/config/<name>/

Writes outputs to:

./out/data/<name>/

Uses the executable:

./build/MEHCscheme
Raises:

RuntimeError – If the executable cannot be found or executed.

Returns:

None

Return type:

None

run_all_configs_euler(time='24:00:00', mempercpu='512G', cpuspertask='1')

Submit all generated configuration files to the Euler cluster via SSH.

This method: - connects to euler.ethz.ch using an SSH key, - updates/builds the project on the cluster, - uploads config files to the remote config directory, - submits one SLURM job per config via sbatch.

Parameters:

  • time (str) – SLURM wall time limit (e.g. "24:00:00").

  • mempercpu (str) – SLURM memory per CPU (e.g. "512G").

  • cpuspertask (str) – SLURM CPUs per task (e.g. "1").

Returns:

None

Return type:

None

Raises:

RuntimeError – If remote commands fail (via internal helpers).

run_convergence(T: float, dt: Callable[[int, int], float], refinements: Callable[[int], Iterable[int]], orders: List[int], tol=1e-08)

Generate configs and run all cases locally.

This is a convenience wrapper calling generate_config_files() followed by run_all_configs().

Parameters:

  • T (float) – Final simulation time.

  • dt (collections.abc.Callable[[int, int], float]) – Timestep rule dt(order, refinement) -> float.

  • refinements (collections.abc.Callable[[int], collections.abc.Iterable[int]]) – Refinement rule refinements(order) -> Iterable[int].

  • orders (list[int]) – Polynomial orders to run.

  • tol (float) – Solver tolerance written to the config file.

Returns:

None

Return type:

None

run_remote_command(client: SSHClient, command: str) → str

Run a command on a connected SSH client and return stdout.

Parameters:

  • client (paramiko.SSHClient) – Connected SSH client.

  • command (str) – Command string to run remotely.

Returns:

Command stdout (decoded).

Return type:

str

Raises:

RuntimeError – If the command cannot be started, the connection drops, or the remote command returns a non-zero exit code.

sp_vector_to_str(expr)

Convert a SymPy vector expression into C assignment statements.

Produces a single string of the form:

out[0] = <expr0>; out[1] = <expr1>; out[2] = <expr2>;
Parameters:

expr (sympy.Matrix) – Vector-valued SymPy expression.

Returns:

Concatenated C assignment statements.

Return type:

str

class SimIO.StokesSimulationHelper(exact_equations: manufacturedStokes, name: str)

Bases: SimulationHelper

Simulation helper specialized for manufactured Stokes problems.

This subclass converts symbolic manufactured Stokes data into a JSON configuration dictionary understood by the external solver.

Parameters:

  • exact_equations (manufacturedStokes) – Manufactured Stokes data container providing mass, viscosity, forcing, and exact velocity.

  • name (str) – Simulation/benchmark name used for config and output folder names.

base_config_file()

Build the base JSON configuration dictionary for the external Stokes solver.

The returned dictionary includes:

  • mesh path (default MFEM reference cube)

  • solver type, visualisation, and printlevel

  • mass and viscosity parameters

  • forcing term (right-hand side)

  • exact velocity solution

Returns:

Configuration dictionary to be written as JSON.

Return type:

dict

class SimIO.manufacturedNavierStokes(u: MutableDenseMatrix, p: Expr, nu: float, coords, t)

Bases: IBVPNavierStokesSolution

Manufactured solution for the incompressible Navier–Stokes equations.

Given symbolic velocity and pressure fields \((\mathbf{u}, p)\), this class constructs the corresponding forcing term so that \((\mathbf{u}, p)\) is an exact solution of

\[\partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p - \nu \Delta \mathbf{u} = \mathbf{f}.\]

The right-hand side is computed symbolically.

Parameters:

  • u (sympy.Matrix) – Time-dependent velocity field.

  • p (sympy.Expr) – Time-dependent pressure field.

  • nu (float) – Kinematic viscosity.

  • coords (list[sympy.Symbol]) – Spatial coordinate symbols.

  • t (sympy.Symbol) – Time symbol.

navier_stokes_rhs(u, p, nu, t)

Compute the Navier–Stokes forcing term corresponding to a manufactured solution.

The forcing is constructed as

\[\mathbf{f} = \partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p - \nu \Delta \mathbf{u}.\]
Parameters:

  • u (sympy.Matrix) – Velocity field.

  • p (sympy.Expr) – Pressure field.

  • nu (float) – Kinematic viscosity.

  • t (sympy.Symbol) – Time symbol.

Returns:

Symbolic forcing term.

Return type:

sympy.Matrix

class SimIO.manufacturedStokes(u, p, mass, nu, coords)

Bases: ExactManipulationsSpace

Manufactured solution for the steady Stokes problem.

This class represents a symbolic manufactured solution of the form

\[m \mathbf{u} - \nu \Delta \mathbf{u} + \nabla p = \mathbf{f},\]

using the vector identity

\[\Delta \mathbf{u} = - \nabla \times (\nabla \times \mathbf{u}) \quad \text{(for divergence-free fields)}.\]

The right-hand side is constructed symbolically from the provided velocity and pressure fields.

Parameters:

  • u (sympy.Matrix) – Prescribed velocity field.

  • p (sympy.Expr) – Prescribed pressure field.

  • mass (sympy.Expr or float) – Mass coefficient (reaction term).

  • nu (sympy.Expr or float) – Kinematic viscosity.

  • coords (list[sympy.Symbol]) – Spatial coordinate symbols (e.g. [x, y, z]).

get_mass()

Return the mass coefficient.

Returns:

Mass coefficient.

Return type:

sympy.Expr or float

get_rhs()

Compute the forcing term corresponding to the manufactured solution.

The right-hand side is constructed as

\[\mathbf{f} = m \mathbf{u} + \nu \nabla \times (\nabla \times \mathbf{u}) + \nabla p.\]
Returns:

Symbolic right-hand side vector.

Return type:

sympy.Matrix

get_u()

Return the prescribed velocity field.

Returns:

Velocity field.

Return type:

sympy.Matrix

get_viscosity()

Return the viscosity parameter.

Returns:

Kinematic viscosity.

Return type:

sympy.Expr or float

get_vorticity()

Compute the vorticity field.

\[\boldsymbol{\omega} = \nabla \times \mathbf{u}\]
Returns:

Vorticity vector.

Return type:

sympy.Matrix