Diese Seite wurde maschinell aus dem Englischen übersetzt und kann Fehler enthalten. Englisches Original

OpenFOAM Backend

axqua.solvers.openfoam erstellt, führt und bewertet OpenFOAM-Modelle. aXqua schreibt die Rechennetz-Dateien selbst, denn eine Flusssohle ist eine Fläche mit einer Höhe pro Ort, der ein Rechennetz von Zellspalten direkt folgen kann.

Rechennetz

A self-contained OpenFOAM constant/polyMesh writer.

This is the OpenFOAM counterpart of axqua.core.selafin: axqua writes the mesh itself rather than driving a mesher, so the mesh is exactly what the case config asks for and nothing is left to a black box.

Why write the mesh directly instead of using snappyHexMesh

snappyHexMesh castellates a background block against an STL and then snaps the cut cells onto it. On a river bed - a single-valued, gently varying surface with no thin features - that machinery buys nothing and costs a great deal: the snapped layer is where the mesh degenerates (slivers, failed layer addition, negative volumes), and every fix is a search through snappyHexMeshDict for the knob that happens to help. A river bed does not need to be snapped to; it can be followed, because it is a height field z_b(x, y). Extruding a plan grid between the bed and a lid gives a mesh that is:

  • all-hexahedral - every cell is a hex, no polyhedra, no cut cells;

  • structured - a plain (i, j, k) grid (with the columns outside the domain simply absent), so the addressing is deterministic and the cell count is known before anything is built;

  • conformal by construction - the vertical levels are computed at the plan vertices, so neighboring columns share their corner points exactly and no matching or merging step is needed;

  • boundary-fitted at the bed - the bottom face of the lowest cell is the DEM surface, so the wall functions see the real bed, not a staircase.

The price is non-orthogonality wherever the bed is steep, which is measured (axqua.solvers.openfoam.quality) rather than assumed away.

Format notes

polyMesh requires, and this writer guarantees:

  • internal faces first, ordered upper-triangular: ascending owner, and within one owner ascending neighbour; owner < neighbour on every one;

  • boundary faces after them, contiguous per patch, in the boundary file order;

  • every face’s normal (right-hand rule over its point list) pointing from its owner cell to its neighbour - out of the domain for a boundary face.

Files are written in ASCII. The mesh is written once and read once; the run’s own output is binary (writeFormat binary in controlDict), and decomposePar writes the per-processor meshes in that format too, so ASCII here costs one slow read and buys a mesh that can be inspected and diffed.

axqua.solvers.openfoam.polymesh.foam_header(cls, obj, *, location=None, fmt='ascii')[Quellcode]

Render the FoamFile banner every OpenFOAM file starts with.

Parameter:
Rückgabetyp:

str

Rückgabetyp:

str

class axqua.solvers.openfoam.polymesh.Patch(name, type='patch', n_faces=0, start_face=0, in_groups=())[Quellcode]

Bases: object

One boundary patch: a contiguous run of faces in the global face list.

Parameter:
name: str
type: str = 'patch'
n_faces: int = 0
start_face: int = 0
in_groups: tuple[str, ...] = ()
render()[Quellcode]
Rückgabetyp:

str

class axqua.solvers.openfoam.polymesh.PolyMesh(points, faces, owner, neighbour, patches=<factory>, n_cells=0)[Quellcode]

Bases: object

An assembled OpenFOAM mesh, ready to write.

faces is an (NF, 4) integer array because every face of an all-hex mesh is a quadrilateral; the writer emits the general 4(a b c d) face syntax, so OpenFOAM neither knows nor cares that they are all quads.

Parameter:
points: numpy.ndarray
faces: numpy.ndarray
owner: numpy.ndarray
neighbour: numpy.ndarray
patches: list[Patch]
n_cells: int = 0
property n_points: int
property n_faces: int
property n_internal_faces: int
patch(name)[Quellcode]
Parameter:

name (str)

Rückgabetyp:

Patch

patch_face_ids(name)[Quellcode]

Global face indices belonging to name (the order the BC list must use).

Parameter:

name (str)

Rückgabetyp:

numpy.ndarray

face_centres(ids=None)[Quellcode]
Parameter:

ids (np.ndarray | None)

Rückgabetyp:

np.ndarray

face_normals(ids=None)[Quellcode]

Area-weighted face normals (Newell over the quad).

Parameter:

ids (np.ndarray | None)

Rückgabetyp:

np.ndarray

face_areas(ids=None)[Quellcode]
Parameter:

ids (np.ndarray | None)

Rückgabetyp:

np.ndarray

write(case_dir)[Quellcode]
Parameter:

case_dir (str | Path)

Rückgabetyp:

Path

axqua.solvers.openfoam.polymesh.assemble(points, internal, boundary, cell_centres)[Quellcode]

Order and orient raw face lists into a valid PolyMesh.

internal is (owner, neighbour, quads) for the internal faces in any order; boundary is a list of (patch_name, patch_type, owner, quads). cell_centres supplies the owner-to-neighbour direction used to orient every face - the structured builder knows them exactly, so orientation needs no topological reasoning.

Faces whose owner exceeds their neighbour are swapped (both the pair and the point order) rather than rejected, so a caller may emit either sense.

Parameter:
Rückgabetyp:

PolyMesh

axqua.solvers.openfoam.polymesh.write_polymesh(mesh, case_dir)[Quellcode]

Write <case_dir>/constant/polyMesh/{points,faces,owner,neighbour,boundary}.

Parameter:
Rückgabetyp:

Path

axqua.solvers.openfoam.polymesh.read_patch_names(case_dir)[Quellcode]

Patch names from a written constant/polyMesh/boundary, in file order.

Lets the dictionaries be regenerated against an existing mesh without rebuilding it - which is what keeps a run honoring the config it was launched with rather than the one it was built with.

Parameter:

case_dir (str | Path)

Rückgabetyp:

list[str]

axqua.solvers.openfoam.polymesh.validate(mesh)[Quellcode]

Return the structural problems that would make OpenFOAM reject the mesh.

Cheap, exact checks on the addressing only - geometric quality (orthogonality, skewness, volume) is axqua.solvers.openfoam.quality’s job.

Parameter:

mesh (PolyMesh)

Rückgabetyp:

list[str]

Terrain-following, all-hexahedral mesh generation for OpenFOAM river cases.

The mesh is a structured (i, j, k) grid, warped onto the bed:

  • in plan, a uniform Cartesian lattice of spacing dx, optionally rotated so its axes line up with the reach (align_to_flow), with the columns whose center falls outside the modeled area simply left out - the survivors form a single 4-connected block, so the mesh is structured with a hole-free footprint;

  • in the vertical, n_layers sigma layers per column running from the DEM bed to a lid that follows the free surface at a fixed clearance.

The lid is the part that matters for the air phase. A flat lid high above the terrain fills the domain with air that has nothing to do and every opportunity to misbehave: large recirculating air cells drive the Courant number, the atmosphere patch has to swallow the resulting in/outflow, and the time step collapses on a phase nobody is interested in. Clamping the lid to freeboard meters above the water surface (from the TELEMAC 2D result - see axqua.solvers.openfoam.hotstart) typically removes 60-90% of the air cells, keeps the atmosphere patch close to the interface where its inletOutlet/totalPressure pair behaves, and leaves the free surface room to move.

Every vertical level is computed at the plan vertices, never at cell centers, so two neighboring columns of different height still share their corner points exactly and the mesh is conformal by construction - no merging, no snapping, no hanging nodes.

axqua.solvers.openfoam.mesh.STRUCTURES_PATCH = 'structures'

Lateral faces that bound a SOLID STRUCTURE rather than the edge of the domain. Kept apart from banks because the two answer opposite questions. Water against the domain edge means the footprint was drawn too tight and the result is wrong; water against a baffle is what a fish pass IS. Merged, the wallContact monitor reported 44.8% on munich-vsf and was read as „the domain stopped the flow spreading“ - on a geometry whose 20 baffles were already 41.5% wet in the seed, at t=0, before a single 3D step.

class axqua.solvers.openfoam.mesh.PlanGrid(dx, angle, origin, nx, ny, col_id, vert_id, cell_xy, vert_xy)[Quellcode]

Bases: object

The uniform lattice of columns the 3D mesh is extruded from.

col_id / vert_id are the structured lookup tables (-1 where the cell or vertex is absent); everything downstream addresses the mesh through them.

Parameter:
dx: float
angle: float
origin: numpy.ndarray
nx: int
ny: int
col_id: numpy.ndarray
vert_id: numpy.ndarray
cell_xy: numpy.ndarray
vert_xy: numpy.ndarray
property n_columns: int
property n_vertices: int
property keep: numpy.ndarray
property cell_area: float
axqua.solvers.openfoam.mesh.flow_angle(cfg)[Quellcode]

Principal direction of the reach [rad], from the channel centerline.

Aligning the lattice with it keeps the bulk flow along a grid axis, which is where a first/second-order upwind-biased scheme is least diffusive; on a diagonal reach it also shrinks the bounding box the lattice has to cover. Falls back to 0 (grid axes = easting/northing) with no centerline.

Parameter:

cfg (Config)

Rückgabetyp:

float

axqua.solvers.openfoam.mesh.build_plan_grid(polygon, dx, *, angle=0.0, max_columns=4000000, blocked=None, blocked_by_structures=True)[Quellcode]

Lay a uniform lattice of spacing dx over polygon (a shapely geometry).

A column is kept when its center lies inside the polygon; the kept set is then hole-filled and reduced to its largest 4-connected component, so the footprint is a single block with no interior voids (an interior void would become a spurious internal wall, and a detached island a disconnected mesh region that decomposePar and the pressure solve both dislike).

blocked is geometry to cut out of the footprint. Under a rigid lid that is not only the solid structures: the dry part of the reach is cut out too, because a column with no water in it has nothing to mesh. blocked_by_structures says which, and exists only so the diagnostics name the right cause - a braided reach routinely loses a few pockets to the dry trim, and reporting that as „a solid structure has cut the domain in two“ sends the reader looking for a wall that was never drawn. Measured on inn-KB15, which has no structures at all: 2 blocks, 51 columns dropped, blamed on a structure.

It defaults to True so that a caller who cuts out geometry and says nothing still gets the warning - being told about a wall that turns out to be a dry bar is a smaller failure than not being told the domain was severed.

Parameter:
Rückgabetyp:

PlanGrid

axqua.solvers.openfoam.mesh.sigma_levels(height, n_layers, *, bed_layer=None, expansion=1.0)[Quellcode]

Per-vertex relative level positions s in [0, 1], shape (NV, n+1).

Uniform by default. expansion > 1 grows the layers from the bed upwards (ratio of top to bed layer thickness). bed_layer pins the absolute thickness of the first (bed) layer - the knob that matters here, because OpenFOAM’s rough wall functions need the first cell center to stand clear of the sand-grain roughness Ks, and on a gravel bed Ks is a sizeable fraction of the flow depth. It is capped at half the local column height so a shallow column still gets a mesh.

Because s is evaluated per vertex, neighboring columns still share their corner points exactly: the mesh stays conformal even where the layer distribution differs.

Parameter:
Rückgabetyp:

numpy.ndarray

class axqua.solvers.openfoam.mesh.OpenFoamMesh(polymesh, grid, n_layers, z, bed, lid, cell_centres, cell_column, cell_layer, column_bed, column_wse=None, column_depth=None, column_uv=None, bed_ks=None, rigid_lid=False, lid_step_median=None, lid_step_p99=None, inlet_patches=<factory>, outlet_patches=<factory>, inlet_discharge=<factory>, notes=<factory>)[Quellcode]

Bases: object

The built 3D mesh plus the bookkeeping the field writers need.

Parameter:

  • polymesh (PolyMesh)

  • grid (PlanGrid)

  • n_layers (int)

  • z (np.ndarray)

  • bed (np.ndarray)

  • lid (np.ndarray)

  • cell_centres (np.ndarray)

  • cell_column (np.ndarray)

  • cell_layer (np.ndarray)

  • column_bed (np.ndarray)

  • column_wse (np.ndarray | None)

  • column_depth (np.ndarray | None)

  • column_uv (np.ndarray | None)

  • bed_ks (np.ndarray | None)

  • rigid_lid (bool)

  • lid_step_median (float | None)

  • lid_step_p99 (float | None)

  • inlet_patches (list[str])

  • outlet_patches (list[str])

  • inlet_discharge (dict[str, float])

  • notes (list[str])

polymesh: PolyMesh
grid: PlanGrid
n_layers: int
z: np.ndarray
bed: np.ndarray
lid: np.ndarray
cell_centres: np.ndarray
cell_column: np.ndarray
cell_layer: np.ndarray
column_bed: np.ndarray
column_wse: np.ndarray | None = None
column_depth: np.ndarray | None = None
column_uv: np.ndarray | None = None
bed_ks: np.ndarray | None = None
rigid_lid: bool = False
lid_step_median: float | None = None

lid_steps() of the lid this mesh was built with - the rigid lid’s applicability test, kept rather than logged and forgotten. None under a two-phase lid, where the surface is solved and the question does not arise.

lid_step_p99: float | None = None
inlet_patches: list[str]
outlet_patches: list[str]
inlet_discharge: dict[str, float]
notes: list[str]
property n_cells: int
property wall_patches: list[str]

the bed, the domain edge, the structure sides.

Read off the mesh that was actually built rather than listed as a constant, because structures exists only when the case has solid structures to blank. The top patch is excluded even under a rigid lid, where it IS a wall: an air-water interface carries negligible shear, so it takes slip.

Type:

Every no-slip wall patch

property top_patch: str

lid under a rigid lid, atmosphere for the two-phase case.

property first_layer_height: numpy.ndarray

Thickness of the bed-adjacent layer at every plan vertex [m].

property min_layer_height: float

Thinnest layer anywhere in the mesh [m] - what limits the time step.

Not the bed layer: that one is pinned (to clear the roughness) and is usually the thickest, so the layers above it are thinner wherever a column is shallow. The Courant condition sees the thinnest of them.

property cell_bounds: tuple[numpy.ndarray, numpy.ndarray]

