Cette page a été traduite automatiquement de l’anglais et peut contenir des erreurs. Original en anglais
Code source de axqua.solvers.openfoam.mesh
"""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 :mod:`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.
"""
from __future__ import annotations
import logging
from dataclasses import dataclass, field
from pathlib import Path
import numpy as np
from axqua.config import Config
from axqua.core import geodata
from axqua.core.geodata import dataset
from axqua.core.raster import sample_raster_at
from axqua.core.structures import (
SOLID, apply_to_bed, load_structures, solid_footprint,
)
from axqua.solvers.openfoam.polymesh import PolyMesh, assemble, validate
log = logging.getLogger("axqua")
BED_PATCH = "bed"
ATMOSPHERE_PATCH = "atmosphere"
LID_PATCH = "lid" # rigid-lid mode: the free surface, as a slip wall
BANKS_PATCH = "banks"
#: 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.
STRUCTURES_PATCH = "structures"
INLET_PREFIX = "inlet"
OUTLET_PREFIX = "outlet"
# --------------------------------------------------------------------------- #
# plan grid
# --------------------------------------------------------------------------- #
[docs]
@dataclass
class PlanGrid:
"""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.
"""
dx: float
angle: float # grid rotation [rad], CCW from east
origin: np.ndarray # (2,) world xy of the rotated lattice's node (0, 0)
nx: int
ny: int
col_id: np.ndarray # (ny, nx) int, -1 where the column is absent
vert_id: np.ndarray # (ny+1, nx+1) int, -1 where the vertex is unused
cell_xy: np.ndarray # (n_columns, 2) world coords of the column centers
vert_xy: np.ndarray # (n_vertices, 2) world coords of the used vertices
@property
def n_columns(self) -> int:
return int(self.cell_xy.shape[0])
@property
def n_vertices(self) -> int:
return int(self.vert_xy.shape[0])
@property
def keep(self) -> np.ndarray:
return self.col_id >= 0
@property
def cell_area(self) -> float:
return float(self.dx * self.dx)
def _rotation(angle: float) -> np.ndarray:
c, s = np.cos(angle), np.sin(angle)
return np.array([[c, -s], [s, c]])
[docs]
def flow_angle(cfg: Config) -> float:
"""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.
"""
if cfg.geodata.channel_centerline is None:
return 0.0
gdf = dataset(cfg).centerline_frame()
pts = []
for geom in gdf.geometry.values:
parts = geom.geoms if geom.geom_type == "MultiLineString" else [geom]
for part in parts:
pts.append(np.asarray(part.coords)[:, :2])
if not pts:
return 0.0
xy = np.concatenate(pts)
xy = xy - xy.mean(axis=0)
# principal axis = leading eigenvector of the coordinate covariance
_, _, vt = np.linalg.svd(xy, full_matrices=False)
return float(np.arctan2(vt[0, 1], vt[0, 0]))
[docs]
def build_plan_grid(polygon, dx: float, *, angle: float = 0.0,
max_columns: int = 4_000_000, blocked=None,
blocked_by_structures: bool = True) -> PlanGrid:
"""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.
"""
import shapely
from scipy import ndimage
rot = _rotation(angle)
coords = np.asarray(polygon.exterior.coords)[:, :2] if hasattr(polygon, "exterior") \
else np.asarray(polygon.convex_hull.exterior.coords)[:, :2]
centre = coords.mean(axis=0)
local = (coords - centre) @ rot # rotate into grid frame
lo = local.min(axis=0) - dx
hi = local.max(axis=0) + dx
nx = int(np.ceil((hi[0] - lo[0]) / dx))
ny = int(np.ceil((hi[1] - lo[1]) / dx))
if nx * ny > max_columns:
raise ValueError(
f"a {dx:g} m lattice over this domain needs {nx * ny:,} plan columns "
f"(cap {max_columns:,}). Coarsen openfoam.cell_size or shrink the domain "
f"(openfoam.domain: wetted)."
)
origin = centre + lo @ rot.T
i = np.arange(nx)
j = np.arange(ny)
uu, vv = np.meshgrid((i + 0.5) * dx, (j + 0.5) * dx)
world = np.stack([uu.ravel(), vv.ravel()], axis=1) @ rot.T + origin
keep = shapely.contains_xy(polygon, world[:, 0], world[:, 1]).reshape(ny, nx)
keep = ndimage.binary_fill_holes(keep)
if blocked is not None:
# AFTER the hole fill, never before: a building punched out of the domain is
# exactly the interior void that binary_fill_holes exists to close, so
# blanking first would simply be undone.
# A column is blanked when the solid INTERSECTS ITS CELL SQUARE, not when it
# contains the cell center.
#
# Testing the center lets a wall thinner than dx - or one that merely passes
# between two rows of centers - block nothing at all, and the mesh then joins
# the two sides of it. A structure is there precisely because its two sides
# are at different levels, so such a hole carries the whole head difference:
# on munich-vsf (3 cm lattice, sheet-steel baffles) 8,222 columns straddled a
# baffle and the fish pass short-circuited itself.
#
# The first fix here dilated the solid by half a cell diagonal and tested the
# center against that. It seals, but it also fattens every structure whether
# or not it was at risk, and on this geometry that was not a rounding error:
# it cut the narrowest opening across the fish pass from 2.25 m to 0.09 m,
# three cells, and the run built on it pressurised its upstream half to 20 m
# of head because the water had nowhere to go. A 5 m weir crest on a 1 m
# lattice would likewise have been widened 28%, silently changing a physical
# parameter.
#
# Intersection gives the guarantee without the fattening: a solid cannot cross
# a cell without intersecting it, so nothing can leak, while a structure that
# already spans several cells grows by strictly less than a half-diagonal. The
# cheap dilation is kept only as a CANDIDATE FILTER - exact intersection is
# then evaluated on those few thousand cells rather than on all of them.
half = float(dx) / 2.0
near = shapely.contains_xy(blocked.buffer(float(dx) * np.sqrt(2.0) / 2.0),
world[:, 0], world[:, 1])
cand = np.flatnonzero(near)
flat = np.zeros(world.shape[0], dtype=bool)
if cand.size:
# Axis-aligned in the GRID frame, then rotated, so the test stays exact
# when align_to_flow has turned the lattice off the world axes.
corners = np.array([[-half, -half], [half, -half],
[half, half], [-half, half]]) @ rot.T
squares = shapely.polygons(
world[cand][:, None, :] + corners[None, :, :])
flat[cand] = shapely.intersects(blocked, squares)
solid = flat.reshape(ny, nx)
leaky = int((keep & solid).sum()) - int(
(keep & shapely.contains_xy(blocked, world[:, 0],
world[:, 1]).reshape(ny, nx)).sum())
removed = int((keep & solid).sum())
keep = keep & ~solid
if removed:
# Same count, honest attribution: under a rigid lid `blocked` carries the
# dry part of the reach as well as the structures, and calling a trimmed
# gravel bar a "solid structure" is how a reader goes looking for a wall
# that was never drawn.
what = ("solid structures" if blocked_by_structures
else "the blanking footprint")
log.info("plan grid: %d columns removed by %s (%.0f m2), "
"of which %d come from sealing them to the lattice (a column is "
"blanked where the geometry crosses its %.3f m cell, so none of "
"them leak)",
removed, what, removed * dx * dx, leaky, float(dx))
labels, n = ndimage.label(keep) # 4-connectivity (the default structure)
if n > 1:
sizes = ndimage.sum(keep, labels, range(1, n + 1))
largest = labels == (int(np.argmax(sizes)) + 1)
dropped = int(keep.sum() - largest.sum())
keep = largest
by_structure = blocked_by_structures and dropped
level = log.warning if by_structure else log.info
level("plan grid: kept the largest of %d connected blocks (%d columns, "
"%d dropped)%s", n, int(keep.sum()), dropped,
" - a solid structure has cut the domain in two; check that a wall was "
"not drawn right across the channel" if by_structure else "")
if not keep.any():
raise ValueError("the plan lattice has no column inside the domain polygon; "
"check openfoam.cell_size against the domain size")
col_id = np.full((ny, nx), -1, dtype=np.int64)
col_id[keep] = np.arange(int(keep.sum()))
used = np.zeros((ny + 1, nx + 1), dtype=bool)
for dj in (0, 1):
for di in (0, 1):
used[dj:dj + ny, di:di + nx] |= keep
vert_id = np.full((ny + 1, nx + 1), -1, dtype=np.int64)
vert_id[used] = np.arange(int(used.sum()))
cj, ci = np.nonzero(keep)
cell_xy = np.stack([(ci + 0.5) * dx, (cj + 0.5) * dx], axis=1) @ rot.T + origin
vj, vi = np.nonzero(used)
vert_xy = np.stack([vi * dx, vj * dx], axis=1) @ rot.T + origin
log.info("plan grid: %d x %d lattice at %.3f m, %d columns kept (%.0f m2)",
nx, ny, dx, int(keep.sum()), keep.sum() * dx * dx)
return PlanGrid(dx=dx, angle=angle, origin=origin, nx=nx, ny=ny,
col_id=col_id, vert_id=vert_id,
cell_xy=cell_xy, vert_xy=vert_xy)
# --------------------------------------------------------------------------- #
# vertical distribution
# --------------------------------------------------------------------------- #
[docs]
def sigma_levels(height: np.ndarray, n_layers: int, *,
bed_layer: float | None = None,
expansion: float = 1.0) -> np.ndarray:
"""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.
"""
height = np.asarray(height, dtype=float)
nv = height.size
if n_layers < 1:
raise ValueError("n_layers must be >= 1")
if expansion <= 0:
raise ValueError("expansion must be > 0")
if np.isclose(expansion, 1.0):
base = np.ones(n_layers)
else:
q = expansion ** (1.0 / max(n_layers - 1, 1))
base = q ** np.arange(n_layers)
base = base / base.sum() # (n_layers,)
s = np.tile(np.concatenate([[0.0], np.cumsum(base)]), (nv, 1))
if bed_layer is not None and n_layers > 1:
with np.errstate(divide="ignore", invalid="ignore"):
frac = np.where(height > 0, np.minimum(bed_layer, 0.5 * height) / height, 0.0)
# keep the requested distribution above the pinned bed layer
upper = s[:, 1:] - s[:, 1:2]
upper = upper / np.where(upper[:, -1:] > 0, upper[:, -1:], 1.0)
s = np.concatenate([np.zeros((nv, 1)),
frac[:, None] + upper * (1.0 - frac)[:, None]], axis=1)
s[:, -1] = 1.0
return s
# --------------------------------------------------------------------------- #
# the extruded mesh
# --------------------------------------------------------------------------- #
[docs]
@dataclass
class OpenFoamMesh:
"""The built 3D mesh plus the bookkeeping the field writers need."""
polymesh: PolyMesh
grid: PlanGrid
n_layers: int
z: np.ndarray # (NV, n_layers+1) vertex levels
bed: np.ndarray # (NV,) bed elevation at the plan vertices
lid: np.ndarray # (NV,) lid elevation at the plan vertices
cell_centres: np.ndarray # (NC, 3)
cell_column: np.ndarray # (NC,) column index of each cell
cell_layer: np.ndarray # (NC,) 0-based layer index of each cell
column_bed: np.ndarray # (n_columns,) bed at the column centers
column_wse: np.ndarray | None = None # (n_columns,) hotstart free surface
column_depth: np.ndarray | None = None # (n_columns,) hotstart water depth
column_uv: np.ndarray | None = None # (n_columns, 2) hotstart depth-averaged U
bed_ks: np.ndarray | None = None # (n bed faces,) Nikuradse ks per bed face
rigid_lid: bool = False # water-only domain under a slip lid
#: :func:`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_median: float | None = None
lid_step_p99: float | None = None
inlet_patches: list[str] = field(default_factory=list)
outlet_patches: list[str] = field(default_factory=list)
inlet_discharge: dict[str, float] = field(default_factory=dict)
notes: list[str] = field(default_factory=list)
@property
def n_cells(self) -> int:
return self.polymesh.n_cells
@property
def wall_patches(self) -> list[str]:
"""Every no-slip wall patch: 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.
"""
top = self.top_patch
return [p.name for p in self.polymesh.patches
if p.type == "wall" and p.name != top]
@property
def top_patch(self) -> str:
"""``lid`` under a rigid lid, ``atmosphere`` for the two-phase case."""
return LID_PATCH if self.rigid_lid else ATMOSPHERE_PATCH
@property
def first_layer_height(self) -> np.ndarray:
"""Thickness of the bed-adjacent layer at every plan vertex [m]."""
return (self.z[:, 1] - self.z[:, 0])
@property
def min_layer_height(self) -> 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.
"""
return float(np.diff(self.z, axis=1).min())
@property
def cell_bounds(self) -> tuple[np.ndarray, np.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.
"""
corners = column_corner_vertices(self.grid)
z_corner = self.z[corners] # (n_columns, 4, n_layers+1)
level = z_corner.mean(axis=1) # (n_columns, n_layers+1)
return level[:, :-1].ravel(), level[:, 1:].ravel()
[docs]
def column_corner_vertices(grid: PlanGrid) -> np.ndarray:
"""(n_columns, 4) vertex ids of each column's plan corners, CCW in grid frame."""
cj, ci = np.nonzero(grid.keep)
return np.stack([grid.vert_id[cj, ci],
grid.vert_id[cj, ci + 1],
grid.vert_id[cj + 1, ci + 1],
grid.vert_id[cj + 1, ci]], axis=1)
def _levels_to_points(grid: PlanGrid, z: np.ndarray) -> np.ndarray:
"""Flatten the per-vertex level table into the OpenFOAM point list.
Point ordering is ``vertex * (n_layers + 1) + k``; every downstream index
derives from that one rule.
"""
nk = z.shape[1]
xy = np.repeat(grid.vert_xy, nk, axis=0)
return np.column_stack([xy, z.ravel()])
[docs]
def extrude(grid: PlanGrid, bed: np.ndarray, lid: np.ndarray, n_layers: int, *,
bed_layer: float | None = None, expansion: float = 1.0,
min_height: float = 0.05) -> tuple[PolyMesh, np.ndarray, np.ndarray,
np.ndarray, np.ndarray, dict]:
"""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.
"""
nk = n_layers + 1
height = np.maximum(lid - bed, min_height)
s = sigma_levels(height, n_layers, bed_layer=bed_layer, expansion=expansion)
z = bed[:, None] + height[:, None] * s
points = _levels_to_points(grid, z)
corners = column_corner_vertices(grid) # (nc, 4)
nc = grid.n_columns
def pid(vert: np.ndarray, k) -> np.ndarray:
return vert * nk + k
# ---- cell centers (mean of the eight corner points) ----------------------
kk = np.arange(n_layers)
zc = 0.5 * (z[corners][:, :, kk] + z[corners][:, :, kk + 1]) # (nc, 4, n_layers)
cell_centres = np.empty((nc * n_layers, 3))
cell_centres[:, 0] = np.repeat(grid.cell_xy[:, 0], n_layers)
cell_centres[:, 1] = np.repeat(grid.cell_xy[:, 1], n_layers)
cell_centres[:, 2] = zc.mean(axis=1).ravel()
cell_column = np.repeat(np.arange(nc), n_layers)
cell_layer = np.tile(kk, nc)
keep = grid.keep
col = grid.col_id
ny, nx = keep.shape
lay = np.arange(n_layers)
def _stack(v_lo: np.ndarray, v_hi: np.ndarray) -> np.ndarray:
"""Vertical quads spanning plan vertices *v_lo*..*v_hi* over every layer."""
vlo = np.repeat(v_lo, n_layers)
vhi = np.repeat(v_hi, n_layers)
k = np.tile(lay, v_lo.size)
return np.stack([pid(vlo, k), pid(vhi, k),
pid(vhi, k + 1), pid(vlo, k + 1)], axis=1)
def _cells(cols: np.ndarray) -> np.ndarray:
return np.repeat(cols, n_layers) * n_layers + np.tile(lay, cols.size)
internal_quads, internal_own, internal_nei = [], [], []
# faces between columns adjacent in x (grid frame)
mx = keep[:, :-1] & keep[:, 1:]
if mx.any():
j, i = np.nonzero(mx)
internal_quads.append(_stack(grid.vert_id[j, i + 1], grid.vert_id[j + 1, i + 1]))
internal_own.append(_cells(col[j, i]))
internal_nei.append(_cells(col[j, i + 1]))
# faces between columns adjacent in y
my = keep[:-1, :] & keep[1:, :]
if my.any():
j, i = np.nonzero(my)
internal_quads.append(_stack(grid.vert_id[j + 1, i], grid.vert_id[j + 1, i + 1]))
internal_own.append(_cells(col[j, i]))
internal_nei.append(_cells(col[j + 1, i]))
# horizontal faces inside a column
if n_layers > 1:
klev = np.tile(np.arange(1, n_layers), nc)
cc = np.repeat(corners, n_layers - 1, axis=0)
internal_quads.append(np.stack([pid(cc[:, 0], klev), pid(cc[:, 1], klev),
pid(cc[:, 2], klev), pid(cc[:, 3], klev)], axis=1))
base = np.repeat(np.arange(nc), n_layers - 1) * n_layers
internal_own.append(base + klev - 1)
internal_nei.append(base + klev)
internal = (np.concatenate(internal_own), np.concatenate(internal_nei),
np.concatenate(internal_quads, axis=0))
# ---- boundary faces ------------------------------------------------------
zero = np.zeros(nc, dtype=np.int64)
bed_quads = np.stack([pid(corners[:, 0], zero), pid(corners[:, 1], zero),
pid(corners[:, 2], zero), pid(corners[:, 3], zero)], axis=1)
bed_owner = np.arange(nc) * n_layers
top = np.full(nc, n_layers, dtype=np.int64)
top_quads = np.stack([pid(corners[:, 0], top), pid(corners[:, 1], top),
pid(corners[:, 2], top), pid(corners[:, 3], top)], axis=1)
top_owner = np.arange(nc) * n_layers + (n_layers - 1)
padded = np.zeros((ny + 2, nx + 2), dtype=bool)
padded[1:-1, 1:-1] = keep
side_quads, side_owner, side_mid = [], [], []
# (neighbor offset, the two plan vertices of the exposed face)
exposures = [
((0, -1), lambda j, i: (grid.vert_id[j, i], grid.vert_id[j + 1, i])),
((0, +1), lambda j, i: (grid.vert_id[j, i + 1], grid.vert_id[j + 1, i + 1])),
((-1, 0), lambda j, i: (grid.vert_id[j, i], grid.vert_id[j, i + 1])),
((+1, 0), lambda j, i: (grid.vert_id[j + 1, i], grid.vert_id[j + 1, i + 1])),
]
for (dj, di), verts in exposures:
exposed = keep & ~padded[1 + dj:1 + dj + ny, 1 + di:1 + di + nx]
if not exposed.any():
continue
j, i = np.nonzero(exposed)
v_lo, v_hi = verts(j, i)
side_quads.append(_stack(v_lo, v_hi))
side_owner.append(_cells(col[j, i]))
side_mid.append(0.5 * (grid.vert_xy[v_lo] + grid.vert_xy[v_hi]))
sides = {
"quads": np.concatenate(side_quads, axis=0) if side_quads
else np.zeros((0, 4), dtype=np.int64),
"owner": np.concatenate(side_owner) if side_owner
else np.zeros(0, dtype=np.int64),
# one midpoint per exposed *column edge*; each spawns n_layers faces
"midpoint": np.repeat(np.concatenate(side_mid), n_layers, axis=0)
if side_mid else np.zeros((0, 2)),
}
return (points, internal, (bed_quads, bed_owner), (top_quads, top_owner), sides,
z, cell_centres, cell_column, cell_layer)
# --------------------------------------------------------------------------- #
# patch classification
# --------------------------------------------------------------------------- #
def _tailwater_for(cfg, patch: str, tailwater: float | None = None):
"""The level this outlet holds, or None when the case prescribes none.
``openfoam.outlet_stage`` first, then *tailwater*, the level the build resolved
from the outflow condition of the case (under ``stage_discharge`` the rating
curve at the discharge of the build, which ``prescribed_elevation`` does not
hold), then ``boundaries.prescribed_elevation``.
"""
stage = getattr(cfg.openfoam, "outlet_stage", None)
if isinstance(stage, dict):
stage = stage.get(patch)
if stage is None:
stage = tailwater
if stage is None:
stage = getattr(cfg.boundaries, "prescribed_elevation", None)
return None if stage is None else float(stage)
[docs]
def split_outlet_above_tailwater(owner, cell_centres, stage, margin, *,
cell_size=0.0):
"""``(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.
"""
z = np.asarray(cell_centres, dtype=float)[np.asarray(owner), 2]
ceiling = float(stage) + float(margin)
walled = z > ceiling
if walled.all():
# every face above the line: the stage is below the whole patch, which is a
# configuration error rather than something to silently act on
log.warning(
"outlet: the tailwater %.3f m + margin %.3f m is below every face of the "
"patch (lowest %.3f m); leaving it fully open rather than sealing the "
"only way out of the domain", stage, margin, z.min())
return np.ones(z.shape, dtype=bool), np.zeros(z.shape, dtype=bool)
open_span = float(np.ptp(z[~walled])) if (~walled).any() else 0.0
if cell_size and (~walled).any() and open_span < 2 * cell_size:
log.warning(
"outlet: only %.3f m of it stays open (tailwater %.3f + margin %.3f), "
"under two cells. Raise openfoam.outlet_air_margin, or the outflow is "
"squeezed through a slot.", open_span, stage, margin)
return ~walled, walled
[docs]
def split_structure_sides(midpoints: np.ndarray, banks: np.ndarray,
blocked, dx: float) -> tuple[np.ndarray, np.ndarray]:
"""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.
"""
import shapely
on_structure = np.zeros(banks.shape, dtype=bool)
if blocked is not None and not blocked.is_empty and banks.any():
mid = midpoints[banks]
near = blocked.buffer(float(dx) * np.sqrt(2.0) / 2.0)
on_structure[banks] = shapely.contains_xy(near, mid[:, 0], mid[:, 1])
return banks & ~on_structure, on_structure
[docs]
def classify_sides(cfg: Config, midpoints: np.ndarray,
tolerance: float) -> tuple[np.ndarray, list[str], dict[str, float]]:
"""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.
"""
import shapely
from axqua.core.boundaries import liquid_line_details, liquid_lines
labels = np.full(midpoints.shape[0], BANKS_PATCH, dtype=object)
names: list[str] = []
discharges: dict[str, float] = {}
cfg = _with_submodel_boundaries(cfg)
if cfg.boundaries.liquid_boundaries is None or midpoints.shape[0] == 0:
return labels, names, discharges
details = liquid_line_details(cfg)
if not details:
merged = liquid_lines(cfg)
details = [{"kind": k, "discharge": None, "geom": g} for k, g in merged.items()]
pts = shapely.points(midpoints[:, 0], midpoints[:, 1])
best = np.full(midpoints.shape[0], np.inf)
counts: dict[str, int] = {}
n_in = n_out = 0
for entry in details:
kind = entry["kind"]
if kind == "inflow":
n_in += 1
name = f"{INLET_PREFIX}-{n_in}"
else:
n_out += 1
name = f"{OUTLET_PREFIX}-{n_out}"
dist = shapely.distance(pts, entry["geom"])
hit = (dist <= tolerance) & (dist < best)
if hit.any():
labels[hit] = name
best[hit] = dist[hit]
names.append(name)
counts[name] = int(hit.sum())
if entry.get("discharge") is not None:
discharges[name] = float(entry["discharge"])
for name, n in counts.items():
log.info("patch %s: %d lateral faces within %.2f m of its liquid line",
name, n, tolerance)
return labels, names, discharges
# --------------------------------------------------------------------------- #
# inputs sampled onto the grid
# --------------------------------------------------------------------------- #
def _roughness_at(cfg: Config, xy: np.ndarray) -> np.ndarray | None:
"""Nikuradse ``ks`` [m] at arbitrary plan points, from the roughness zones.
The same ``geodata.roughness_zones`` + ``roughness_table`` pair the TELEMAC build
turns into ``FRIC_ID``/``.tbl`` rows, so the 2D and 3D models are rough in the
same places; ``None`` when the case has no zonation (the caller then falls back
to a single ``friction`` value).
"""
if cfg.geodata.roughness_zones is None or cfg.geodata.roughness_table is None:
return None
import shapely
from scipy.spatial import cKDTree
table = dataset(cfg).roughness_table()
zones = dataset(cfg).roughness_zones()
field_name = cfg.mesh.roughness_zone_field
if field_name not in zones.columns:
log.warning("roughness zones have no %r column; skipping per-face ks",
field_name)
return None
ids = zones[field_name].astype(int).to_numpy()
ks = np.array([table.get(int(i), np.nan) for i in ids], dtype=float)
out = np.full(xy.shape[0], np.nan)
pts = shapely.points(xy[:, 0], xy[:, 1])
tree = shapely.STRtree(zones.geometry.values)
hits = tree.query(pts, predicate="within")
if hits.size:
out[hits[0]] = ks[hits[1]]
missing = np.isnan(out)
if missing.any(): # outside every polygon -> nearest zone centroid
cent = np.array([[g.x, g.y] for g in zones.geometry.centroid.values])
_, idx = cKDTree(cent).query(xy[missing])
out[missing] = ks[idx]
return out
def _smooth_lid(grid: PlanGrid, values: np.ndarray, *, passes: int,
dilate: int) -> np.ndarray:
"""Gray-dilate then average the lid over the plan lattice.
A lid taken verbatim from a 2D result is as ragged as that result: local dips
over shallow cells would pinch the air layer to nothing and make thin, badly
shaped cells. Dilating first (a running maximum) guarantees the smoothed lid is
never *below* the surface it came from, then averaging removes the steps.
"""
from scipy import ndimage
used = grid.vert_id >= 0
field2d = np.full(grid.vert_id.shape, np.nan)
field2d[used] = values
if dilate > 0:
filled = np.where(np.isnan(field2d), -np.inf, field2d)
field2d = np.where(used, ndimage.grey_dilation(filled, size=2 * dilate + 1),
np.nan)
for _ in range(max(passes, 0)):
filled = np.where(np.isnan(field2d), 0.0, field2d)
weight = ndimage.uniform_filter(used.astype(float), size=3, mode="nearest")
avg = ndimage.uniform_filter(filled, size=3, mode="nearest")
with np.errstate(invalid="ignore", divide="ignore"):
field2d = np.where(used & (weight > 0), avg / weight, field2d)
return field2d[used]
[docs]
def resolve_headroom(cfg: Config, state, *, notes: list[str] | None = None) -> float:
"""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 (:meth:`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.
"""
of = cfg.openfoam
if state is None or of.headroom_mode != "auto" or of.mode == "rigid-lid":
return float(of.freeboard)
derived = state.headroom(of.freeboard)
if notes is not None and derived > of.freeboard + 1e-9:
notes.append(f"freeboard raised {of.freeboard:g} -> {derived:.2f} m: at "
f"{state.velocity_scale():.2f} m/s the flow can convert "
f"{state.velocity_scale() ** 2 / 19.62:.2f} m of velocity head, "
"and the surface must be able to rise into it")
return float(derived)
[docs]
def resolve_margin(cfg: Config, state) -> float:
"""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 :meth:`State2D.lateral_margin`).
"""
of = cfg.openfoam
if state is None or of.headroom_mode != "auto":
return float(of.wet_margin)
return float(state.lateral_margin(of.wet_margin, bank_slope=of.bank_slope))
[docs]
def resolve_lid(cfg: Config, grid: PlanGrid, bed: np.ndarray,
wse_v: np.ndarray | None, freeboard: float, *,
notes: list[str] | None = None) -> np.ndarray:
"""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.
"""
of = cfg.openfoam
rigid = of.mode == "rigid-lid"
notes = notes if notes is not None else []
if wse_v is not None:
lid = wse_v + freeboard
else:
lid = np.full_like(bed, float(np.nanmax(bed)) + of.freeboard)
notes.append("no 2D hotstart: flat lid, so the domain carries the full air column")
if rigid:
# The lid IS the free surface. No freeboard, because a freeboard would be air.
lid = np.where(np.isfinite(wse_v), wse_v, bed)
notes.append("lid placed on the 2D free surface (freeboard ignored)")
if of.lid == "flat" and not rigid:
lid = np.full_like(lid, float(np.nanmax(lid)))
notes.append(f"flat lid at {lid[0]:.3f} m a.s.l.")
if of.lid_elevation is not None:
lid = np.full_like(lid, float(of.lid_elevation))
notes.append(f"lid pinned to the configured {of.lid_elevation:g} m a.s.l.")
floor = of.min_water_depth if rigid else of.min_column_height
lid = np.maximum(lid, bed + floor)
if of.lid == "follow" and not rigid:
lid = _smooth_lid(grid, lid, passes=of.lid_smoothing, dilate=of.lid_dilation)
lid = np.maximum(lid, bed + of.min_column_height)
# Under a rigid lid the lid IS the converged 2D free surface: already smooth, and
# smoothing it again only pushes it into the bed in steep places, where the
# re-clamp then leaves a kink that inverts cells.
return lid
[docs]
def resolve_layers(cfg: Config, bed_ks: np.ndarray | None,
column_height: np.ndarray | None = None, *,
notes: list[str] | None = None) -> int:
"""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.
"""
of = cfg.openfoam
notes = notes if notes is not None else []
n_layers = int(of.n_layers)
if of.mode != "rigid-lid":
return n_layers
thinnest = float(np.nanmedian(bed_ks)) if bed_ks is not None else 0.05
thinnest = max(thinnest, 0.05)
fit = max(2, int(of.min_water_depth / thinnest))
if fit >= n_layers:
return n_layers
if of.auto_layers is not False:
notes.append(
f"{n_layers} layers reduced to {fit}: the shallowest meshed column is "
f"{of.min_water_depth:g} m, and {n_layers} layers would make cells "
f"thinner than the bed roughness ({thinnest:.3f} m)")
return fit
typical = (float(np.nanmedian(column_height)) if column_height is not None
else float(of.min_water_depth))
notes.append(
f"auto_layers: false - keeping {n_layers} layers where the roughness would "
f"allow {fit}. Over the median {typical:.2f} m column that is "
f"{typical / n_layers:.3f} m cells, but in the shallowest meshed column "
f"({of.min_water_depth:g} m) it is {of.min_water_depth / n_layers:.3f} m "
f"against a {thinnest:.3f} m bed roughness - check the quality report for "
"folded cells there")
return n_layers
def _domain_polygon(cfg: Config, state, *, domain: str, wet_margin: float,
wet_depth: float):
"""The plan footprint the lattice covers: the ROI, or the wetted corridor."""
from shapely.ops import unary_union
roi, roi_label = _openfoam_roi(cfg)
if domain == "roi" or state is None:
return roi, roi_label
wet = state.wet_footprint(wet_depth=wet_depth, buffer=wet_margin)
if wet is None or wet.is_empty:
log.warning("no wetted footprint in the 2D result; meshing the full ROI")
return roi, "the full ROI boundary (no wetted footprint found)"
clipped = unary_union([wet]).intersection(roi)
if clipped.is_empty:
return roi, "the full ROI boundary (wetted footprint missed the ROI)"
if clipped.geom_type == "MultiPolygon":
clipped = max(clipped.geoms, key=lambda g: g.area)
return clipped, (f"the 2D wetted extent (H > {wet_depth:g} m) buffered by "
f"{wet_margin:g} m, clipped to {roi_label}")
def _with_submodel_boundaries(cfg: Config) -> Config:
"""*cfg*, with the sub-model's own liquid lines in place of the case's.
A cropped OpenFOAM domain no longer touches the boundaries the reach enters and
leaves by, so it carries its own - swapped in here rather than in the case config,
which the 2D run is converged against.
"""
override = getattr(cfg.openfoam, "liquid_boundaries", None)
if not override:
return cfg
import copy
local = copy.copy(cfg)
local.boundaries = copy.copy(cfg.boundaries)
local.boundaries.liquid_boundaries = Path(override)
return local
def _openfoam_roi(cfg: Config):
"""(polygon, label) for the ROI the OpenFOAM domain is cut from.
``openfoam.roi`` overrides the case boundary so a sub-model - a structure, a fish
pass - can be meshed in 3D without the uniform channel above and below it, while
the 2D case keeps the full reach it was built and converged on.
"""
override = getattr(cfg.openfoam, "roi", None)
if not override:
return dataset(cfg).roi_polygon(), "the full ROI boundary"
import geopandas as gpd
from shapely.ops import polygonize, unary_union
path = Path(override)
if not path.exists():
raise FileNotFoundError(
f"openfoam.roi points at {path}, which does not exist. It is written by "
"the case's own ROI helper; run that first, or unset openfoam.roi to mesh "
"the whole reach.")
gdf = gpd.read_file(path)
geoms = [g for g in gdf.geometry.values if g is not None and not g.is_empty]
if not geoms:
raise ValueError(f"openfoam.roi {path.name} holds no geometry")
if set(gdf.geom_type) & {"LineString", "MultiLineString"}:
polys = list(polygonize(unary_union(geoms)))
if not polys:
raise ValueError(f"{path.name}: boundary lines do not close into a polygon")
return max(polys, key=lambda p: p.area), f"the OpenFOAM ROI {path.name}"
return unary_union(geoms), f"the OpenFOAM ROI {path.name}"
# --------------------------------------------------------------------------- #
# the public builder
# --------------------------------------------------------------------------- #
[docs]
def plan_spacing(cfg: Config) -> float:
"""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.
"""
factor = cfg.openfoam.cell_size_factor
if not factor:
return float(cfg.openfoam.cell_size)
return float(factor) * float(geodata.nominal_channel_size(cfg))
[docs]
def lid_steps(grid: PlanGrid, lid: np.ndarray, bed: np.ndarray, *,
radius: int = 2) -> tuple[float, float]:
"""``(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.
"""
from scipy.ndimage import maximum_filter, minimum_filter
corners = column_corner_vertices(grid)
col_lid = lid[corners].mean(axis=1)
col_depth = np.maximum((lid - bed)[corners].mean(axis=1), 1e-3)
field = np.full(grid.col_id.shape, np.nan)
inside = grid.col_id >= 0
field[inside] = col_lid[grid.col_id[inside]]
size = 2 * radius + 1
span = (maximum_filter(np.nan_to_num(field, nan=-1e9), size)
- minimum_filter(np.nan_to_num(field, nan=+1e9), size))
per_col = np.full(grid.n_columns, np.nan)
per_col[grid.col_id[inside]] = span[inside]
ratio = per_col / col_depth
ratio = ratio[np.isfinite(ratio)]
if not ratio.size:
return 0.0, 0.0
return float(np.median(ratio)), float(np.percentile(ratio, 99))
[docs]
def build_mesh(cfg: Config, *, state=None, dem: str | Path | None = None,
tailwater: float | None = None) -> OpenFoamMesh:
"""Build the terrain-following hex mesh for *cfg*'s OpenFOAM case.
*state* is an optional :class:`~axqua.solvers.openfoam.hotstart.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.
"""
of = cfg.openfoam
rigid = of.mode == "rigid-lid"
if rigid and state is None:
# Said here rather than discovered at the lid, where it used to surface as a
# bare NameError: under a rigid lid the lid IS the 2D free surface, so without
# a 2D result there is no lid to place and nothing to fall back on. A flat one
# would be a horizontal water surface, which is the thing this mode exists to
# avoid.
raise ValueError(
"openfoam.mode: rigid-lid needs a 2D result to take the free surface from, "
"and none was given. Run the pre-run (openfoam.pre_run.enabled) or pass a "
"--hotstart, or switch to mode: vof to solve the surface instead.")
dx = plan_spacing(cfg)
dem = Path(dem) if dem is not None else Path(cfg.geodata.dem_initial)
notes: list[str] = []
lid_step_median = lid_step_p99 = None
if of.cell_size_factor:
notes.append(f"plan spacing {dx:.2f} m = {of.cell_size_factor:g}x the 2D "
f"channel size ({geodata.nominal_channel_size(cfg):.2f} m)")
if rigid:
notes.append(
"RIGID LID: the domain is water only, the lid sits on the 2D free surface "
"as a slip wall, and there is no air phase. The free surface is an input, "
"not a result - it cannot rise, overtop or wet a dry bar.")
polygon, how = _domain_polygon(cfg, state, domain=of.domain,
wet_margin=resolve_margin(cfg, state),
wet_depth=cfg.hydrodynamics.wet_depth)
notes.append(f"footprint: {how}")
structures = load_structures(cfg)
blocked = solid_footprint(structures)
dry_footprint = None
if rigid and state is not None:
# Water only: a dry column has no water column to mesh, and keeping it would
# add a sliver of cells doing nothing. The wetted edge becomes a vertical wall
# at the waterline - the standard rigid-lid domain, and an approximation worth
# naming: the shoreline cannot move.
# The trim depth is the LAYER thickness, not the reporting wet_depth: a
# column 1 cm deep cannot be meshed into n_layers cells without collapsing
# into slivers, and checkMesh rejects those as incorrectly oriented face
# pyramids. One cell of dry ground beyond that line keeps the liquid-boundary
# lines (which sit at the water's edge) inside the domain.
trim_depth = max(cfg.hydrodynamics.wet_depth, of.min_water_depth)
dry_footprint = state.dry_footprint(wet_depth=trim_depth, margin=dx)
notes.append(f"columns shallower than {trim_depth:g} m left out: a water "
"column thinner than one cell cannot be meshed")
if dry_footprint is not None and not dry_footprint.is_empty:
from shapely.ops import unary_union
blocked = (dry_footprint if blocked is None
else unary_union([blocked, dry_footprint]))
n_solid = sum(1 for s in structures if s.mode == SOLID)
if n_solid:
notes.append(f"{n_solid} solid structure(s) removed from the domain: their "
"sides become no-slip walls from bed to lid")
angle = flow_angle(cfg) if of.align_to_flow else 0.0
if of.align_to_flow:
notes.append(f"lattice rotated {np.degrees(angle):+.1f} deg onto the reach axis")
grid = build_plan_grid(polygon, dx, angle=angle,
max_columns=of.max_plan_columns, blocked=blocked,
blocked_by_structures=bool(n_solid))
if rigid:
notes.append("domain trimmed to the wetted body; the waterline is a vertical "
"wall and cannot move")
bed = sample_raster_at(dem, grid.vert_xy[:, 0], grid.vert_xy[:, 1])
column_bed = sample_raster_at(dem, grid.cell_xy[:, 0], grid.cell_xy[:, 1])
if structures:
# overflow structures ARE terrain the DEM did not have; raising the bed keeps
# the plan mesh unchanged and lets water pass over the crest
bed, touched = apply_to_bed(structures, grid.vert_xy, bed)
column_bed, _ = apply_to_bed(structures, grid.cell_xy, column_bed)
if touched.any():
notes.append(f"{int(touched.sum())} mesh vertices raised to an overflow "
"structure crest")
# ---- lid ---------------------------------------------------------------
column_wse = column_depth = column_uv = None
wse_v = None
if state is not None:
wse_v = state.sample_surface(grid.vert_xy)
column_wse, column_depth, column_uv = state.sample_columns(grid.cell_xy)
wse_v = np.where(np.isfinite(wse_v), wse_v, bed)
wse_v = np.maximum(wse_v, bed)
freeboard = resolve_headroom(cfg, state, notes=notes) if state is not None \
else float(of.freeboard)
lid = resolve_lid(cfg, grid, bed, wse_v, freeboard, notes=notes)
floor = of.min_water_depth if rigid else of.min_column_height
if rigid:
# Asked here, at build time, because the answer decides whether the mode is
# applicable at all - and the run itself will not tell you. munich-vsf
# finished, balanced its discharge to -0.000% and reported a healthy Courant
# number for 61 hours while its velocity cap held 3,719 cells together.
lid_step_median, lid_step_p99 = step_med, step_p99 = lid_steps(grid, lid, bed)
notes.append(f"lid steps: the prescribed surface moves {100 * step_med:.1f}% "
f"of the local depth within 2 cells (p99 {100 * step_p99:.0f}%)")
# from `report`, which owns the two grades, so the refusal here and the
# verdict a reader sees afterwards are the same numbers rather than copies
from axqua.solvers.openfoam.report import LID_STEP_SEVERE, LID_STEP_WARN
# LEG C IS EXEMPT, and the exemption is the whole reason it exists. The
# refusal below rests on the top boundary being a WALL: a wall cannot answer
# a drop by plunging, so it converts the head into velocity instead. Under
# potentialFreeSurfaceFoam the top is not a wall - the surface is free to
# move there, carried as `zeta` - so the premise does not hold and the step
# is a mesh-quality question (which `quality` reports) rather than a mode
# error. Keeping the lid ON the stepped 2D surface is also what keeps leg C's
# OWN validity condition satisfied: its linearisation needs |zeta - z_lid|
# small, and a lid that already follows the surface leaves only the 3D
# correction for zeta to carry.
from . import potential
if potential.is_potential(cfg):
notes.append(
"lid steps are not a refusal here: the top boundary is a free "
"surface (zeta), not a wall, so it answers a drop by moving")
elif step_p99 > LID_STEP_SEVERE and not of.allow_stepped_lid:
# A refusal, not a warning: above one depth the surface steps further than
# the water is deep, which is a weir, a drop, or the slots of a fish pass.
# A lid is a WALL there, and the velocity it produces is an artifact of
# the mode. Two munich-vsf runs finished, balanced their discharge to
# -0.000%, held Courant at 0.90 for 61 hours, and were invalid - so the
# cost of building anyway is days, and the cost of stopping is a rebuild.
raise ValueError(
f"RIGID LID DOES NOT APPLY HERE: at the 99th percentile the "
f"prescribed free surface steps {100 * step_p99:.0f}% of the local "
f"water depth within two cells (median {100 * step_med:.1f}%), which "
f"is more than the water is deep. A lid is a WALL, so it cannot "
f"answer a drop by plunging the way a free surface does - it converts "
f"the head into velocity instead (sqrt(2 g dz): a 1 m step drives "
f"4.4 m/s), and the fastest cells in the result will be that artifact "
f"rather than the flow.\n"
f" Use openfoam.mode: vof where the surface has to move; or check "
f"that the steps are real rather than the 2D result read across a "
f"structure the lattice is too coarse to seal; or set "
f"openfoam.allow_stepped_lid: true to build anyway (the verdict still "
f"reaches the result either way).")
if step_p99 > LID_STEP_WARN:
severe = step_p99 > LID_STEP_SEVERE
message = (
f"RIGID LID {'DOES NOT APPLY' if severe else 'MAY NOT APPLY'} HERE: "
f"at the 99th percentile the prescribed "
f"free surface steps {100 * step_p99:.0f}% of the local water depth "
f"within two cells. A lid is a WALL, so it cannot answer a drop by "
f"plunging the way a free surface does - it converts the head into "
f"velocity instead (sqrt(2 g dz), so a 1 m step drives 4.4 m/s). "
f"Expect the fastest cells to sit at those steps and the velocity cap "
f"to become load-bearing. Either use mode: vof where the surface has "
f"to move, or check that the steps are real rather than the 2D result "
f"read across a structure the lattice is too coarse to seal."
+ (" Built anyway: openfoam.allow_stepped_lid is set." if severe
else ""))
notes.append(message)
log.warning("%s", message)
# ---- bed layer sized against the roughness ------------------------------
bed_ks = _roughness_at(cfg, grid.cell_xy)
bed_layer = of.bed_layer_height
pin_bed = of.auto_bed_layer if of.auto_bed_layer is not None else (not rigid)
if bed_layer is None and pin_bed and bed_ks is not None:
# OpenFOAM's nutkRoughWallFunction places the first grid point on a log law
# whose origin is displaced into the roughness; that is only meaningful with
# the first cell CENTER above the roughness crests, y1 > ks. The center sits
# at half the layer, so the layer must be ~2*ks - which is why a river mesh
# cannot simply be refined towards the bed the way an aerofoil mesh can.
#
# The reference ks is the median over the WETTED columns, not over the whole
# footprint: the dry floodplain zones carry the vegetation ks (isar: 0.5 m),
# which has no bearing on a wall function under water and would otherwise
# pin the bed layer to half the column height everywhere.
wet_cols = (column_depth > cfg.hydrodynamics.wet_depth) if \
column_depth is not None else np.ones(bed_ks.shape, dtype=bool)
sample = bed_ks[wet_cols] if wet_cols.any() else bed_ks
bed_layer = 2.0 * float(np.nanmedian(sample))
# ...but it must still leave room for the layers above it, or the vertical
# distribution degenerates into one thick cell plus a stack of slivers.
ceiling = float(np.nanmedian(np.maximum(lid - bed, of.min_column_height))) / 3.0
if bed_layer > ceiling:
notes.append(
f"bed layer wanted {bed_layer:.3f} m (2x the wetted median ks) but is "
f"clipped to {ceiling:.3f} m to leave room for the layers above it - "
"the wall function is not admissible at this roughness/depth ratio "
"(see the quality report)")
bed_layer = ceiling
else:
notes.append(
f"bed layer pinned to {bed_layer:.3f} m (2x the wetted median ks, so "
"the first cell center clears the roughness crests)")
n_layers = resolve_layers(cfg, bed_ks, np.maximum(lid - bed, floor), notes=notes)
(points, internal, (bed_quads, bed_owner), (top_quads, top_owner), sides,
z, cell_centres, cell_column, cell_layer) = extrude(
grid, bed, lid, n_layers, bed_layer=bed_layer,
expansion=of.layer_expansion, min_height=floor)
labels, liquid_names, discharges = classify_sides(
cfg, sides["midpoint"], of.boundary_tolerance or 1.5 * of.cell_size)
boundary = [(BED_PATCH, "wall", bed_owner, bed_quads)]
banks = labels == BANKS_PATCH
# the height of each cell: the slot check of an outlet is about its LAYERS, which
# on a river lattice are a tenth of the plan spacing or less
level = z[column_corner_vertices(grid)].mean(axis=1)
cell_height = np.diff(level, axis=1)[cell_column, cell_layer]
for name in liquid_names:
sel = labels == name
if name.startswith(OUTLET_PREFIX) and of.outlet_air_margin is not None:
stage = _tailwater_for(cfg, name, tailwater)
if stage is not None:
owners = sides["owner"][sel]
keep, walled = split_outlet_above_tailwater(
owners, cell_centres, stage, of.outlet_air_margin,
cell_size=float(np.median(cell_height[owners])))
if walled.any():
index = np.flatnonzero(sel)
# the air part is a wall, so it joins the domain edge
banks[index[walled]] = True
sel = np.zeros_like(sel)
sel[index[keep]] = True
notes.append(
f"{name}: {int(walled.sum()):,} of {int(keep.size):,} faces "
f"are above the tailwater ({stage:.3f} m + "
f"{of.outlet_air_margin:g} m) and are walled, so air cannot "
"circulate in through the outlet and out through the lid")
boundary.append((name, "patch", sides["owner"][sel], sides["quads"][sel]))
# Split the solid structures out of the domain edge. A lateral face belongs to a
# structure when its midpoint sits on the blanked footprint; half a cell diagonal
# of slack, because the face lies on the boundary of the blanked columns rather
# than inside them, and the blanking is by cell-square INTERSECTION (above), so a
# face can stand up to that far from the solid it bounds.
edge, on_structure = split_structure_sides(sides["midpoint"], banks, blocked, dx)
boundary.append((BANKS_PATCH, "wall", sides["owner"][edge], sides["quads"][edge]))
if on_structure.any():
boundary.append((STRUCTURES_PATCH, "wall",
sides["owner"][on_structure], sides["quads"][on_structure]))
notes.append(
f"{int(on_structure.sum()):,} lateral faces bound a solid structure and "
f"are patched as '{STRUCTURES_PATCH}', apart from the {int(edge.sum()):,} "
f"on the domain edge - water touching a baffle is the flow doing its job, "
"water touching the edge means the footprint is too tight")
top_name = LID_PATCH if rigid else ATMOSPHERE_PATCH
# A rigid lid IS a boundary the flow cannot cross, so `wall` is right for it; the
# slip condition on U is what keeps it shear-free, as an air-water interface is.
#
# LEG C IS THE EXCEPTION, and getting this wrong cost an evening. Under
# potentialFreeSurfaceFoam the surface MOVES BY FLUX THROUGH THIS PATCH, and the
# `wall` type puts it in the `wall` group - so kOmegaSST applied wall functions
# and a zero wall distance at the free surface. The symptom was not subtle and
# also not diagnostic: 100% of cells velocity-limited at 6 m/s within three time
# steps, and SIGFPE inside kOmegaSSTBase::correct().
from . import potential
top_type = "patch" if (potential.is_potential(cfg) or not rigid) else "wall"
boundary.append((top_name, top_type, top_owner, top_quads))
mesh = assemble(points, internal, boundary, cell_centres)
problems = validate(mesh)
for problem in problems:
log.error("polyMesh: %s", problem)
if problems:
raise ValueError("the generated polyMesh is structurally invalid: "
+ "; ".join(problems))
result = OpenFoamMesh(
polymesh=mesh, grid=grid, n_layers=n_layers, z=z, bed=bed, lid=lid,
cell_centres=cell_centres, cell_column=cell_column, cell_layer=cell_layer,
column_bed=column_bed, column_wse=column_wse, column_depth=column_depth,
column_uv=column_uv, bed_ks=bed_ks,
rigid_lid=rigid,
lid_step_median=lid_step_median, lid_step_p99=lid_step_p99,
inlet_patches=[n for n in liquid_names if n.startswith(INLET_PREFIX)],
outlet_patches=[n for n in liquid_names if n.startswith(OUTLET_PREFIX)],
inlet_discharge=discharges, notes=notes,
)
log.info("OpenFOAM mesh: %d cells (%d columns x %d layers), %d faces",
mesh.n_cells, grid.n_columns, n_layers, mesh.n_faces)
return result