Lower and upper z of every cell [m a.s.l.], averaged over its 4 corners.

This is what turns a free-surface elevation into a VOF alpha field: the submerged fraction of a cell is how much of [z_lo, z_hi] lies below the surface, which puts a one-cell-sharp interface exactly where the 2D model says it is instead of the block approximation setFields would give.

axqua.solvers.openfoam.mesh.column_corner_vertices(grid)[Quellcode]

(n_columns, 4) vertex ids of each column’s plan corners, CCW in grid frame.

Parameter:

grid (PlanGrid)

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.mesh.extrude(grid, bed, lid, n_layers, *, bed_layer=None, expansion=1.0, min_height=0.05)[Quellcode]

Extrude grid between the per-vertex bed and lid into a hex mesh.

Returns (polymesh, z, cell_centres, cell_column, cell_layer, side_faces) where side_faces maps each exposed lateral face group to the arrays the patch classifier needs.

Parameter:
Rückgabetyp:

tuple[PolyMesh, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, dict]

axqua.solvers.openfoam.mesh.split_outlet_above_tailwater(owner, cell_centres, stage, margin, *, cell_size=0.0)[Quellcode]

(open faces, walled faces) for one outlet, split at the tailwater.

An outlet spans bed to lid. Below the tailwater it passes water, which is what it is for; above, it is an open pressure boundary to AIR - and the lid is open too, so the air phase gets a circulation loop through the domain with nothing setting its rate. On munich-vsf that loop grew from 7.0 to 9.7 m3/s, about 7.5 m/s through a 1.28 m2 face, and once it was large enough it reversed the outlet and dragged water in with it. Every leg - 0.045 m, 0.030 m, MULES and isoAdvector alike - lost its result the same way at about t=40, and the leak was at the downstream end where the water had backed up 28 cm above the stage it was supposed to hold.

A face is placed by its owner cell’s center, which is the only elevation a lateral face has here. The margin is what keeps this from becoming a weir: the surface at the crop face must still be free to move.

axqua.solvers.openfoam.mesh.split_structure_sides(midpoints, banks, blocked, dx)[Quellcode]

Split the lateral wall faces into (domain edge, structure sides).

The two answer opposite questions, which is the whole reason to separate them. Water against the domain edge means the footprint was drawn too tight and the result is not trustworthy; water against a baffle is a vertical-slot fish pass working as designed. Merged into one banks patch they are indistinguishable, and on munich-vsf the combined monitor read 44.8% and was reported as „the footprint stopped the flow spreading“ - on a geometry whose twenty baffles were already 41.5% wet in the TELEMAC seed, at t=0, before a single 3D step was taken.

A face belongs to a structure when its midpoint falls on the blanked footprint grown by half a cell diagonal. The slack is needed, not cosmetic: the face lies on the boundary of the blanked columns rather than inside them, and blanking is by cell-square intersection, so a face can legitimately stand that far from the solid it bounds. Growing by more would start claiming the domain edge where a structure happens to sit against it.

Parameter:
Rückgabetyp:

tuple[numpy.ndarray, numpy.ndarray]

axqua.solvers.openfoam.mesh.classify_sides(cfg, midpoints, tolerance)[Quellcode]

Tag each exposed lateral face as an inlet, an outlet, or a bank wall.

Reuses the case’s boundaries.liquid_boundaries layer - the same lines that drive the TELEMAC .cli - so the 2D and 3D models take their flow in and out at exactly the same places. Each line becomes its own patch (inlet-1, inlet-2, …) so a reach with two feeding branches can prescribe each branch’s own discharge, and the per-line Target flow field (when present) supplies it.

Parameter:
Rückgabetyp:

tuple[numpy.ndarray, list[str], dict[str, float]]

axqua.solvers.openfoam.mesh.resolve_headroom(cfg, state, *, notes=None)[Quellcode]

The freeboard this build will use [m], derived or as configured.

Under headroom_mode: auto the configured freeboard is a floor and the flow decides the rest (State2D.headroom()) - so a reach fast enough to convert real velocity head gets the air it needs to show it, and a slow one does not pay for air it will never use. Under a rigid lid there is no freeboard at all.

Parameter:
Rückgabetyp:

float

axqua.solvers.openfoam.mesh.resolve_margin(cfg, state)[Quellcode]

The plan-footprint buffer past the 2D wetted edge [m], derived or configured.

A surface free to rise is also free to spread, and a footprint trimmed to the 2D water line would stop it with a wall (see State2D.lateral_margin()).

Parameter:

cfg (Config)

Rückgabetyp:

float

axqua.solvers.openfoam.mesh.resolve_lid(cfg, grid, bed, wse_v, freeboard, *, notes=None)[Quellcode]

The lid elevation at every plan vertex [m a.s.l.].

Four cases, and which one applies decides what the model can say about the water surface:

  • rigid lid - the lid is wse_v, the converged 2D free surface, sampled per vertex. Sloping, stepped, and prescribed: no freeboard, because a freeboard would be air;

  • ``lid: follow`` with a 2D seed - a freeboard above that same surface, then smoothed, so most air cells never exist but the surface can still move;

  • ``lid: flat``, or no seed at all - a horizontal plane over the whole reach, which is correct and carries the full air column;

  • ``lid_elevation`` - pinned where the case says, in any mode.

The floor keeps the column from closing onto the bed: min_water_depth under a rigid lid, min_column_height otherwise.

Parameter:
Rückgabetyp:

np.ndarray

axqua.solvers.openfoam.mesh.resolve_layers(cfg, bed_ks, column_height=None, *, notes=None)[Quellcode]

Sigma layers per column, after the rigid-lid roughness check.

Under a rigid lid the layers span the WATER only (isar: 0.2-1.0 m), not water plus freeboard (0.7-1.8 m), so a count carried over from a two-phase case gives cells three to four times thinner. Past a point they are thinner than the grains on the bed - meaningless, and thin enough that the bed’s variation within one plan cell folds them, which checkMesh rejects as incorrectly oriented face pyramids.

The check is sized on the SHALLOWEST meshed column (min_water_depth), so a reach whose typical column is much deeper than that floor pays for its own margins. openfoam.auto_layers: false says the case has looked at that trade and wants its layer count anyway; the note then records what it costs in the shallow columns, which is where the folded cells would appear.

Parameter:

  • cfg (Config)

  • bed_ks (np.ndarray | None)

  • column_height (np.ndarray | None)

  • notes (list[str] | None)

Rückgabetyp:

int

axqua.solvers.openfoam.mesh.plan_spacing(cfg)[Quellcode]

The plan lattice spacing [m], from cell_size or cell_size_factor.

The factor form is relative to the TELEMAC channel resolution, so „three times coarser than the 2D mesh“ is a single number that keeps meaning across cases instead of a meter value that has to be recomputed per reach.

Parameter:

cfg (Config)

Rückgabetyp:

float

axqua.solvers.openfoam.mesh.lid_steps(grid, lid, bed, *, radius=2)[Quellcode]

(median, p99) of the lid’s local range as a fraction of the water depth.

The applicability test for a rigid lid, and the one that was missing when munich-vsf spent 61 hours producing a field its velocity cap was holding together.

A rigid lid is a slip wall at the free surface, so wherever the prescribed surface drops sharply the ceiling becomes a near-vertical step and the flow beneath it is squeezed through a sluice that does not exist. A real free surface answers a drop by plunging; a lid cannot, so it converts the head into velocity instead - sqrt(2 g dz), which for a 1 m step is 4.4 m/s in a reach whose own water moves at 0.4.

Small everywhere (a few per cent) means the surface is smooth on the scale of a cell and the lid is a fair approximation. Order 1 means the surface has steps as large as the water is deep - a weir, a drop structure, the slots of a fish pass - and those are exactly the places a rigid lid must not be used.

Parameter:
Rückgabetyp:

tuple[float, float]

axqua.solvers.openfoam.mesh.build_mesh(cfg, *, state=None, dem=None, tailwater=None)[Quellcode]

Build the terrain-following hex mesh for cfg’s OpenFOAM case.

state is an optional State2D (the converged TELEMAC 2D result). With one, the lid follows the free surface and the footprint can be trimmed to the wetted corridor; without one, the lid is flat and the footprint is the whole ROI - correct, but far more air. tailwater is the level of the outflow the build resolved, above which an outlet is walled.

Parameter:
Rückgabetyp:

OpenFoamMesh

Mesh quality and physical validity for the OpenFOAM case, before OpenFOAM runs.

Two separate questions, both answered here because both must be answered before committing solver hours:

Is the mesh numerically sound? Non-orthogonality, skewness, aspect ratio, cell volume - the same measures checkMesh reports, computed in numpy so the answer arrives while the mesh is still in memory and a bad cell_size / n_layers combination can be corrected in seconds rather than after a decomposition and a failed run. checkMesh is still run by axqua.solvers.openfoam.runtime and remains the authority; this is the fast pre-check, exactly as axqua.mesh_quality is for the TELEMAC mesh.

Is the wall treatment admissible? This is the check that matters on a gravel-bed river and the one that has no analogue in an aerodynamics workflow. OpenFOAM’s nutkRoughWallFunction manipulates the log-law E parameter for a sand-grain roughness Ks, and that manipulation is only meaningful while the first cell center stands clear of the roughness elements - y_1 > Ks, i.e. a bed layer at least 2 Ks thick. On the isar reach ks runs 0.085-0.22 m in 0.3-0.9 m of water, so a „well resolved“ bed layer of a few centimeters would put the first cell center inside the gravel and the wall function would return a confident, meaningless answer. Refining towards the bed is therefore not available here, and the report says so per bed face instead of letting the run imply otherwise.

class axqua.solvers.openfoam.quality.MeshReport(n_cells, n_faces, n_points, n_internal_faces, patches=<factory>, min_volume=0.0, max_volume=0.0, total_volume=0.0, max_non_ortho=0.0, mean_non_ortho=0.0, n_non_ortho_severe=0, max_skewness=0.0, n_skew=0, n_bad_pyramid=0, max_aspect_ratio=0.0, n_bad_volume=0, bed_layer_min=0.0, bed_layer_median=0.0, ks_min=None, ks_max=None, n_bed_faces=0, n_wall_function_invalid=0, ks_plus=None, regime='unknown', warnings=<factory>)[Quellcode]

Bases: object

Everything measurable about the mesh without launching OpenFOAM.

Parameter:
n_cells: int
n_faces: int
n_points: int
n_internal_faces: int
patches: dict[str, int]
min_volume: float = 0.0
max_volume: float = 0.0
total_volume: float = 0.0
max_non_ortho: float = 0.0
mean_non_ortho: float = 0.0
n_non_ortho_severe: int = 0
max_skewness: float = 0.0
n_skew: int = 0
n_bad_pyramid: int = 0
max_aspect_ratio: float = 0.0
n_bad_volume: int = 0
bed_layer_min: float = 0.0
bed_layer_median: float = 0.0
ks_min: float | None = None
ks_max: float | None = None
n_bed_faces: int = 0
n_wall_function_invalid: int = 0
ks_plus: float | None = None
regime: str = 'unknown'
warnings: list[str]
property usable: bool

True when the mesh is structurally sound enough to run on.

Deliberately not the same as passing every checkMesh check. A mesh derived from a real DEM will occasionally carry a handful of marginal faces where the survey has a seam or a step - the isar reach produces exactly one skew face over the limit in 4.9 million, and no amount of lid smoothing moves it because it comes from the bathymetry, not from the meshing. Refusing to run on that would be refusing to run on the terrain. What genuinely stops a run is a cell with no volume or an unclosed cell, and those are what this reports.

property ok: bool

True when the mesh would pass every checkMesh check.

property marginal_faces: int

Faces that exceed a checkMesh limit without breaking the mesh.

lines()[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.quality.face_geometry(mesh)[Quellcode]

(face centers, face area vectors) by OpenFOAM’s triangle decomposition.

A quad on a terrain-following mesh is generally not planar - its four corners sit on four columns of different height - so its centroid is not the average of its corners. OpenFOAM fans the face about that average and area-weights the resulting triangles; anything cruder shifts the center by a few percent of the cell, which is enough to move the reported non-orthogonality by several degrees.

Rückgabetyp:

tuple[numpy.ndarray, numpy.ndarray]

axqua.solvers.openfoam.quality.cell_geometry(mesh, face_centres, face_normals)[Quellcode]

(cell centers, cell volumes) by OpenFOAM’s pyramid decomposition.

Each face forms a pyramid with an estimated cell center; the true center is the volume-weighted mean of the pyramid centroids. Reproducing this exactly is what makes the non-orthogonality reported here the same number checkMesh will report, rather than an optimistic approximation of it.

Parameter:
Rückgabetyp:

tuple[numpy.ndarray, numpy.ndarray]

axqua.solvers.openfoam.quality.bad_face_pyramids(mesh, face_centres, face_normals, cell_centres)[Quellcode]

Faces whose pyramid to their cell center has non-positive volume.

checkMesh calls these „incorrectly oriented“ and treats them as fatal, and rightly: the cell is folded even though its total volume can still come out positive - which is exactly how a first version of this report passed a mesh checkMesh rejected. A rigid-lid mesh produces them at the waterline, where the water column thins to nothing and the cells collapse.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.quality.non_orthogonality(mesh, face_normals, cell_centres)[Quellcode]

(per-face angle [deg], mean [deg]) for the internal faces.

The mean is checkMesh’s: it averages the cosines and takes the arccosine of that, not the mean of the angles - the two differ by a couple of degrees on a mesh like this, and quoting the second while calling it the first would make two reports of the same mesh disagree for no reason.

Parameter:
Rückgabetyp:

tuple[numpy.ndarray, float]

axqua.solvers.openfoam.quality.skewness(mesh, face_centres, face_normals, cell_centres)[Quellcode]

Face skewness, matching OpenFOAM’s primitiveMeshTools::faceSkewness.

The skewness vector is the offset of the owner-neighbour line’s crossing point from the face center; what matters is the normalization. OpenFOAM divides it by the face’s own half-width in the skew direction (floored at 0.2 |d| for an internal face, 0.4 |d| for a boundary one), not by the cell-to-cell distance. Dividing by |d| instead - the intuitive choice - understates the result by roughly the face aspect ratio, which on this mesh is a factor of ~7: the number would look reassuring while sitting near checkMesh’s limit.

Boundary faces are included, as checkMesh includes them, by mirroring the owner cell across the face.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.quality.aspect_ratios(mesh, face_normals, volumes)[Quellcode]

Cell aspect ratio, matching OpenFOAM’s primitiveMeshTools::cellClosedness.

max(component ratio, hydraulic ratio) where the component ratio is the largest over the smallest of the per-direction summed face-area magnitudes, and the hydraulic one is (1/6) sum|Sf| / V^(2/3). The component form dominates on a terrain-following mesh - a cell a meter wide in plan and a decimeter tall has ten times the horizontal face area of the vertical - and it is precisely the number that exposes how flat the layers are, so taking only the hydraulic form (which reads ~3 where OpenFOAM reads ~17) would hide it.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.quality.assess(of_mesh, *, state=None, roughness_constant=0.5)[Quellcode]

Measure of_mesh (an OpenFoamMesh).

state is the optional 2D hotstart; with one, the roughness regime (ks+) is evaluated at the reach’s own depth and velocity rather than guessed.

Parameter:

roughness_constant (float)

Rückgabetyp:

MeshReport

Ausgangszustand aus einem TELEMAC-Ergebnis

Seed an OpenFOAM case from the converged TELEMAC 2D result.

Starting a VOF run from a dry bed, or from a flat block of water, means paying for the whole filling transient in the most expensive model in the chain. The 2D result already knows where the water is, how deep it is and how fast it moves, so the 3D run can start from that state and spend its time on what only 3D can give - the vertical structure, the secondary currents, the local surface deformation.

The same state is used for three separate jobs, which is why it lives in one class:

  • the footprint - the 3D mesh only needs to cover the wetted corridor plus a margin, not the whole ROI (State2D.wet_footprint());

  • the lid - clamped a fixed freeboard above the 2D free surface, which is what keeps the air phase small (State2D.sample_surface());

  • the initial fields - alpha.water from the free-surface elevation and U from the depth-averaged velocity (State2D.sample_columns()).

Sampling is nearest-node on a KD-tree rather than P1 interpolation on the element table - see State2D.tree for why a merged parallel SELAFIN leaves no usable triangulation, and why that costs nothing when the 2D mesh is finer than the 3D lattice it seeds.

class axqua.solvers.openfoam.hotstart.State2D(x, y, triangles, depth, surface, bottom, u, v, time=0.0, source=None, u3d=None, v3d=None, z3d=None, _tree=None, _spacing=None)[Quellcode]

Bases: object

A depth-averaged flow state on the TELEMAC mesh.

Parameter:

  • x (np.ndarray)

  • y (np.ndarray)

  • triangles (np.ndarray)

  • depth (np.ndarray)

  • surface (np.ndarray)

  • bottom (np.ndarray)

  • u (np.ndarray)

  • v (np.ndarray)

  • time (float)

  • source (Path | None)

  • u3d (np.ndarray | None)

  • v3d (np.ndarray | None)

  • z3d (np.ndarray | None)

  • _tree (object)

  • _spacing (float | None)

x: np.ndarray
y: np.ndarray
triangles: np.ndarray
depth: np.ndarray
surface: np.ndarray
bottom: np.ndarray
u: np.ndarray
v: np.ndarray
time: float = 0.0
source: Path | None = None
u3d: np.ndarray | None = None
v3d: np.ndarray | None = None
z3d: np.ndarray | None = None
classmethod from_slf(path, *, frame=-1)[Quellcode]

Read path (a TELEMAC r2d.slf), by default its last frame.

Parameter:
Rückgabetyp:

State2D

property has_profile: bool

Whether this seed knows how the velocity varies over the depth.

sample_profile(xy, z)[Quellcode]

(n, 2) horizontal velocity at the 3D points (xy, z).

Nearest plan node on the same KD-tree the 2D path uses, then linear in z between that node’s own levels. Outside the column the end level is held rather than extrapolated: below the bed and above the surface there is no profile to continue, and a linear extrapolation there would invent a jet.

Falls back to the depth-averaged velocity when the seed carries no profile, so a caller never has to ask which kind of seed it has.

Parameter:
Rückgabetyp:

numpy.ndarray

property tree

KD-tree over the 2D nodes, built once.

Sampling is nearest-node, deliberately, not P1 interpolation on the element table. A SELAFIN merged from a parallel run is not a clean triangulation - the isar result carries 4,393 coincident nodes and 66,759 zero-area triangles out of 460,314, artifacts of the subdomain interfaces that TELEMAC ignores but that make matplotlib’s TrapezoidMapTriFinder reject the mesh outright. Repairing that into something a trifinder accepts is guesswork about which overlapping element is authoritative.

Nearest-node needs none of it and costs nothing in accuracy here, because the 2D mesh is finer than the 3D lattice it is being sampled onto (isar: 0.3 m channel cells into a 0.5 m lattice): this is downsampling, and the residual error is the surface slope over one 2D cell - on the order of a millimeter. resolution_check() warns if that relationship ever reverses.

property node_spacing: float

Median distance to a node’s nearest neighbour [m].

resolution_check(cell_size)[Quellcode]

Warn when the 3D lattice is finer than the 2D mesh that seeds it.

Parameter:

cell_size (float)

Rückgabetyp:

str | None

wet(wet_depth=0.01)[Quellcode]
Parameter:

wet_depth (float)

Rückgabetyp:

numpy.ndarray

velocity_scale(wet_depth=0.01)[Quellcode]

Representative water speed [m/s]: the 95th percentile over wet nodes.

Used to size the limitVelocity cap and the turbulence seed. A percentile rather than the maximum, because a single partially-wet cell in a 2D result can carry an absurd velocity that would set a useless cap.

Parameter:

wet_depth (float)

Rückgabetyp:

float

depth_scale(wet_depth=0.01)[Quellcode]

Representative flow depth [m]: the median over wet nodes.

Parameter:

wet_depth (float)

Rückgabetyp:

float

headroom(floor=0.0, *, wet_depth=0.01, depth_allowance=0.25)[Quellcode]

How far above the 2D surface the lid must stand [m].

The seed is an approximation, and the 3D run has to be free to disagree with it - if the surface rises into the lid, the result is quietly constrained by a mesh decision rather than by the flow. The room it needs has two parts:

  • the velocity head V^2/2g the flow could convert at a stagnation point, an obstruction or the outside of a bend. Taken at the 95th percentile speed, since one partially-wet 2D cell can carry an absurd velocity;

  • a fraction of the depth (depth_allowance) for everything 2D cannot see at all - secondary currents, local surface deformation, the standing waves a depth-averaged model has no equation for.

floor (the configured freeboard) is a minimum, never a cap, so a case that already runs keeps at least the air it had. Air is the expensive part of a VOF run, which is why this is derived rather than set generously: on a slow reach it stays at the floor, and on a fast one it grows where it must.

Parameter:
Rückgabetyp:

float

lateral_margin(floor=0.0, *, bank_slope=0.1, wet_depth=0.01)[Quellcode]

How far past the 2D wetted edge the plan footprint must reach [m].

The same argument one dimension over: if the surface may rise by headroom(), the water line may move sideways, and on a bank of slope bank_slope a rise h moves it h / slope. 10% is a plain gravel bank; a steeper one needs less room, and the configured wet_margin floor covers the case where the bank is a wall.

Parameter:
Rückgabetyp:

float

sample_surface(xy)[Quellcode]

Free-surface elevation [m a.s.l.] at arbitrary plan points.

Where the 2D model is dry the free surface is the bed, so the lid derived from this naturally hugs the ground over dry bars instead of flying over them - which is the whole point of following the surface.

Parameter:

xy (numpy.ndarray)

Rückgabetyp:

numpy.ndarray

sample_columns(xy, *, wet_depth=0.01)[Quellcode]

(free surface, depth, [u v]) at the plan points of the mesh columns.

Velocity is zeroed wherever the 2D depth is below wet_depth: a depth-averaged velocity carried by a millimeter of water is numerically large and physically meaningless, and seeding it into a 3D run puts a spurious jet on the bank at t = 0.

Parameter:
Rückgabetyp:

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

wet_footprint(*, wet_depth=0.01, buffer=5.0, resolution=None)[Quellcode]

Plan polygon of the wetted area, dilated by buffer meters.

Built by rasterising the wet nodes and polygonising the mask rather than by unioning wet triangles: this mesh has 460k of them, and a shapely union over that many polygons is minutes of work for a footprint that only needs to be accurate to a cell. Small holes are closed first, so an isolated dry node mid-channel does not punch a hole through the 3D domain.

Returns the largest resulting polygon, or None when nothing is wet.

Parameter:
dry_footprint(*, wet_depth=0.01, margin=0.0, resolution=None)[Quellcode]

Plan polygon of everything that is not wet, for a water-only domain.

The complement of wet_footprint(), so a rigid-lid mesh can blank dry bars the same way it blanks a building. margin grows the water body before taking the complement, keeping a little dry ground inside the domain - without it the domain edge lands exactly on the waterline and the inflow/outflow lines, which sit at that edge, fall outside the mesh and no inlet patch is found.

Parameter:
summary(*, wet_depth=0.01)[Quellcode]
Parameter:

wet_depth (float)

Rückgabetyp:

list[str]

axqua.solvers.openfoam.hotstart.load_hotstart(cfg, path=None)[Quellcode]

Best-effort 2D hotstart for cfg: the case’s own r2d.slf unless given.

Returns None (with a warning) rather than raising when no 2D result exists - a case can still be meshed and run without one, it just carries a flat lid and a cold start, and the caller reports that.

Parameter:

path (str | Path | None)

Rückgabetyp:

State2D | None

The 0/ directory: initial and boundary conditions for the interFoam case.

Randbedingungen

Every type here was checked against the OpenFOAM 9 source and the interFoam/RAS/{waterChannel,weirOverflow} tutorials, because the difference between a river case that runs and one that oscillates is almost entirely in this table.

patch

alpha.water

U

p_rgh

inlet

variableHeightFlowRate

flowRateInletVelocity (mass flow)

fixedFluxPressure

outlet

inletOutlet at the stage

pressureInletOutletVelocity

prghPressure

bed, banks

zeroGradient

noSlip

fixedFluxPressure

atmosphere

inletOutlet (0)

pressureInletOutletVelocity

totalPressure 0

Three of those deserve their reasoning recorded.

The inlet does not prescribe a water level. variableHeightFlowRate lets alpha at the inlet float between 0 and 1, so the water on the patch stands as high as the water behind it. The inlet stage is an outcome, not an input - which is what you want, because forcing both Q and the depth at an inlet over-determines it and shows up as a standing wave a few cells downstream.

The inlet prescribes the mass that enters, not the volume. One velocity for the whole patch, the mass flow of the water divided by the mass per length that stands on the patch (flowRateInletVelocity with massFlowRate and the density field). Whatever part of the patch is wet then carries the discharge, and the air above it enters at the speed of the water, which nothing notices. The condition written for this purpose, variableHeightFlowRateInletVelocity, sets the velocity in proportion to alpha and scales it so that water and air together make up the discharge: the water that enters is Q <alpha^2> / <alpha>, short by whatever the surface is smeared over on the patch. Measured on the Isar example with the same state: 1.434 of 1.6 m3/s and 0.743 of 0.8 m3/s with that condition, 1.596 and 0.798 with the mass flow (on munich-vsf the older runs delivered 82 to 84 %). What is left, 0.2 to 0.3 %, is the share of the air in the mass on the patch.

The outlet holds a stage. The tutorials use a free outfall (zeroGradient everywhere), which lets the model choose its own tailwater. This case has a calibrated one - the same rating curve the TELEMAC run used - so the outlet gets prghPressure with a hydrostatic p profile below the target level and an inletOutlet alpha profile that is 1 below it and 0 above. That is the direct 3D analogue of TELEMAC’s PRESCRIBED ELEVATIONS, so 2D and 3D are driven identically and their results are comparable. boundaries.outflow_condition: free still selects the free outfall.

The air on the outlet is at rest, like the air under the open top. The solver carries p_rgh = p + rho g (z - hRef), and totalPressure 0 on the open top says p_rgh = 0 there, which is air at rest at any elevation. The outlet prescribes p, so above the stage it has to prescribe the pressure of that same air, -rho_air g (z - hRef), and the water below the stage hangs on the air pressure at the stage. A profile that is simply zero above the stage disagrees with the open top by rho_air g (z - hRef): 8 kPa on a reach at 816 m a.s.l. with hRef left at zero, which is a wind of the speed the velocity limit allows between the outlet and the top, and it is that wind, in a handful of cells, that set the time step of the whole run (isar: 0.0013 s; see reference_level()). constant/hRef is written with the lowest outlet stage, so that p_rgh is the piezometric head above the tailwater and stays small.

Ausgangsbedingungen

alpha.water and U come from the 2D hotstart (see axqua.solvers.openfoam.hotstart) and are the only non-uniform internal fields; p_rgh starts at 0 (as in every tutorial - the first pressure solve resolves it) and the turbulence fields start uniform at scales derived from the reach’s own depth and velocity.

axqua.solvers.openfoam.fields.scalar_list(values)[Quellcode]

nonuniform List<scalar> body.

Parameter:

values (numpy.ndarray)

Rückgabetyp:

str

axqua.solvers.openfoam.fields.vector_list(values)[Quellcode]

nonuniform List<vector> body.

Parameter:

values (numpy.ndarray)

Rückgabetyp:

str

axqua.solvers.openfoam.fields.render_field(cls, obj, dimensions, internal, boundary, *, time_dir='0')[Quellcode]

Assemble one OpenFOAM field file.

Parameter:
Rückgabetyp:

str

axqua.solvers.openfoam.fields.initial_alpha(of_mesh, *, wet_depth=0.01)[Quellcode]

Per-cell water fraction from the hotstart free surface.

The submerged fraction of each cell, so the interface is sharp to within one cell wherever the 2D surface crosses a cell. With no hotstart the domain starts dry, which is legitimate but means paying for the whole filling transient in the most expensive model of the chain.

A column that is dry in the 2D result starts dry. The free surface of a dry 2D node is its own bed, and the bed of the column comes from the terrain model, which resolves what the 2D mesh does not. Wherever the terrain at the column lies below the bed of the nearest dry node, „below the 2D free surface“ is a layer of water on dry ground. On the Isar example at 1 m cells that was 583 m3, a sixth of the water of the reach, on 55 % of the dry columns (median 0.03 m, 0.49 m at the 99th percentile), and it ran down the banks from the first time step at the speed the velocity limit allows.

Parameter:

wet_depth (float)

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.fields.initial_velocity(of_mesh, alpha, state=None)[Quellcode]

Per-cell velocity: the seed’s velocity in water, zero in air.

From a 2D seed this is uniform over the depth, deliberately, rather than an assumed log profile: the log profile’s normalization is only meaningful where the wall function is admissible, and on a gravel bed that is often nowhere (see axqua.solvers.openfoam.quality). The vertical structure is what the 3D run is for, so it is left to develop during the spin-up rather than invented here.

From a 3D seed (pre_run.dimension: 3d) the profile is not assumed - it was computed, by a solver, on this reach’s own bathymetry - so each cell is sampled at its own elevation instead. That is the one case where prescribing vertical structure is better than starting flat.

Scaling by alpha keeps an interface cell from carrying a full-strength velocity into the air.

Parameter:

alpha (numpy.ndarray)

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.fields.turbulence_scales(state, cfg)[Quellcode]

(k, epsilon, omega) seeds from the reach’s own depth and velocity.

Rückgabetyp:

tuple[float, float, float]

axqua.solvers.openfoam.fields.outlet_pressure(z, stage, *, water_density, air_density, reference)[Quellcode]

Static pressure [Pa] on outlet faces at the elevations z for a tailwater at stage.

Air at rest above the stage, and still water below it that hangs on the pressure of that air at the stage. reference is the level hRef the solver counts the elevation from. prghPressure turns the profile into p_rgh = p + rho g (z - hRef) with the density that is at the face, which is zero wherever there is air: the value totalPressure 0 gives the open top. With a profile of zero above the stage the two boundaries disagree by rho_air g (z - hRef) and air is driven from one to the other.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.fields.REFERENCE_FILE = 'hRef'

the level the two-phase solver counts the elevation from.

Type:

constant/hRef

axqua.solvers.openfoam.fields.reference_level(of_mesh, stages)[Quellcode]

The level hRef [m a.s.l.] of the two-phase case.

The solver works with p_rgh = p + rho g (z - hRef). Where the water surface stands at eta, p_rgh jumps across it by (rho_water - rho_air) g (eta - hRef), and the pressure solution is only as exact as the tolerance of its solver times that jump. Left at zero on a mesh in meters above sea level the jump is 8 MPa where the flow is driven by a few pascals. So the level is put where the water is: at the lowest stage the outlets hold, which also makes p_rgh / (rho g) the piezometric head above the tailwater. A free outfall prescribes no pressure, and any level of the reach serves: the lowest point of the lid.

Parameter:

stages (Mapping[str, float] | None)

Rückgabetyp:

float

axqua.solvers.openfoam.fields.write_reference_level(case_dir, level)[Quellcode]

Write constant/hRef, or remove it where level is None.

A case that is rebuilt single-phase in the folder of a two-phase one must not keep the file: the solver reads it whenever it is there, and the pressure the rigid lid prescribes at its outlet counts from zero.

Parameter:
Rückgabetyp:

Path | None

axqua.solvers.openfoam.fields.outlet_stages(outflow_stage, patches)[Quellcode]

Per-patch tailwater level [m a.s.l.], broadcasting a scalar.

A single number is the right answer for a whole-reach case, which has one outfall. A sub-model crop does not: cutting a reach mid-stream leaves several open faces at different stations, and the free surface across them spans whatever the parent’s water-surface slope puts there - 1.18 m on the KB15 40 m crop. Prescribing one level on all of them would impose a tailwater the parent model never had on at least one face.

A mapping must name every outlet patch. Silently defaulting the ones it omits would reintroduce exactly the bug this exists to prevent, in the case where it is hardest to notice.

Rückgabetyp:

dict[str, float] | None

axqua.solvers.openfoam.fields.rating_entry(table)[Quellcode]

A rating curve as the table of a Function1: table ((Q h) ...).

Rückgabetyp:

str

axqua.solvers.openfoam.fields.write_fields(of_mesh, cfg, case_dir, *, state=None, outflow_stage=None, discharges=None, uniform_bed_ks=False, rating=None)[Quellcode]

Write 0/{alpha.water, U, p_rgh, nut, k, omega|epsilon}.

outflow_stage is the water level [m a.s.l.] the outlet patches hold; None selects the free-outfall boundary instead (outflow_condition: free). discharges maps each inlet patch to its own Q [m3/s]. uniform_bed_ks writes one global bed roughness instead of the per-face list (see _wall_entries()); set it for a calibration template. rating is a stage-discharge relation [(Q, h), ...]: with it, every outlet of a two-phase case holds the level that the relation gives at the discharge that leaves through it, with the boundary condition waterLevel_alpha_prgh of the BAW (openfoam.outflow_boundary: rating-curve). outflow_stage is then the level the run starts from.

Parameter:
Rückgabetyp:

list[Path]

Wörterbücher und Fall Assembly

constant/ and system/ dictionaries for the interFoam case.

This module is where the air phase is dealt with. interFoam in OpenFOAM 9 is already a PIMPLE solver, so the question is never „which solver“ but „which settings“, and the settings that matter are these four:

  1. ``MULESCorr yes`` (semi-implicit MULES). With the default explicit MULES the phase-fraction equation carries its own Courant limit and the tutorials run at maxAlphaCo 0.2. The semi-implicit limiter removes that constraint, so the time step is set by the flow rather than by the interface and maxCo / maxAlphaCo can go to ~0.9 - a three- to four-fold larger step for the same physics.

  2. ``limitVelocity`` in ``system/fvOptions``. Air is a thousand times lighter than water, so a small pressure error accelerates it enormously; the resulting air jet then sets the Courant number and the time step collapses on a phase nobody is interested in. Capping |U| at a few times the expected water velocity stops that dead. Water never reaches the cap, so the water solution is untouched.

  3. ``momentumPredictor no``. Standard for gravity-driven free-surface flow: the momentum predictor buys nothing when buoyancy dominates and costs an extra solve.

  4. ``nNonOrthogonalCorrectors 2``. A terrain-following mesh is non-orthogonal wherever the bed slopes; without the correctors the pressure solution leaks that as spurious velocity near the bed. Together with nCorrectors 3 and the outer correctors this is 18 pressure solutions per time step, which is most of what a step costs. Both counts are settings of the case (openfoam.pressure_correctors, non_orthogonal_correctors): on the Isar example, whose worst face is 57 degrees off, 2 and 1 gave the same depths, velocities and discharges to three digits after 40 s at 0.6 of the cost.

The lid clamping that removes most air cells altogether is in axqua.solvers.openfoam.mesh, and matters more than all four.

Staging

Two dict sets are written, under system/stage1-spinup/ and system/stage2-run/; axqua.solvers.openfoam.runtime copies one into system/ to select it. Stage 1 settles the interface from the 2D hotstart at a tight Courant number with diffusive schemes; stage 2 restarts from it and runs the production settings. The split exists because the first seconds after a hotstart are the least like the converged flow - the 2D state is depth-averaged, so there is no vertical structure at all at t = 0 - and paying for robustness there is far cheaper than paying for it over the whole run.

class axqua.solvers.openfoam.dicts.Stage(name, start_from, end_time, max_courant, alpha_scheme, c_alpha, velocity_scheme, n_outer_correctors, purpose)[Quellcode]

Bases: object

One run stage’s numerics.

Parameter:
name: str
start_from: str
end_time: float
max_courant: float
alpha_scheme: str
c_alpha: float
velocity_scheme: str
n_outer_correctors: int
purpose: str
axqua.solvers.openfoam.dicts.stages(cfg, *, single=False)[Quellcode]

The run’s stages, in order.

single collapses a two-phase case to ONE full-length production stage. It exists for the Bayesian campaign: HydroBayesCal copies a case template and makes a single solver call, so a two-stage case would run only its spin-up - 30 s at half interface compression and upwinded momentum - while the config claimed the full end_time. The interface is settled once up front instead (calibration._prespin_template) and promoted into the template’s 0/.

Parameter:

single (bool)

Rückgabetyp:

list[Stage]

axqua.solvers.openfoam.dicts.transport_properties(cfg)[Quellcode]
Rückgabetyp:

str

axqua.solvers.openfoam.dicts.turbulence_properties(cfg)[Quellcode]

constant/turbulenceProperties - the ESI name for the RAS settings.

Foundation 8+ called this momentumTransport and keyed the model as model; ESI v2406 knows neither, so a Foundation-shaped file is read as an empty dictionary and the run silently falls back to the built-in defaults.

Under kEpsilon the closure coefficients are written explicitly, at the model’s own defaults. Two reasons, both about calibration: OpenFOAM falls back to built-in values for any coefficient it does not find in the file, so a calibration that perturbs one which is absent would leave every run using the same value and never say so; and HydroBayesCal’s update_dictionary_entry requires the key to exist, raising Key 'Cmu' not found otherwise. Other models (kOmegaSST, the default) are written exactly as before.

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.gravity()[Quellcode]

constant/g. z is up: the mesh is in a projected CRS with real elevations.

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.BAW_LIBRARY = 'lib_BAW_public_BCs_v2412_20260813.so'

The library of the boundary conditions of the BAW that the OpenFOAM installation of aXqua builds (BAW_LIBRARY of its installer), with waterLevel_alpha_prgh.

axqua.solvers.openfoam.dicts.uses_rating_outlet(cfg)[Quellcode]

Whether the outlets follow a stage-discharge relation (needs a free surface).

Rückgabetyp:

bool

axqua.solvers.openfoam.dicts.control_dict(cfg, stage, *, patches, boundary_patches=None)[Quellcode]
Parameter:
Rückgabetyp:

str

axqua.solvers.openfoam.dicts.stage_write_interval(of, stage)[Quellcode]

Simulated time between two saved results of stage [s].

write_interval, adjusted so that the end of the stage is a write time. The solver saves at whole intervals counted from the start of its run and at no other time, also not at its end. A first stage of 2 s with results every 50 s therefore ended without a result, the second stage started from the latest one it found, which was the start, and the first stage had been computed for nothing: seen on the Isar example, where nothing but the time in the listing said so. With the defaults (30 s and 10 s) the end happened to be a write time. The same holds for the end of the run. The second stage starts at spinup_time, so its results come at spinup_time plus whole intervals: with 2 s, 3000 s and 100 s the last one would be at 2902 s and the rest of the run would not be saved. The interval becomes 99.93 s.

Parameter:

stage (Stage)

Rückgabetyp:

float

axqua.solvers.openfoam.dicts.monitor_interval(of)[Quellcode]

Simulated-time spacing of the discharge monitors [s].

A tenth of the field write interval, so a run writing ten result frames leaves a hundred discharge samples - enough for the ten-sample steady window axqua.solvers.openfoam.report looks for, without a line per time step.

Rückgabetyp:

float

axqua.solvers.openfoam.dicts.fv_schemes(cfg, stage)[Quellcode]
Parameter:

stage (Stage)

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.ISO_RECONSTRUCTION = ('isoAlpha', 'gradAlpha', 'plicRDF')

isoAdvector reconstruction schemes ESI v2406 ships, in increasing order of fidelity. isoAlpha is the original isosurface cut; the two PLIC schemes fit a plane per interface cell, plicRDF iterating it against a reconstructed distance function. On aXqua’s own structured hex lattice - isoAdvector’s best case - PLIC is the one worth having.

axqua.solvers.openfoam.dicts.fv_solution(cfg, stage)[Quellcode]
Parameter:

stage (Stage)

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.fv_options(cfg, velocity_cap)[Quellcode]
Parameter:

velocity_cap (float)

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.decompose_par_dict(cfg)[Quellcode]

system/decomposeParDict for the processor count of the case.

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.decomposition_dict(n_processors)[Quellcode]

A decomposeParDict that splits a case into n_processors domains.

Parameter:

n_processors (int)

Rückgabetyp:

str

axqua.solvers.openfoam.dicts.decomposition_of(case_dir)[Quellcode]

numberOfSubdomains of the case’s system/decomposeParDict, None without one.

Parameter:

case_dir (str | Path)

Rückgabetyp:

int | None

axqua.solvers.openfoam.dicts.decomposition_for(case_dir, n_processors)[Quellcode]

The decomposition dictionary of a run on n_processors ranks, for -decomposeParDict.

None where the case’s own system/decomposeParDict names that number: the run uses it, with its method and coefficients, which may be the user’s. Otherwise aXqua’s scotch decomposition for the number is written to system/decomposeParDict-<n> and that file is returned; decomposePar and the solver are both given it, so the case’s own file is never changed (a run on other ranks, killed or not, leaves nothing behind in it).

Parameter:
Rückgabetyp:

Path | None

axqua.solvers.openfoam.dicts.write_decomposition(case_dir, n_processors)[Quellcode]

Make the case’s own system/decomposeParDict name n_processors domains, for a case that others decompose with that file (the calibration template, which HydroBayesCal copies). A file that names the number already is left as it is. Returns whether the file was written.

Parameter:
Rückgabetyp:

bool

axqua.solvers.openfoam.dicts.write_dicts(of_mesh, cfg, case_dir, *, velocity_cap)[Quellcode]

Write constant/ and system/ (both staged dict sets).

Parameter:
Rückgabetyp:

list[Path]

axqua.solvers.openfoam.dicts.activate(case_dir, stage, cfg=None, *, single=False)[Quellcode]

Make one stage’s dictionaries the active system/ set.

With cfg the three staged files are regenerated from that config before being activated, so a run always honors the configuration it was launched with. Without it they would come from the build, and changing end_time or max_courant and re-running openfoam_run.py would silently run the old values while every message reported the new ones - a trap precisely because nothing about it looks wrong.

The patch names are read back from the written mesh (axqua.solvers.openfoam.polymesh.read_patch_names()), so regenerating needs no remeshing and no state carried over from the build.

Parameter:
Rückgabetyp:

None

Assemble the complete OpenFOAM case: mesh, fields, dictionaries.

The counterpart of axqua.pipeline.run() for the 3D free-surface path. It draws on exactly the same configuration the TELEMAC case was built from - the same ROI, the same liquid boundaries, the same roughness zones, the same discharge and the same outflow rating - so the two models are driven identically and their results can be compared without arguing about the forcing.

What it produces, under <sim_dir>/openfoam/:

0/                alpha.water, U, p_rgh, nut, k, omega   (seeded from r2d.slf)
constant/
  polyMesh/       points, faces, owner, neighbour, boundary
  g, transportProperties, turbulenceProperties
system/
  controlDict, fvSchemes, fvSolution                     (the ACTIVE stage)
  fvOptions, decomposeParDict
  stage1-spinup/, stage2-run/                            (both dict sets)
case.foam                                                (ParaView marker)
class axqua.solvers.openfoam.case.OpenFoamArtifacts(case_dir, polymesh_dir, field_files=<factory>, dict_files=<factory>, mesh=None, report=None, discharge=0.0, inlet_discharges=<factory>, outflow_stage=None, velocity_cap=0.0, hotstart=None, notes=<factory>, cfg=None, state=None)[Quellcode]

Bases: object

Everything the build produced, mirroring axqua.pipeline.Artifacts.

Parameter:
case_dir: Path
polymesh_dir: Path
field_files: list[Path]
dict_files: list[Path]
mesh: object = None
report: object = None
discharge: float = 0.0
inlet_discharges: dict[str, float]
outflow_stage: float | dict[str, float] | None = None
velocity_cap: float = 0.0
hotstart: Path | None = None
notes: list[str]
cfg: object = None
state: object = None
lines()[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.case.RATING_SHARES = (0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 3.0)

Discharges, as shares of the case’s own, at which a derived rating curve is tabulated for an outlet that follows it (openfoam.outflow_boundary: rating-curve). The outlet of the BAW refuses to extrapolate, so the range has to cover what the outflow does while the run settles.

axqua.solvers.openfoam.case.rating_table(cfg, discharge, folder, *, derive=None)[Quellcode]

The stage-discharge relation an outlet follows: [(Q, h), ...] sorted by Q.

boundaries.stage_discharge where it has at least two discharges; otherwise a curve derived from the terrain over RATING_SHARES of discharge, written to folder (one pair is enough for a fixed level, not for a relation).

Parameter:

discharge (float)

Rückgabetyp:

list[tuple[float, float]]

axqua.solvers.openfoam.case.build_case(cfg, *, state=None, hotstart=None, case_dir=None, discharge=None, uniform_bed_ks=False)[Quellcode]

Build the OpenFOAM case for cfg.

state short-circuits reading the 2D result (used by tests); hotstart points at a different r2d.slf than the config’s own; discharge overrides boundaries.prescribed_flowrate. uniform_bed_ks writes one global bed roughness rather than the per-face list, which is what a calibration template needs (see axqua.solvers.openfoam.fields._wall_entries()).

Parameter:
Rückgabetyp:

OpenFoamArtifacts

axqua.solvers.openfoam.case.BUILD_RECORD = 'axqua-build.json'

Written beside case.foam so the build’s own findings survive the terminal. A build-time warning that only reaches a log is the failure this exists to stop: munich-vsf ran 61 hours, balanced its discharge to -0.000%, held Courant at 0.90, and carried no trace of the warning that said the answer could not be trusted.

axqua.solvers.openfoam.case.water_speed(of_mesh, cfg, state=None)[Quellcode]

The reach’s own water speed [m/s]: the 95th percentile of the 2D result.

Parameter:

cfg (Config)

Rückgabetyp:

float

axqua.solvers.openfoam.case.limiting_velocity(of_mesh, cfg, velocity_cap, state=None)[Quellcode]

The velocity that sets the Courant number of the flow.

Under a rigid lid that is the water’s own speed. In the two-phase case it is the air next to the water, which moves a little faster than the water does (AIR_SPEED_FACTOR) and never faster than the limitVelocity cap. It is not the cap, although this function said so until October 2026: the runs it was calibrated on had air at the cap because the outflow drove it there (see axqua.solvers.openfoam.fields.outlet_pressure()), and the estimate built on them was eight times too small once that was corrected.

Parameter:
Rückgabetyp:

float

axqua.solvers.openfoam.case.flow_courant(cfg, courant=None)[Quellcode]

The Courant limit of the flow in the production stage (maxCo).

Parameter:
Rückgabetyp:

float

axqua.solvers.openfoam.case.estimate_time_step(of_mesh, cfg, velocity_cap, *, courant=None, state=None)[Quellcode]

Time step [s] the Courant target will settle on.

Two phases. The limiting length is the thinnest layer in the mesh, not the bed layer and not cell_size: the bed layer is pinned thick enough to clear the roughness, so the layers above it are thinner wherever a column is shallow, and what crosses a layer is the vertical motion at the waterline. Two limits apply, the Courant number of the flow, which the air next to the water reaches first, and that of the water surface; the smaller step is the estimate. Within 12 % of the four measured steps the factors come from.

Rigid lid. No air, and taking the two-phase model across gave an estimate nineteen times too pessimistic (0.019 s against a measured 0.369 s). There the water moves horizontally across the full plan spacing, and only its vertical component - which a terrain-following mesh gets from the bed slope, w ~ u.S - crosses a thin layer. So the two directions are summed on their own terms instead of the worst one being applied to both.

Parameter:
Rückgabetyp:

float

axqua.solvers.openfoam.case.terrain_slope(of_mesh)[Quellcode]

Median bed slope [-] across one plan cell, over the meshed columns.

On a terrain-following mesh this is what turns horizontal flow into vertical flux: water running along a bed that falls by S per cell crosses the layers at roughly u.S. It is read off the structured lattice rather than assumed, because it varies by an order of magnitude between a smooth sand reach and the bar-and-pool bed here (isar-2025: ~0.08).

Rückgabetyp:

float

axqua.solvers.openfoam.case.reach_flush_time(of_mesh, state)[Quellcode]

Rough time [s] for water to traverse the modeled reach once.

The yardstick for „how long must this run be“: a free-surface run cannot be called steady before the domain has been flushed at least once, and comparing that against the configured end_time is the difference between a run that means something and one that merely finished.

Rückgabetyp:

float

axqua.solvers.openfoam.case.cost_lines(of_mesh, state, cfg, velocity_cap)[Quellcode]

What this run will cost, stated before anyone commits to it.

A 3D VOF model of a whole water body is expensive in a way a 2D model is not, and the expense is not obvious from the cell count alone: it is the cell count times the number of time steps, and the time step is set by the thinnest cell. Spelling both out here is the difference between choosing a resolution and discovering one.

Parameter:
Rückgabetyp:

list[str]

axqua.solvers.openfoam.case.estimate_cells(cfg, *, state=None)[Quellcode]

Cell count a build would produce, without building it.

Cheap enough to call before committing to a resolution: it lays the plan lattice and multiplies by the layer count, skipping the DEM sampling and the extrusion.

Parameter:
Rückgabetyp:

int

axqua.solvers.openfoam.case.summarise(art)[Quellcode]

Printable build summary: what was made, and what to run next.

Parameter:

art (OpenFoamArtifacts)

Rückgabetyp:

list[str]

Lauf und Bewertung

Driving the OpenFOAM tools: checkMesh, decomposePar, interFoam.

Mirrors axqua.env.TelemacRuntime and shares its ShellRuntime base: axqua never imports OpenFOAM’s environment, it sources etc/bashrc in a subshell per command. That keeps OpenFOAM’s $PATH, $LD_LIBRARY_PATH and MPI wrappers out of axqua’s own process, where they would collide with TELEMAC’s.

The run is streamed, not captured: an interFoam run on a million cells is hours of work, and a silent terminal gives no way to tell a healthy run from one whose time step has collapsed. OpenFoamProgress renders that as a bar of simulated time with the live Courant number and time step beside it - the two numbers that say whether the air phase is behaving.

class axqua.solvers.openfoam.runtime.OpenFoamProgress(end_time, *, start_time=0.0, out=None, sink=None, echo=True)[Quellcode]

Bases: object

Echo streamed OpenFOAM output with a simulated-time progress bar.

Shares axqua.progress.ProgressBar with the TELEMAC renderer, so a 3D OpenFOAM run and a TELEMAC run look the same in the terminal. The extra quantities shown - Co and dt - are the diagnostic that matters for this solver: a healthy interFoam run holds dt roughly constant near the Courant target, while a run whose air phase is misbehaving shows dt falling by orders of magnitude with the simulated clock barely advancing.

Parameter:
feed(line)[Quellcode]
Parameter:

line (str)

Rückgabetyp:

None

close()[Quellcode]
Rückgabetyp:

None

class axqua.solvers.openfoam.runtime.OpenFoamRuntime(of_config)[Quellcode]

Bases: ShellRuntime

Run OpenFOAM utilities and solvers for a case.

check_available()[Quellcode]

Return the OpenFOAM version string, raising if the environment is broken.

Rückgabetyp:

str

check_mesh(case_dir, *, on_line=None)[Quellcode]

Run checkMesh; return (usable, failed_checks, output).

checkMesh exits 0 even when it fails checks, so the verdict is read from its own report rather than from the exit status.

usable is not „passed every check“. A mesh cut from a real DEM routinely carries a few faces past a skewness or non-orthogonality limit where the survey has a step or a seam - the isar reach has exactly one in 4.9 million, and it comes from the bathymetry, not from the meshing. Treating that as a blocking failure would mean refusing to model the terrain. Only the checks in _FATAL_CHECKS - a cell with no volume, an unclosed cell, broken addressing - actually prevent a run.

Parameter:
Rückgabetyp:

tuple[bool, list[str], str]

decompose(case_dir, *, force=True, n_processors=None)[Quellcode]

decomposePar for n_processors ranks (the case’s own count when not given). Where the case’s system/decomposeParDict names another count, a dictionary of the run’s own is used (dicts.decomposition_for()), so the case’s file is never changed.

Where the processor folders hold results that the case itself does not (a stage run without reconstructing), nothing is split: decomposePar -force would delete them, and whether they are continued or replaced is for the next stage to decide (run_stage()). Returns None then.

Parameter:
reconstruct(case_dir, *, latest=False)[Quellcode]
Parameter:
to_vtk(case_dir, *, latest=False, fields=None)[Quellcode]

Convert the reconstructed case to VTK for QGIS / ParaView.

Part of the runner’s lifecycle, not the plugin’s (plan §13): converting a multi-gigabyte case is exactly the heavy postprocessing QGIS must never do, and it belongs on the machine that ran the solver.

Parameter:
run_stage(case_dir, stage, *, end_time, start_time=0.0, n_processors=None, cfg=None, show_progress=True, sink=None, should_stop=None)[Quellcode]

Activate stage’s dictionaries and run the solver over them.

Pass cfg so the dictionaries are regenerated from the configuration this run was launched with; without it the build-time dictionaries are used and a changed end_time would be reported but not honored (see axqua.solvers.openfoam.dicts.activate()).

Returns the CompletedProcess; a non-zero exit is returned rather than raised, so a caller running both stages can report which one failed.

The case is made ready for the stage first (_prepare()): a stage that starts afresh moves what earlier runs left into previous-run/ and is split anew for its ranks, a stage that continues runs on what the stage before it left, also where only the processor folders hold it. A parallel run gets the decomposition dictionary of its ranks (a parallel solver reads it as well and refuses one that names more domains than it has ranks).

Parameter:
latest_time(case_dir)[Quellcode]

Largest numeric time directory in the case (0 when only 0/ exists).

Parameter:

case_dir (str | Path)

Rückgabetyp:

float

Has the OpenFOAM run reached a steady discharge balance?

The 3D analogue of axqua.flux_convergence, and deliberately the same judgment: a free-surface run is converged when the water entering equals the water leaving, measured as a relative imbalance against the prescribed inflow. The tolerance grades carry over unchanged (hydrodynamics.flux_tolerance, 1e-3 for a result read as discharge/depth/velocity), so a 2D and a 3D run of the same reach are held to the same standard and the numbers are directly comparable.

The data come from the surfaceFieldValue monitors axqua.solvers.openfoam.dicts._flux_functions() writes into controlDict: sum(phi) weighted by alpha.water over each inlet and outlet patch, i.e. the water discharge alone rather than the water-plus-air flux phi would give. OpenFOAM writes one surfaceFieldValue.dat per monitor under postProcessing/<name>/<startTime>/.

Sign convention: phi is positive out of the domain, so an inlet reports a negative value. This module reports magnitudes, as the TELEMAC side does.

axqua.solvers.openfoam.report.DELIVERY_TOLERANCE = 0.02

How far the WATER actually entering may fall short of the prescribed discharge before it is a finding. variableHeightFlowRateInletVelocity, which the case used until October 2026, sets the velocity proportional to alpha, so the MIXTURE flux integrates to the prescription exactly while the water flux integrates to Q * <alpha^2>/<alpha> - short by whatever the interface is smeared by. On munich-vsf that ratio was 0.821, steady from t=10 to t=97: the case believed 0.135 m3/s and delivered 0.111. Nothing in the output said so, because the balance only ever compared the inflow against the OUTflow, and both were equally short. The inlet prescribes the mass flow now (see fields) and delivers within 0.3 %, and the check stays: 2% is tight on purpose - this is the one number a calibration cannot absorb.

class axqua.solvers.openfoam.report.DischargeHistory(time, per_patch=<factory>, inflow=None, outflow=None, target=0.0, tolerance=0.001, steady_time=None, converged=False, rigid_lid=False, top=None)[Quellcode]

Bases: object

Per-patch water discharge against time, and the balance between them.

Parameter:

  • time (np.ndarray)

  • per_patch (dict[str, np.ndarray])

  • inflow (np.ndarray | None)

  • outflow (np.ndarray | None)

  • target (float)

  • tolerance (float)

  • steady_time (float | None)

  • converged (bool)

  • rigid_lid (bool)

  • top (np.ndarray | None)

time: np.ndarray
per_patch: dict[str, np.ndarray]
inflow: np.ndarray | None = None
outflow: np.ndarray | None = None
target: float = 0.0
tolerance: float = 0.001
steady_time: float | None = None
converged: bool = False
rigid_lid: bool = False

True when the case runs under a rigid lid. The balance is then STRUCTURAL, not evidence of anything - see lines().

top: np.ndarray | None = None

Water leaving through the open top (the lidLeak monitor), positive out; None for a run without that monitor. It is part of the balance: a run that loses water through the top is steady when inflow = outflow + that loss, and whether the loss is acceptable is the surface verdict’s question.

property imbalance: numpy.ndarray

Relative |Q_in| - |Q_out| over |Q_in|, the same measure as in 2D.

Less the water that leaves through the top, where the run has that monitor.

property final_imbalance: float | None

|relative imbalance| the run settled at - the number a study minimizes.

The mean of the last STEADY_WINDOW finite samples, not the single last one: the same window already decides convergence, and one final sample on a still-oscillating run is noise rather than a settled value.

None when there is nothing to read, which is not the same as zero - a caller ranking runs must be able to tell „balanced“ from „never reported“.

property delivered: float | None

Water actually entering, over the discharge the case prescribes.

The inflow monitor is alpha-weighted, so it is the WATER flux; the prescription is handed to a boundary condition that controls the MIXTURE. Those are not the same number whenever the inlet interface is smeared over more than nothing, and the gap is a straight bias on every result.

lines()[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.report.read_monitors(case_dir)[Quellcode]

Every Q_* monitor in the case, concatenated across restart directories.

A two-stage run writes postProcessing/Q_inlet_1/0/ for the spin-up and .../30/ for the production stage; both are read and joined in time order so the history spans the whole run rather than only its last leg.

Parameter:

case_dir (str | Path)

Rückgabetyp:

dict[str, tuple[numpy.ndarray, numpy.ndarray]]

axqua.solvers.openfoam.report.analyse(cfg, case_dir=None, *, target=None)[Quellcode]

Build the discharge history and decide whether the run has settled.

Parameter:
Rückgabetyp:

DischargeHistory

axqua.solvers.openfoam.report.LID_LEAK_TOLERANCE = 0.01

Water leaving through the lid, as a fraction of the prescribed discharge, above which the run is not a result. Area is the wrong scale for a breach and this is the right one: on munich-vsf the lid was wet over 0.47% of its area - a figure that reads as noise - while venting FIFTEEN TIMES the discharge. The pressure outlet silently replaced what left, so the domain ran a steady spurious circulation and every other diagnostic stayed healthy. 1% is already generous: a lid is a boundary water has no business crossing at all. Judged on the mean over the second half of the run, not on one sample: the Isar example at 1 m lost 0.8 % on average with single samples of 1.7 %.

class axqua.solvers.openfoam.report.SurfaceFreedom(lid_area=0.0, wall_area=0.0, lid_total=0.0, wall_total=0.0, lid_leak=0.0, lid_leak_mean=0.0, discharge=0.0, prescribed=False, measured=False, settled=True)[Quellcode]

Bases: object

Whether the free surface was ever stopped by the mesh rather than the flow.

Parameter:
lid_area: float = 0.0
wall_area: float = 0.0
lid_total: float = 0.0
wall_total: float = 0.0
lid_leak: float = 0.0
lid_leak_mean: float = 0.0
discharge: float = 0.0
prescribed: bool = False
measured: bool = False
settled: bool = True

False while the run has not yet passed its spin-up. The monitors exist and are being written, but every sample so far belongs to a transient that is violent by design, so there is no verdict to give. Distinguished from free deliberately: „nothing to report yet“ and „the surface was free“ look identical in a dataclass of zeros, and one of them is a false all-clear.

property lid_fraction: float
property wall_fraction: float
property leak_fraction: float

Water lost through the lid in the second half of the run, over the discharge the case prescribes.

property peak_leak_fraction: float
property touched_lid: bool

the freeboard was too small.

Type:

The surface itself reached the lid

property leaking: bool
property hit_lid: bool
property hit_wall: bool
property free: bool
lines(cfg)[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.report.surface_freedom(cfg, case_dir=None)[Quellcode]

Read the lidContact / wallContact monitors into a verdict.

The seed fixes two things the flow has no say in - how high the lid stands and how far the footprint reaches - and both are sized from a TELEMAC result that is itself a model. If the 3D surface reached either, the answer was set by a meshing decision rather than by the hydraulics, and nothing else in the output would say so: the run converges, the discharge balances, the picture looks like a river.

Zero on both is the result you want. Reported per run rather than checked at build time, because it is the only moment the question can actually be answered.

Parameter:

case_dir (str | Path | None)

Rückgabetyp:

SurfaceFreedom

axqua.solvers.openfoam.report.LID_STEP_WARN = 0.5

p99 lid step, as a fraction of the local depth, above which the mode is doubtful: half the water depth lost within two cells is a drop, not a slope. A heuristic, which is why the build warns here rather than refusing.

axqua.solvers.openfoam.report.LID_STEP_SEVERE = 1.0

…and above this the prescribed surface steps further than the water is deep - a weir, a drop structure, or the slots of a fish pass. The build REFUSES here unless openfoam.allow_stepped_lid is set; the verdict reaches the result either way.

Both live in this module, not in mesh, so that reading a finished run’s verdict costs nothing heavier than numpy - and so the build-time grade and the reported grade are the same two numbers rather than two copies that can drift.

class axqua.solvers.openfoam.report.LidApplicability(rigid_lid=False, median=None, p99=None, recorded=False)[Quellcode]

Bases: object

The build’s rigid-lid verdict, carried forward to the finished run.

A rigid lid is a slip wall on the prescribed free surface, so the applicability question - does that surface step sharply on the scale of a cell? - can only be asked of the input, at build time. The result cannot answer it afterwards: a run whose lid is wrong converges, balances its discharge and looks like a river. munich-vsf did exactly that for 61 hours, twice, and in both cases the output carried no trace of the build-time warning.

So the numbers are written into the case at build time (axqua.solvers.openfoam.case.BUILD_RECORD) and read back here, which is the same role SurfaceFreedom plays for the two-phase case: a statement about whether the answer was set by a meshing decision rather than by the flow.

Parameter:
rigid_lid: bool = False
median: float | None = None
p99: float | None = None
recorded: bool = False
property doubtful: bool
property severe: bool
lines(cfg)[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.report.lid_applicability(cfg, case_dir=None)[Quellcode]

Read the build record’s lid-step test back out of a built case.

Parameter:

case_dir (str | Path | None)

Rückgabetyp:

LidApplicability

axqua.solvers.openfoam.report.verdict_lines(cfg, case_dir=None, history=None)[Quellcode]

Everything a finished run has to be read against, in one place.

Parameter:
Rückgabetyp:

list[str]

axqua.solvers.openfoam.report.write_report(cfg, case_dir=None, *, out_dir=None, target=None)[Quellcode]

Analyze the run and write discharge-convergence.csv + .png.

Parameter:
Rückgabetyp:

tuple[DischargeHistory, list[Path]]

Kalibrierung

Surrogate-assisted Bayesian calibration of the OpenFOAM case.

The OpenFOAM half of aXqua’s HydroBayesCal integration; the TELEMAC half is axqua.bayescal, and the plumbing both share is axqua.hbc. What is calibrated here is the bed roughness ks and, optionally, the k-epsilon closure coefficients, against velocity components measured at points - U_x, U_y, U_z (and TKE) - which is what an ADV campaign such as a FlowTracker survey actually produces, and what a depth-averaged 2D model cannot be calibrated against at all.

Run the campaign coarse and under a rigid lid. mode: rigid-lid removes the air phase by construction (about 90% of a two-phase case’s cells) and fits the layers to the water alone, and cell_size_factor coarsens the plan lattice; a surrogate needs tens of runs, and a production-resolution two-phase case cannot supply them. campaign_config() builds that variant from the case’s own config without disturbing it. Two consequences worth being explicit about:

  • the free surface is prescribed from the 2D result rather than solved, so any roughness error that would have shown up as a surface-slope error is absorbed into the velocity field instead. The posterior should be verified once at production resolution with mode: vof before it is believed.

  • rigid-lid has exactly one run stage (axqua.solvers.openfoam.dicts.stages()), so the case HydroBayesCal copies per run is correct straight out of build_case with no stage to activate incorrectly.

One global ks. HydroBayesCal’s OpenFOAM binding writes a single Ks uniform to the bed patch. There is no per-zone roughness here, so this posterior is not comparable with the zoned zone<N> posterior of a TELEMAC calibration on the same reach, and the two must not be reported as if they were.

axqua.solvers.openfoam.calibration.CAMPAIGN_SUBDIR = 'openfoam'

Sub-directory of calibration_dir holding everything this module produces.

axqua.solvers.openfoam.calibration.VERTICAL_ERROR_FLOOR = 0.03

Error floor [m/s] on the vertical velocity component. Deliberately coarser than the horizontal floor: a side-looking ADV’s w at 0.6*h is dominated by mounting tilt and by the probe’s own beam geometry, so the raw per-point error claims a precision the measurement does not have. This keeps U_z in the likelihood as a weak constraint on the secondary circulation instead of as the loudest voice in it. It is a judgment, not a measurement - drop U_z from calibration_quantities (it is still extracted) if that trade is not wanted.

axqua.solvers.openfoam.calibration.campaign_config(cfg, *, cell_size_factor, end_time, write_interval, n_processors, n_avg_timesteps=5, turbulence='kEpsilon', mode='rigid-lid', n_layers=None, spinup_time=None, air_velocity_cap=None, min_column_height=None, auto_bed_layer=None)[Quellcode]

A copy of cfg describing the campaign case, not the production case.

A deep copy rather than a mutation: the caller’s config still describes the full-resolution two-phase case that the verification run rebuilds, so the two never have to be reconciled and running a campaign cannot quietly change what openfoam_preprocessing.py would build afterwards.

mode defaults to rigid-lid, which is what makes a surrogate affordable: about 90% of a two-phase case’s cells are air and a campaign needs tens of runs. Pass vof where the rigid lid is not admissible - it is a slip WALL, so wherever the prescribed surface steps it converts head into velocity, and on a shallow reach that is most of the domain. The remaining overrides exist because a two-phase campaign has to be cheapened along different axes than a rigid-lid one (fewer layers, a tighter air cap) and those are exactly the numbers a production config should not carry.

Parameter:

  • cfg (Config)

  • cell_size_factor (float)

  • end_time (float)

  • write_interval (float)

  • n_processors (int)

  • n_avg_timesteps (int)

  • turbulence (str)

  • mode (str)

  • n_layers (int | None)

  • spinup_time (float | None)

  • air_velocity_cap (float | None)

  • min_column_height (float | None)

  • auto_bed_layer (bool | None)

Rückgabetyp:

Config

axqua.solvers.openfoam.calibration.stage_case_template(cfg, *, state=None, rebuild=False, parameters=(), n_avg_timesteps=5)[Quellcode]

Build (or reuse) the case HydroBayesCal copies once per calibration run.

Parameter:
Rückgabetyp:

Path

axqua.solvers.openfoam.calibration.openfoam_parameters(cfg, *, names=None)[Quellcode]

The OpenFOAM-legal subset of the case’s declared calibration parameters.

axqua.calibration.merged_parameters() also collects TELEMAC roughness zones and .cas keywords, every one of which the OpenFOAM binding rejects - and rejects after the experimental design has been sampled and a case copied. Filtering here, with a warning naming each dropped parameter, is the difference between a mis-declared parameter costing a log line and costing a setup.

Parameter:
Rückgabetyp:

list

axqua.solvers.openfoam.calibration.check_turbulence(cfg, parameters)[Quellcode]

Refuse a k-epsilon coefficient on a case that does not solve k-epsilon.

Checked here rather than in OpenFoam.validate(), which runs on every load_config and cannot see the calibration block. The failure it prevents is specific and expensive: with kOmegaSST the case carries no kEpsilonCoeffs, so HydroBayesCal raises on the first perturbed run - after the design has been sampled.

Parameter:
Rückgabetyp:

None

axqua.solvers.openfoam.calibration.emit_openfoam_config(cfg, *, template, csv, parameters, calibration_quantities, extraction_quantities, out_dir, n_avg_timesteps=5, n_processors=8, init_runs=None, max_runs=None, only_bal_mode=False)[Quellcode]

Write config_OpenFOAM.py for the staged bal_openfoam.py.

Deliberately not a parameterisation of axqua.calibration.emit_hbc_config(): the OpenFOAM schema uses simulation where TELEMAC uses hydrodynamic_simulation, puts n_cpus inside sampling, carries an extra interfoam block, and treats case_template_dir as a directory copied per run rather than referenced in place.

Parameter:
Rückgabetyp:

Path

axqua.solvers.openfoam.calibration.stage_openfoam_driver(out_dir, *, checkout=None)[Quellcode]

Stage bal_openfoam.py beside the emitted config.

Parameter:
Rückgabetyp:

Path

axqua.solvers.openfoam.calibration.run_openfoam_calibration(cfg, *, calibration_quantities=('U_x', 'U_y', 'U_z'), extraction_quantities=('U_x', 'U_y', 'U_z', 'TKE'), parameters=None, cell_size_factor=4.0, end_time=600.0, write_interval=20.0, n_avg_timesteps=5, n_processors=None, init_runs=None, max_runs=None, vel_err_floor=0.01, vertical_err_floor=0.03, launch_mode='prepare', rebuild=False, checkout=None, state=None, force=False, ctx=None)[Quellcode]

Calibrate the OpenFOAM case against measured velocity components.

launch_mode:

prepare

build the template, write the CSV and config, stage the driver, print the command and stop. Nothing is run.

smoke

the same, then a tiny isolated 3-run campaign on a short template. This is the gate to run before committing solver hours: it exercises the whole chain - parameter routing into the dictionaries, the solver launch, the field extraction, the per-run accumulation - in minutes.

run

launch the full campaign.

resume

launch with only_bal_mode, reusing a completed initial design.

Parameter:
Rückgabetyp:

int

axqua.solvers.openfoam.calibration.read_posterior(campaign_dir, *, statistic='map')[Quellcode]

The calibrated parameter values from a finished campaign.

Reads BAL_dictionary.pkl, which carries the posterior sample of every BAL iteration alongside calibration_parameters and param_values - so the file is self-describing and the parameter order cannot be mismatched here.

statistic is map (the highest-density point of the posterior sample, approximated by the sample mode over a histogram) or mean. The MAP is the default because a roughness posterior is often skewed, and its mean can sit where the sample has little mass.

Parameter:

  • campaign_dir (Path)

  • statistic (str)

Rückgabetyp:

dict[str, float]

axqua.solvers.openfoam.calibration.apply_posterior(cfg, values, *, bound_margin=0.05)[Quellcode]

Write calibrated values into cfg, and report what changed.

Refuses when a calibrated value sits within bound_margin of its own prior bound: that is the signature of an answer lying outside the range that was allowed, and adopting it silently would present a clipped estimate as a result.

Parameter:
Rückgabetyp:

list[str]

axqua.solvers.openfoam.calibration.model_relative_z_resolver(campaign, *, out_dir, margin=0.05, min_depth=0.05, cache=None, state=None)[Quellcode]

A z_resolver that puts each target at its measured RELATIVE depth.

z = bed_model(x, y) + f * depth_model(x, y), with f derived from the measurement (axqua.ground_truth.relative_height()). Nothing but x, y and a ratio comes from the survey, so no vertical survey error can reach the target - which is the whole point, since on this reach the survey carried an un-subtracted rover-pole height, the measurement height was never added, and the DEM has a depth-proportional bias.

Writes target-placement.csv beside the calibration inputs so every target’s position is auditable rather than implicit.

Parameter:

Einphasige freie Oberfläche

Leg C: a free surface with no air phase at all (potentialFreeSurfaceFoam).

Why it is worth the trouble on this case. The two-phase runs spend four cells in five on air: munich-vsf leg A is 582,924 cells of which 22% hold water, and the time step is set by the air’s velocity cap, not by the water. lww-134 measured the cost of making that case correct - doubling the freeboard so the lid stops clipping the surface, and adding the layers that go with it - at 3.3x the per-step cost for 1.5x the cells. Resolving the air better is the expensive direction. This solver deletes it.

HOW IT WORKS. The domain is water only, as under aXqua’s rigid lid, and the mesh top is a real boundary. The free surface is not meshed: it is carried as a displacement field zeta on that boundary, entering the pressure equation through waveSurfacePressure as p_gh = -(g & zeta). So the surface moves while the mesh does not - which is the whole difference from mode: rigid-lid, where the surface is frozen at whatever the 2D seed said.

THE TRAP, AND WHY 0/zeta IS WRITTEN. p_gh is an ABSOLUTE gauge: createFields sets p_gh = p - (g & mesh.C()), so hydrostatic equilibrium with the surface at z_s needs p_gh = |g| z_s on the top. The boundary condition instead imposes p_gh = |g| zeta_z, and zeta defaults to zero - which asserts the undisturbed surface lies at z = 0. On a wave tank or a ship hull that is true by construction and nobody notices. On this fish pass the surface falls 1.4 m over 25 m, so the default datum would impose a pressure error of |g| * 1.4 = 13.7 m2/s2 along the pass and drive a flow that is pure gauge.

zeta is READ_IF_PRESENT, so the datum is fixable without touching OpenFOAM: seed its boundary values on the top patch with the local lid elevation. zeta then means „where the surface actually is“ rather than „how far it has moved from z=0“, the linearisation is about the 2D surface, and the quantity that has to stay small is the 3D correction to it rather than the whole fall. Verified on the v2406 oscillatingBox tutorial, which ships no 0/zeta: seeding one reports min/max zetap = 0.5, 0.5 on the first step and evolves from there, so the value is read and used rather than reset.

WHAT IS STILL AN APPROXIMATION, stated rather than buried. The surface condition is linearised: the pressure is applied on the mesh top, not on the displaced surface, so it is only right while |zeta - z_lid| is small against the depth. A plunging jet, an overturning wave and air entrainment are all outside it, and a vertical-slot pass has jets through every slot. That is a question for the measurement, not for this docstring - which is why leg C is a leg of a comparison rather than a replacement.

axqua.solvers.openfoam.potential.TOP_PATCH = 'lid'

Patch name the water-only mesh gives its top (see mesh.LID_PATCH).

axqua.solvers.openfoam.potential.is_potential(cfg)[Quellcode]

True when this case is leg C rather than a two-phase or rigid-lid run.

Rückgabetyp:

bool

axqua.solvers.openfoam.potential.STEEP_FACE_FRACTION = 0.005

Fraction of top faces that may be steeper than 45 degrees before the surface condition stops meaning anything. waveSurfacePressure advances the surface as zeta += dt * nf * phi/|Sf| - ALONG THE FACE NORMAL - and then imposes p_gh = -(g & zeta), which reads only the vertical component. On a near-horizontal top face those are the same thing. On a steep one the displacement is mostly HORIZONTAL, which is not a free-surface motion at all, and the vertical part that survives is whatever the slope leaves behind.

axqua.solvers.openfoam.potential.check_lid_is_flat_enough(of_mesh)[Quellcode]

(steep fraction, flattest |n_z|) for the top patch, and log the verdict.

Measured on munich-vsf, whose lid follows a 14-slot stepped 2D surface: 4.1% of top faces steeper than 45 degrees, the worst at |n_z| = 0.068 - 86 degrees from horizontal, effectively a vertical wall being asked to behave as a free surface. The consequence was not subtle once found: zeta collapsed from its seeded 0.33 m to -2.9 m within one second of simulated time, monotonically, on a column 0.2 m deep. The giveaway was that the first reported min/max zetap was 0.0479 in every run regardless of what else changed, and z * |n_z| over the lid reproduces 0.0479 / 2.7512 exactly - so what the BC reports, and works with, is the projection onto the face normal.

Rückgabetyp:

tuple[float, float]

axqua.solvers.openfoam.potential.check_applicable(cfg)[Quellcode]

Refuse the combinations that would run and be wrong.

A two-phase mesh has an air region this solver has no equation for, and it would fill with water and report a plausible answer. Better to stop here.

Rückgabetyp:

None

axqua.solvers.openfoam.potential.lid_elevation_by_column(of_mesh)[Quellcode]

Free-surface elevation per plan COLUMN, indexed by column id.

Built from the top faces themselves, via each face’s owner cell, so it needs no assumption about patch face ordering. The first version indexed the lid patch’s face list as though it were in column order; that happens to be true for this mesher, but it is not stated anywhere and p, zeta and the outlet condition all depend on it agreeing. Reading the owner is one array lookup and removes the question.

Not of_mesh.lid (per plan VERTEX) and not column_wse (the 2D seed, which the mesher may have clipped): the reference has to be the elevation of the very faces waveSurfacePressure acts on.

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.potential.top_face_elevation(of_mesh)[Quellcode]

Top-patch face-centre elevations, in the patch’s own face order.

This is what 0/zeta needs: the boundary value list is read in patch order.

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.potential.surface_pgh_on_patch(of_mesh, patch, g)[Quellcode]

p_gh per face on a lateral patch, from the free surface above it.

p_gh = |g| z_s IS the statement „the free surface here is at z_s“, so a lateral boundary that wants a prescribed tailwater has to say it face by face. Writing ONE value over the whole patch asserts a single level across a face whose own lid varies - on munich-vsf the outlet spans z 0.321 to 0.823 - and the mismatch is a sustained head error that drains the domain. That is what collapsed zeta from 0.33 to -3.0 m in the first second of the first run that got as far as running.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.potential.hydrostatic_pressure(of_mesh, g=9.81)[Quellcode]

Kinematic p whose p_gh puts the surface on the mesh top.

createFields recomputes p_gh = p - (g & C) from the p on disk, so p is what actually sets the initial state - writing p_gh values alone would be silently overwritten and only its boundary TYPES would survive. p = |g| (z_lid - z) gives p_gh = |g| z_lid: hydrostatic, with the surface exactly where the 2D seed put it.

Parameter:

g (float)

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.potential.write_potential_fields(of_mesh, cfg, case_dir, *, state=None)[Quellcode]

0/ for leg C. Returns the files written.

Rückgabetyp:

list[Path]

axqua.solvers.openfoam.potential.transport_properties(cfg)[Quellcode]

Single phase: one kinematic viscosity, no phase pair, no surface tension.

Rückgabetyp:

str

axqua.solvers.openfoam.potential.fv_schemes()[Quellcode]

No alpha equation, so no compression scheme and no MULES limiter.

ddtSchemes default Euler matters more than it looks: waveSurfacePressure reads the ddt scheme registered for zeta to decide how to integrate the surface displacement, and aborts on anything it does not recognize.

Rückgabetyp:

str

axqua.solvers.openfoam.potential.fv_solution(cfg)[Quellcode]

p_gh in place of p_rgh, and nothing at all for alpha.

Rückgabetyp:

str

Feinsediment über einem Festbett

Fine sediment over a bed that stays in place, with Olsen’s sediDriftFoam2.

What is computed

sediDriftFoam2 (Nils Reidar B. Olsen, NTNU) is simpleFoam with one more equation: the concentration of one grain size, which is carried with the flow, settles with its fall velocity and is exchanged with the bed. The flow is steady and has one phase under a fixed water surface, so the case is built on the mesh of openfoam.mode: rigid-lid whatever mode the free-surface case of the reach uses.

The solver can also move the bed and the water surface. aXqua does not use that: the number of iterations before the first bed change (initialIterations) is set beyond the end of the run. Three rules in the source of the solver are the reason, and each of them would change a natural reach in a way that has nothing to do with sediment. A bed point is lowered until the water is ten bottom cells deep. The points of a vertical are spread evenly after every bed change, which undoes the bottom layer described below. And the inflow and the last meter before the largest x coordinate are treated as the ends of a laboratory flume.

What comes out is therefore the sediment budget of a bed that stays where it is: the concentration at the bed, the rate at which grains settle onto it, the rate at which the flow picks them up, and the difference of the two, per bed area. Over a period that is a mass of fine sediment the bed takes up or gives off, and for a bed layer of given thickness a change of its share of fines (evaluate()).

What the solver requires of a mesh, and how it is checked

The requirements are read from sediDriftFoam2.C and verified on the mesh that is written (check_columns()), because the solver itself checks none of them:

  • patches named bedWall and freeSurface whose face k belong to the same column of hexahedra;

  • bed faces with one vertex order throughout, so that the solver finds the neighbours of a bed face by comparing vertex numbers;

  • a positive height of the first cell center above the bed face center, in the solver’s own definition mesh.C().z() - patch.Cf().z(). It enters a logarithm and ends the run with a floating point exception where it is zero or negative. On a terrain-following mesh that happens wherever neighboring verticals have layers of different thickness on a sloping bed, since the cell center then lies towards the thick side, which is the low side. The layer at the bed is therefore given one thickness everywhere (pin_bottom_layer()). A steep and twisted bed face still shortens the distance, so it is computed for every face with OpenFOAM’s own definitions of the two centers (axqua.solvers.openfoam.quality) and a mesh with a face at or below the roughness limit of the logarithm is refused.

How deposition and pick-up are booked

Read from the source as well, because the budget in evaluate() has to be the one the solver solves. Grains leave the water through the convection term of the bed cell (the bounded scheme turns the inflow of the fall velocity into a sink), which makes the settling rate fall velocity x concentration in the bed cell x plan area. Pick-up is a volume source in the same cell, fall velocity x reference concentration / (2 x height of the cell center) x cell volume. Both are evaluated here per plan area of the column, and the difference of what enters and what leaves the reach is compared with their sum over the bed.

axqua.solvers.openfoam.sediment.GRAIN_DENSITY = 2650.0

the formulas of the solver carry the submerged density 1650 kg/m3 and the density of water 1000 kg/m3 as numbers.

Type:

Density of the grains [kg/m3]. Not a setting

axqua.solvers.openfoam.sediment.VISCOSITY = 1e-06

Kinematic viscosity of water [m2/s], which the formulas of the solver assume.

axqua.solvers.openfoam.sediment.PARTICLE_PARAMETER = 25296.0

the grain size times this number is van Rijn’s particle parameter.

Type:

(1.65 g / nu^2)^(1/3) [1/m] as the solver writes it

axqua.solvers.openfoam.sediment.TAN_REPOSE = 0.7

Tangent of the angle of repose in Dey’s correction of the critical shear stress on a sloping bed, the smallest factor of that correction, and the cosine used by the sand slide. The values of Olsen’s own example; the sand slide never runs here.

axqua.solvers.openfoam.sediment.NEVER = 2000000000

no bed change before this iteration, so none at all.

Type:

An int of the solver

axqua.solvers.openfoam.sediment.LOG_LIMIT = 1.0

The argument of the solver’s logarithm, 30 x height / roughness. At 1 the shear velocity is infinite; below e it exceeds 0.4 times the velocity it is derived from.

axqua.solvers.openfoam.sediment.fall_velocity(grain_size)[Quellcode]

Settling velocity of a natural grain in still water [m/s].

After Ferguson and Church (2004), with their constants for natural sand of sieve diameter (18 and 1.0): one formula that is Stokes‘ law for silt and the square-root law for coarse sand, and rises with the grain size throughout. The three formulas of van Rijn (1984), which go with the transport formulas of the solver, were tried first and not kept: where Stokes‘ law hands over at 0.1 mm they step down by 13 %, so that a slightly coarser grain would settle more slowly. The two agree within 10 % over the range, and both give the 0.2 m/s that Olsen’s example states for 2 mm.

Parameter:

grain_size (float)

Rückgabetyp:

float

axqua.solvers.openfoam.sediment.settling_velocity(cfg)[Quellcode]

The fall velocity of the case [m/s], given or computed.

Rückgabetyp:

float

axqua.solvers.openfoam.sediment.inflow_concentrations(cfg, inlets)[Quellcode]

Concentration by mass [kg/m3] at every inlet patch.

One number serves all of them. A mapping names the lines, as inlet-1 or as the number alone, and has to name each of them: a line left out would get clear water without anybody having said so.

Parameter:

inlets (list[str])

Rückgabetyp:

dict[str, float]

axqua.solvers.openfoam.sediment.sediment_config(cfg)[Quellcode]

A copy of cfg whose OpenFOAM mesh is the one this solver needs.

Water only under a lid on the surface of the 2D result, whatever the free-surface case of the reach is set to, and judged by the rules of the rigid lid.

axqua.solvers.openfoam.sediment.pin_bottom_layer(of_mesh, thickness)[Quellcode]

Give the layer at the bed one thickness at every vertical; returns it.

Changes the levels of the points and nothing else. The lattice, the faces and the numbering stay as they were built, which is why the solver’s rules about them, verified for the mesher, still hold afterwards.

Parameter:

thickness (float)

Rückgabetyp:

float

axqua.solvers.openfoam.sediment.rename_patches(polymesh)[Quellcode]

bed and lid under the names the solver looks up.

Rückgabetyp:

None

axqua.solvers.openfoam.sediment.bed_neighbours(bed_faces)[Quellcode]

(east, west, north, south, count) of every bed face, as the solver finds them.

sediDriftFoam2.C compares vertex numbers: the east neighbour of a face has its vertices 0 and 1 where the face has 3 and 2, the north neighbour its vertices 0 and 3 where the face has 1 and 2. A face without a neighbour keeps 0 there, as the arrays of the solver do, and count says how many it has.

Parameter:

bed_faces (numpy.ndarray)

Rückgabetyp:

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]

axqua.solvers.openfoam.sediment.bed_vectors(points, bed_faces, centres, neighbours)[Quellcode]

Unit vector normal to the bed at every bed face, as the solver computes it.

From the centers of the four neighboring faces where a face has all four (computevx3.h), from its own corners otherwise (computevx4.h).

Parameter:
Rückgabetyp:

numpy.ndarray

class axqua.solvers.openfoam.sediment.BedGeometry(cells, corners, centre, plan_area, height, volume, vector, depth, roughness_zone, outlets=<factory>, inlets=<factory>)[Quellcode]

Bases: object

What the solver sees of the bed, per bed face in the order of the patch.

Parameter:
cells: numpy.ndarray
corners: numpy.ndarray
centre: numpy.ndarray
plan_area: numpy.ndarray
height: numpy.ndarray
volume: numpy.ndarray
vector: numpy.ndarray
depth: numpy.ndarray
roughness_zone: numpy.ndarray
outlets: dict[str, numpy.ndarray]
inlets: dict[str, float]
property n: int
save(path)[Quellcode]
Parameter:

path (Path)

Rückgabetyp:

Path

classmethod load(path)[Quellcode]
Parameter:

path (Path)

Rückgabetyp:

BedGeometry

axqua.solvers.openfoam.sediment.check_columns(of_mesh)[Quellcode]

Verify what the solver assumes about the mesh; return what it sees of the bed.

Raises MeshError for a mesh the solver would compute nonsense on without saying so. The height of the first cell center is left to check_heights(), which knows the roughness.

Rückgabetyp:

BedGeometry

axqua.solvers.openfoam.sediment.required_bottom_layer(geometry, roughness, bottom)[Quellcode]

Thickness of the bottom layer [m] with which every bed face suits the solver.

The height of the first cell center above a bed face is half the bottom layer plus an offset of that face, which does not depend on the layer: the cell is the bed face moved up, so its center lies over the middle of the face in plan, while the center of the face lies towards its steeper part. On a steep and twisted face the offset is negative (-0.026 m on the worst face of the Isar example). This returns the layer with which the worst face stands clear of the solver’s logarithm by the factor e, at which the shear velocity it computes is 0.4 times the flow velocity.

Parameter:
Rückgabetyp:

float

axqua.solvers.openfoam.sediment.check_heights(geometry, roughness, bottom, *, largest=None)[Quellcode]

Refuse a mesh on which the solver’s logarithm fails; return what to note.

The solver computes the shear velocity as 0.4 U / ln(30 h / ks) with h the height of the cell center above the bed face center. largest is the thickest bottom layer the mesh allows, which decides what the refusal advises.

Parameter:
Rückgabetyp:

list[str]

class axqua.solvers.openfoam.sediment.SedimentArtifacts(case_dir, mesh=None, geometry=None, record=<factory>, notes=<factory>)[Quellcode]

Bases: object

What the build of the fine sediment case produced.

Parameter:
case_dir: Path
mesh: object = None
geometry: BedGeometry | None = None
record: dict[str, Any]
notes: list[str]
lines()[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.sediment.write_fields(of_mesh, cfg, case_dir, *, state, discharges, concentrations, roughness)[Quellcode]

Write the nine fields of 0/.

concentrations are by mass [kg/m3]; the solver carries the share of the volume. The flow starts from the velocity of the 2D result, the same in all layers, and the concentration from the mean of what enters.

The bed has one roughness, roughness, in the wall function of the flow as in the shear stress on the grains. The free-surface case gives every bed face the roughness of its zone, and that does not work here: on the Isar the zones reach 0.5 m, ten times the bottom layer, and with them the residuals of the pressure, the vertical velocity and the turbulence stay at 5e-3 from the fiftieth iteration on, where one roughness lets them fall below 1e-5 within 400 iterations.

Parameter:
Rückgabetyp:

list[Path]

axqua.solvers.openfoam.sediment.RESIDUAL = 1e-05

Initial residual of every equation below which the run counts as converged.

axqua.solvers.openfoam.sediment.sediment_properties(cfg, *, roughness, velocity)[Quellcode]

constant/sedimentProperties.

The first block is the sediment of the case. The second keeps the bed and the water surface where they are; the third is read by the solver and has no effect without a bed change.

Parameter:
Rückgabetyp:

str

axqua.solvers.openfoam.sediment.write_dicts(cfg, case_dir, *, roughness, velocity, iterations)[Quellcode]
Parameter:
Rückgabetyp:

list[Path]

axqua.solvers.openfoam.sediment.build_case(cfg, *, state=None, hotstart=None, case_dir=None, discharge=None)[Quellcode]

Build the fine sediment case of cfg from the 2D result of the reach.

state is the 2D result as the OpenFOAM mesher reads it; without one the case’s own steady result is read, and without that there is no water surface to put the lid on and the build says so.

Parameter:
Rückgabetyp:

SedimentArtifacts

axqua.solvers.openfoam.sediment.is_built(case_dir)[Quellcode]
Parameter:

case_dir (Path)

Rückgabetyp:

bool

axqua.solvers.openfoam.sediment.latest_time(case_dir)[Quellcode]

Name of the last time folder a run wrote, None before the first.

Parameter:

case_dir (Path)

Rückgabetyp:

str | None

axqua.solvers.openfoam.sediment.read_internal(path, n_cells)[Quellcode]

The cell values of an ASCII field file, scalars or vectors.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.sediment.read_patch(path, patch, n_faces)[Quellcode]

The value entry of one patch in an ASCII field file.

Parameter:
Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.sediment.read_log(path)[Quellcode]

Iterations done, the initial residuals of the last one and how the run ended.

Parameter:

path (Path)

Rückgabetyp:

dict[str, Any]

axqua.solvers.openfoam.sediment.shields(reynolds)[Quellcode]

Critical Shields parameter as the solver fits the curve to the particle Reynolds number: 0.6 below 0.2, 0.06 from 500 on, a polynomial in the logarithm between.

Parameter:

reynolds (numpy.ndarray)

Rückgabetyp:

numpy.ndarray

axqua.solvers.openfoam.sediment.bed_exchange(geometry, velocity, *, grain_size, fall, roughness, pick_up)[Quellcode]

Shear stress and pick-up at every bed face, by the formulas of the solver.

velocity is the velocity vector in the bed cells. Follows sediDriftFoam2.C line by line, including what is particular to it: the shear velocity from the logarithmic law at the height of the cell center, Dey’s factor on the critical shear stress of a sloping bed, van Rijn’s transport parameter, and the reference concentration limited to 0.1.

Returns the shear stress and its critical value [N/m2], the transport parameter T [-], the reference concentration [-] and the volume source of the bed cell S_explicit [1/s].

Parameter:
Rückgabetyp:

dict[str, numpy.ndarray]

class axqua.solvers.openfoam.sediment.BedResult(table, budget, layer=None, files=<factory>)[Quellcode]

Bases: object

The fine sediment at the bed, per bed face, and the budget of the reach.

Parameter:
table: dict[str, numpy.ndarray]
budget: dict[str, Any]
layer: Path | None = None
files: list[Path]
lines()[Quellcode]
Rückgabetyp:

list[str]

axqua.solvers.openfoam.sediment.evaluate(cfg, case_dir=None, *, write=True)[Quellcode]

Read the last result of the case and book the sediment at the bed.

Per bed face, in the plan of the column: the concentration in the cell at the bed, the settling and the pick-up per bed area and their difference; the shear stress with its critical value, as the solver computed them; and, where the case names a period, the mass the bed takes up or gives off in it and the share of fines of the bed layer at its end. The budget compares what enters and leaves the reach with what the bed exchanges.

Parameter:
Rückgabetyp:

BedResult

axqua.solvers.openfoam.sediment.BUDGET_TOLERANCE = 0.02

Share of the moving sediment the budget may leave open before the result says so.

axqua.solvers.openfoam.sediment.bed_layer(net, *, period, thickness, porosity, initial)[Quellcode]

What a net rate does to a bed layer whose surface stays where it is.

net is the mass the bed takes up per area and time [kg/m2/s], negative where it gives off. The layer consists of coarse grains that stay, grain density x thickness x (1 - porosity) per area, and of fines in the share initial of its dry mass. All fines the bed takes up are counted to the layer.

Two limits are marked in limit and only the first changes a number: the bed gives off no more fines than the layer holds (-1), and where the fines of the layer no longer fit into the pores of the coarse grains (+1) the surplus lies on the bed. The pores take fines that pack as the coarse grains do.

Parameter:
Rückgabetyp:

dict[str, numpy.ndarray]

axqua.solvers.openfoam.sediment.write_bed_layer(path, geometry, table, epsg)[Quellcode]

The bed as polygons in the plan, one per column, with the table as attributes.

Parameter:
Rückgabetyp:

Path

axqua.solvers.openfoam.sediment.set_iterations(case_dir, iterations)[Quellcode]

Change the largest number of iterations of a built case.

Parameter:
Rückgabetyp:

None

axqua.solvers.openfoam.sediment.run_case(cfg, *, case_dir=None, runtime=None, iterations=None, check_mesh=True, sink=None, show_progress=True, should_stop=None)[Quellcode]

Run sediDriftFoam2 on the built case, from its first iteration.

One process: Olsen’s solvers walk through the points of the mesh in the order of their numbers, which a mesh divided among processes does not keep, and the launchers of the installation refuse a parallel run.

Returns the finished process. A run that ends with an error raises SolverError with what the solver said last.

Parameter:

  • case_dir (str | Path | None)

  • iterations (int | None)

  • check_mesh (bool)

  • show_progress (bool)

Backend Adapter

The OpenFOAM backend: build, run, postprocess, export.

Like its TELEMAC counterpart this is an adapter, not new physics - build_case, the two-stage run and the discharge report already exist and are what openfoam_preprocessing.py / openfoam_run.py call.

One thing is deliberately absent: the TELEMAC pre-run. An OpenFOAM case is seeded from a converged 2D result, and prerun.ensure_seed produces one - but calling it from here would make the OpenFOAM backend reach into the TELEMAC backend, which tests/test_capabilities.py forbids by reading the source, and which would be wrong even if it were not checked: the two backends are siblings. Orchestration that drives both belongs above both, so the executor obtains the seed and hands it in - exactly what cli._run_openfoam already does.

class axqua.solvers.openfoam.backend.OpenFoamBackend[Quellcode]

Bases: BaseBackend

interFoam two-phase free-surface runs (Foundation OpenFOAM, openfoam.org).

name: str = 'openfoam'
check_environment(cfg)[Quellcode]
Rückgabetyp:

bool

build(cfg, capability=None, *, ctx=None, hotstart=None, cell_size=None, layers=None, seed=None, **_ignored)[Quellcode]

Mesh the reach and write the case.

seed is the 2D result to hotstart from, supplied by the caller. When it is absent the case still builds - cold, under a flat lid - which is the documented --no-pre-run behavior rather than a failure.

Parameter:
run(cfg, capability=None, *, ctx=None, n_processors=None, reconstruct=True, check_mesh=True, iterations=None, **_ignored)[Quellcode]
Parameter:
postprocess(cfg, ctx=None)[Quellcode]

Reports, and the optional VTK conversion.

Both are guarded: the run itself has already succeeded by the time this is called, and a missing figure is not worth discarding hours of solver time.

Parameter:

ctx (Any)

study(cfg, capability=None, *, ctx=None, launch_mode='prepare', checkout=None, rebuild=False, seed=None, **kwargs)[Quellcode]

Multi-run investigations: the calibration campaign and the mesh convergence study.

The 2D seed is passed in rather than obtained here: axqua.prerun imports the TELEMAC backend transitively, and this module must not reach into it (tests/test_capabilities.py enforces that by reading the source).

Parameter:
extract_objective(cfg, ctx=None)[Quellcode]

The relative inlet/outlet discharge imbalance - what the run is judged on.

Parameter:

ctx (Any)

Rückgabetyp:

float | None

describe_results(cfg)[Quellcode]

What a visualizer can open for this case.

Cheap by contract: directory listings and the 0/ field names, no pyvista and no solver launch, because a caller may ask this just to decide which figures are available. case.foam is handed over as-is - VisIt and ParaView both read an OpenFOAM case directory natively, so nothing is converted here.

Rückgabetyp:

list

export_qgis_results(cfg, ctx=None)[Quellcode]

List what QGIS can open. No PyQGIS, no copies of the case.

Parameter:

ctx (Any)