diff --git a/docs/README.md b/docs/README.md index 854bc68..de745f9 100644 --- a/docs/README.md +++ b/docs/README.md @@ -21,6 +21,7 @@ and lands in `design/` with the code that implements it. | Note | Status | |---|---| | [design/time-outside-space.md](design/time-outside-space.md) | **Implemented** 2026-10-26: archive format 2, one list of instants with the spaces inside each. | +| [design/domains-beyond-box.md](design/domains-beyond-box.md) | **Implemented** 2026-11-04 (C27 Phase 5): boundary regions by tag, periodic on a region whose faces match, `file` and `discrete` domains; decisions D1-D4. | ## History diff --git a/docs/archive-format.md b/docs/archive-format.md index 5cddd67..a8192be 100644 --- a/docs/archive-format.md +++ b/docs/archive-format.md @@ -398,6 +398,89 @@ per-rank grid metadata as these strings and stores them in fixed-width HDF5 datasets. (Parallel HDF5 rejects variable-length strings, so the width is agreed across ranks with an `allreduce` of the maximum length.) +## Meshes: facets, provenance and the cell type's order + +A mesh is a grid: a document with one instant whose one space is the mesh +(proteus names it `Mesh Spatial_Domain`, as in its run archives). Two +optional keys, on a uniform space or on any grid of a spatial collection, +carry what only a mesh needs: + +```yaml +topology: {type: Triangle, number_of_elements: 2, data_item: ...} +geometry: {type: XYZ, data_item: ...} +attributes: [...] +facets: + topology: {type: Polyline, number_of_elements: 4, nodes_per_element: 2, + data_item: ...} # node numbers of each facet + tags: {data_item: ...} # one integer per facet + view: global # or subdomain +provenance: + tool: a mesh converter # free text: the archive names no tools + version: '0.1' + files: [{name: square.msh, sha256: 3f5a...}] + extra: {...} # anything else the writer keeps +``` + +**Facets** are the tagged faces of the mesh (edges in 2D, points in 1D), +listed by their node numbers, with an integer tag each: **0 untagged, >0 a +condition**. The names of the conditions belong to the specification; the +mesh holds integers. A tag per facet, not per node, because a node at a +corner belongs to two faces that may carry different conditions. The facet +type follows the cell type and has the same order: `Polyline` facets for +`Triangle` and `Quadrilateral` cells, `Edge_3` for `Tri_6`, `Triangle` for +`Tetrahedron`, `Tri_6` for `Tet_10`, `Quadrilateral` for `Hexahedron` +(wedges and pyramids may use `Mixed`). +**-1 marks a face on the interface between two ranks.** It is the +partitioner's product, so a mesh as a whole (`view: global`, the default) +never stores it; one rank's piece (`view: subdomain`) may. `write_ymf` +refuses a -1 in a global view, and any tag below -1, when the tags are +inline; `ymf.mesh.check_facets(grid, base_dir)` reads tags held in an HDF5 +file and also checks that every facet node is a node of the grid and every +facet a face of some cell. + +**Provenance** says where the mesh came from: the tool, its version, and +each source file's name and SHA-256 (`ymf.mesh.file_provenance(tool, paths)` +computes them), so a reader can tell whether the mesh still matches the +files beside it. It takes no other keys; `extra` holds the rest. + +**Geometry order is the cell type's.** `Tri_6` is a quadratic triangle and +`Tet_10` a quadratic tetrahedron; there is no separate order field to +disagree with the type (`geometry_order(type)` reads it off). A quadratic +cell lists its vertices first, in the linear type's order, then one node on +each edge, in the order of `QUADRATIC_EDGES`: + +| type | edge nodes after the vertices sit on | +|---|---| +| `Edge_3` | (0,1) | +| `Tri_6` | (0,1) (1,2) (2,0) | +| `Quad_8` | (0,1) (1,2) (2,3) (3,0) | +| `Tet_10` | (0,1) (1,2) (2,0) (0,3) (1,3) (2,3) | +| `Pyramid_13` | (0,1) (1,2) (2,3) (3,0) (0,4) (1,4) (2,4) (3,4) | +| `Wedge_15` | (0,1) (1,2) (2,0) (3,4) (4,5) (5,3) (0,3) (1,4) (2,5) | +| `Hex_20` | (0,1) (1,2) (2,3) (3,0) (4,5) (5,6) (6,7) (7,4) (0,4) (1,5) (2,6) (3,7) | + +This is VTK's order, which is what a viewer draws: ParaView's XDMF readers +pass the connectivity through unpermuted. Sources do not all agree with it. +A common mesh generator numbers the last two edge nodes of a 10-node +tetrahedron (2,3) then (1,3); proteus's own quadratic tetrahedron numbers +them (0,1) (1,2) (2,3) (0,2) (1,3) (0,3). A writer states its source's +convention as data and permutes: + +```python +from ymf.mesh import reorder_cells, edge_midpoint_offsets +source = ((0, 1), (1, 2), (2, 0), (0, 3), (2, 3), (1, 3)) +cells = reorder_cells(cells, "Tet_10", source) # columns 8 and 9 swap +assert edge_midpoint_offsets("Tet_10", nodes, cells).max() < 1e-12 # straight-sided sample +``` + +`edge_midpoint_offsets` is the check: on a straight-sided mesh every edge +node sits at its edge's midpoint once the order is right, and about half an +edge away when two are swapped. + +`ymf2xmf` writes the grid as before; the facet list and the provenance are +not written to the `.xmf` (XDMF's `Set` is not supported, see +`docs/xdmf-model.yaml`). + ## What the writer checks, and what it doesn't YMF splits checking across a trust boundary: diff --git a/docs/changes.md b/docs/changes.md index 5bcd05e..eefa6c8 100644 --- a/docs/changes.md +++ b/docs/changes.md @@ -7,6 +7,165 @@ rediscover why a thing is the way it is. Open gaps are in ## 2026-11 +- **2026-11-04: `discrete` domains are posed (D4).** A mesh as ymf stores + one (`topology`, `geometry`, and the tagged `facets` and `provenance` of + cekees/ymf#6) is a domain: `adr_problem` poses it whole, with bounds from + its inline nodes (or `bounds`, or the solver's, for nodes in a file), its + boundary regions by facet tag and `mesh: {levels}` only; periodic is + refused. The design note, decided 2026-11-04, lands as + `docs/design/domains-beyond-box.md`. +- **2026-11-04: `file` domains are posed (D2, D3).** `adr_problem` carries + a file domain whole (`domain.region`) for the solver to read. Optional + `bounds: {lower, upper}` set the box; without them `lower` and `upper` + are None and the solver takes them from the mesh (posing follows meshing + for this kind alone), and the dimension must be given. A mesh file + (`2dm`, `3dm`, `hex`, `tetgen`, `ymf`) fixes the mesh: its discretization + gives `mesh: {levels}` only, keyed `levels=L`, and a mesh size is + reported as lacking; `gmsh_geo` is meshed by `h`. `ymf.discretizations`, + `ymf.closure.realizations` and `ymf.state` take the domain for this + (optional; without it a mesh needs a size, as before). Periodic + conditions on a file domain are refused. +- **2026-11-04: periodic conditions on a `piecewise_linear` region (D3).** + Refused before ("periodic faces are a box's"). Allowed now when, on each + periodic axis, the region's boundary segments (2-D) or facets (3-D) on the + lower face translate exactly onto those on the upper face (as point sets, + to 1e-8 of the box); otherwise refused with the count of pieces without a + partner. `ymf.symbolic.problem.periodic_faces_match`, used by both + readers. The solver's part (no boundary points added by the mesher, a + partner check after meshing) is Proteus's companion PR. +- **2026-11-04: boundary regions by tag (C27 Phase 5, decision D1).** + `boundary_regions` entries take `tags: [ints or names]` instead of + `where`, exactly one of the two; names resolve through the domain's + `boundary_tags`. `adr_problem` and `split_problem` give a tag region as + `flags` and no `where`. Checked on load: one of the two, every name + resolves, no tags on a box, a periodic region uses `where`. +- **2026-11-04: the shallow-water family.** `tests/fixtures/shallow_water`: + the model (depth and discharge; porosity, Manning friction and + Darcy-Forchheimer drag as coefficients of the state; `wall` and + `outflow` condition sets), Proteus's numerics (SW2DCV as one sub-problem, + its algebraic scheme in the free-form `discretization` block), a lake at + rest over a bump, and the Torino flume (dam break down a 3 % slope into a + rod canopy, from `torino_flume/dswe/torino.py`). The second family + needed no new schema: what two-phase flow added covers it. +- **2026-11-04: boundary conditions in `split_problem`.** Each sub-problem + carries the conditions it imposes, in `adr_problem`'s plain-data shape + (`where` code, per-component values, series, robin parts, periodic axes), + for every kind: those on its own unknowns, and those on a coefficient its + introduced unknown carries. In the dam break: flow gets the free-slip + walls, `p = 0` and an open `v` on top; VOF gets air flowing in through + `χh`; the level-set solves get none, as Proteus's `setAtmosphere` sets + none for them. +- **2026-11-03: gauges.** `model.gauges` in the shape the force-history + plan (C22, step 4) set out: `point`, `line` and `line_integral` gauges of + named fields, and `boundary_integral` gauges of wetted area, forces and + moments per boundary region, with `sample_rate` and `file`. What each kind + needs is checked on load. Also: a reference error now says "an analytical + solution", not "a analytical solution". +- **2026-11-03: condition sets.** `strong_form.condition_sets` names bundles + of conditions without a region (a model's `free_slip`, `atmosphere`), and + `strong_form.boundary_sets` applies them to regions. They are a separate + key rather than a `set:` entry among `boundary_conditions` because + strictyaml cannot accept either of two mapping shapes in one list, and + making a condition's keys optional would lose the schema error that names + the file holding an incomplete one. `ymf.boundary_conditions` expands + them, and `adr_problem`, `ymf.units` and the load-time checks all read + through it; a condition written out replaces a set's member for the same + region and variable. +- **2026-11-03: `slip`, `open`, `inflow`, `characteristic`; conditions on + a coefficient.** The boundary kinds two-phase flow and shallow water need: + a free-slip wall, an open boundary (C12's declared open boundary), a value + imposed where the flow enters, a stage imposed through the incoming + characteristics. A condition may name a coefficient of the state (`χ`). + Conditions are now checked on load: the region exists, the variable is an + unknown or a coefficient, `slip` is on a vector, and a value is given + exactly when the kind takes one. `ymf.units` checks `slip`, `inflow` and + `characteristic` like `dirichlet`. `adr_problem` refuses the four kinds + and conditions on coefficients, naming the solver they are for. +- **2026-11-03: each sub-problem classified on its own + (`ymf.symbolic.split.split_problem`).** The plain data of a split + discretization, one ADR form per sub-problem. Its own unknowns are + unknowns; what it reads from the others is a known field (a function of + the coordinates and time, written by name, its gradient + `grad__`), and must be declared in `lags`; a coefficient + lagged via the sub-problem that carries it is replaced by the carrier. + `binds` are substituted (a symbol with only a meaning, `h`, stays a + parameter the solver supplies), and `regularize` smooths `heaviside` and + `sign` exactly as Proteus does (`gf.H`, `gf.D`, `2 H_ε − 1`). Schema: + a sub-problem's `for` and `regularize`, a lag's `of`. The two-phase dam + break classifies in all five sub-problems: flow as variable-density + Navier-Stokes with ρ and μ from VOF's carrier, VOF and the level set as + transport in the lagged velocity (`v·∇φ` is a Hamiltonian in `φ` alone + there), redistancing as an eikonal Hamiltonian, mass correction as a + regularized reaction-diffusion. +- **2026-11-03: coefficients of the state keep their names.** + `adr_problem` lists each coefficient that follows the state under + `coefficients`: its numpy code, `depends_on`, and its derivative in each + unknown it depends on (`sharp: true` instead, for `heaviside` or `sign` + of an unknown). The ADR form keeps them substituted, so a solver of one + system is unaffected; the names are for a generator or a dispatcher. A + problem without such coefficients has no `coefficients` key, so its + record is unchanged. +- **2026-11-03: a discretization may be split into sub-problems.** + `subproblems:` (solved in order each step: what each `solves`, + `enforces`, `introduces`, `lags`, `writes_back`, with its own space, + time integration, stabilization and a `solver` that refines the + entry's), `binds:` (symbols such as `h` and `ε = 1.5 h` that the + discretization supplies), `regularize:` (sharp functions made smooth) and + `step_control:`. The structure and every name are checked on load + (`ymf.references`); `ymf.discretizations.subproblem_solver` gives a + sub-problem's effective solver. Classification per sub-problem is not + done yet. Splitting is not a new `method`: a split discretization is + still a mesh one, and its studies and output keys work as before. +- **2026-11-03: `heaviside`, `sign`, `min`, `max`, `norm`, `sym`.** The + two-phase model writes `χ = heaviside(φ)`, its initial signed distance + with `min`, `max` and `sqrt`, its viscous stress with `sym(∇v)`; + shallow-water friction uses `norm(q)`. They parse, print back + (`to_notation`), and become numpy code (`min`/`max` as nested + `numpy.minimum`/`numpy.maximum`, not `functools.reduce`). The derivative + of `heaviside` or `sign` is a Dirac delta, which has no value at a point: + where the ADR form needs one (a sharp function of an unknown), it is + refused with a message saying the discretization must regularize it. +- **2026-11-03: named equations.** `strong_form.equations` takes + `{name, for, equation}` beside bare strings. Bare strings are still + ordered like the unknowns; with `for` on every equation the order written + stops mattering. Names merge across `extends` like every other named + list, and are what a discretization's sub-problems will cite. Checked on + load: a name used twice, `for` naming something that is not an unknown, + two equations for one unknown, `for` on some but not all, and (outside an + abstract file) an unknown no equation is for. `ymf.equations_in_order` + is the one reader, shared by `ymf.symbolic` and `ymf.units`. +- **2026-11-03: a coefficient's formula may name the unknowns.** A + coefficient could be a formula in the coordinates, time and the other + coefficients; one that named an unknown was "could not be resolved". Now + it may: a state-dependent density, conductivity or friction law is a + coefficient whose formula names the unknowns, substituted into the + equations so the ADR form carries the dependence and its derivatives. + Decided by cekees for the model-family work (two-phase level-set flow: + `χ = heaviside(φ)`, `ρ(χ)`, `μ(χ)`; shallow water: Manning friction, + canopy drag), in place of a separate `relations:` block: one kind of + thing, a coefficient, whether it is constant, varies in space or follows + the state. What came with it: boundary data, initial conditions and exact + solutions that lean on such a coefficient are refused (they are evaluated + at points in x, y, z, t); `ymf.units` warns when a characteristic scale or + a dimensionless number is built on one; a coefficient may not take an + unknown's name; a formula may not take a derivative of an unknown (yet). + `ymf.symbolic.state_coefficients` lists which coefficients follow the + state, and on what. +- **2026-11-03: a mesh carries its tagged facets, its provenance, and its + order in its cell type.** Proteus's mesh campaign (C7.1, step A) needs a + mesh to say which faces carry which boundary condition, and where it came + from, so converters from other codes' inputs can write ymf. A grid may + now carry `facets` (node tuples and one integer tag per facet: 0 + untagged, >0 a condition, -1 a rank interface, allowed only in a + `subdomain` view) and `provenance` (tool, version, files with SHA-256; + no list of tools anywhere). Facets, not node tags, because a corner node + belongs to two faces. Geometry order stays the topology type's + (`geometry_order`), and the node order of each quadratic type is now + written down (`QUADRATIC_EDGES`, VTK's, what viewers draw); the + 10-node tetrahedra of one widely used generator and of proteus number + their edge nodes differently, and `ymf.mesh.edge_node_permutation` + maps any stated convention onto ymf's. Older documents are unchanged. + - **2026-11-01: one name per thing: `model` and `approximation`.** The spec's top-level key has been `model` since a66ea1c, and the archive of a spec stores the realized spec under the key the spec uses, so it has written diff --git a/docs/design/domains-beyond-box.md b/docs/design/domains-beyond-box.md new file mode 100644 index 0000000..04e019b --- /dev/null +++ b/docs/design/domains-beyond-box.md @@ -0,0 +1,179 @@ +# Domains beyond a box: `file`, `discrete`, and periodic on a non-box + +Status: proposal 2026-10-26 (branch `domains-beyond-box`, bed1608); decided by +cekees 2026-11-04 (§ 5); built in cekees/ymf#8 and, for Proteus, +cekees/proteus#188, except Proteus's reading of a `discrete` domain, which +waits on cekees/proteus#186 (§ 3). Sections 2-4 are the proposal as written, +each followed by what was decided and built. Scope is what was left of C27 +Phase 5: Neumann, Robin, tabulated hydrographs and `piecewise_linear` regions +already ran. + +## 1. Where things stand + +`adr_problem` (`symbolic/problem.py:_posed`) solves a `box`, or a +`piecewise_linear` region, and raises `SymbolicError` for any other kind, +naming it. It refuses periodic conditions when the domain is a region +(`problem.py`, "periodic faces are a box's"). Three facts shape the options. + +- **Proteus can already build these domains.** `DomainYMF.from_ymf` turns a + `kind: file` block into the matching `Domain` class (`gmsh_geo` and the + mesh-file formats in `_FILE_CLASSES`). A `discrete` block becomes a + `TriangulatedSurfaceDomain` only, and only for `Triangle` topology: a + surface, not a volume to solve on. The refusals are in ymf's layer. +- **Everything downstream of the domain assumes a box's bounds.** + `ADRProblem.physics` sets `L` and `x0` from `domain.lower/upper` (for a + region, the bounding box of its vertices), and that box sets the boundary + tolerance. A `file` domain has no vertices to take a box from until it is + meshed, and `_posed` runs before meshing. +- **Boundary regions are coordinate predicates only.** + `boundary_regions: [{name, where}]`, evaluated on the point. `getDBC(x, flag)` + is already handed the mesh's boundary flag, but nothing in the spec can + select by it. For a region that is fine (the predicate holds on any mesh of + it). For a CAD or Gmsh model, whose physical groups are the natural names, + writing coordinates is the wrong tool. + +## 2. `file` domains + +**Proposal: solve them, with boundary regions selectable by tag.** + +1. `boundary_regions` entries take `tags: [ints or names]` as an alternative to + `where` (exactly one of the two). Names resolve through the domain's + `boundary_tags` (Gmsh physical groups come through + `getGmshPhysicalGroups`, which `Domain.gmsh2proteus` already uses). Tag + selection is by the `flag` argument Proteus passes, so it needs no new + plumbing in `ADRProblem._dirichlet/_neumann/_robin`: a test of + `flag in tags`, beside the predicate. + This is a change to the schema and to composition (a region now has two + forms), so it belongs in the same PR as the first solve. +2. Bounds. A `file` block gains optional `bounds: {lower, upper}`. With it, + the box machinery is unchanged. Without it, `adr_problem` refuses with the + fix ("give bounds, or run `pro msh` first") rather than guessing. Reading + the bounds back from the generated mesh is nicer but moves a step that is + pure today (`_posed` builds the problem from the document alone) to after + meshing; do not do that until someone needs it. +3. The mesh. `discretizations[].mesh` is `cells` (box) or `h` (region). A + file's mesh is set by the file, so `mesh: {levels}` only, and `h`/`cells` + on a `file` domain is an error that says so. Multilevel from a file is + whatever the Gmsh route supports; check before promising `levels > 1`. + +Out of scope: formats Proteus has no `Domain` class for (`from_ymf` already +reports these by name), and STL surfaces, which are boundaries, not solvable +volumes. + +**Decided and built.** D1: `tags` beside `where`, exactly one per region, +names through `boundary_tags`; a tag region reaches the solver as `flags` and +no `where` (an older solver fails on it, `KeyError: 'where'`, rather than +taking it for the whole boundary). D2: `bounds` is optional; without it the +solver takes the bounds from the mesh, so posing follows meshing for this +kind alone (Proteus reads them from a tetgen `.node` or 2dm/3dm `ND` lines; +other formats give `bounds`). The mesh: a mesh file (`2dm`, `3dm`, `hex`, +`tetgen`, `ymf`) takes `mesh: {levels}` only, keyed `levels=L`; a geometry +(`gmsh_geo`) takes `h`, which Proteus's gmsh route already uses (`-clmax +0.5 he`), so item 3 above was half wrong. `levels > 1` works from a file: +every mesh-file class refines its coarse mesh uniformly +(`generateFromExistingCoarseMesh`), and a tetgen cube read as a file +converges at rate 1.94, 1.96 over levels 1-3. + +## 3. `discrete` domains + +**Superseded by D4 (cekees, 2026-11-04): a `discrete` domain is solved.** It +is a mesh as ymf stores one (C7.1, mesh decision 3: a mesh is a type of +domain): `topology`, `geometry` and, from cekees/ymf#6, tagged `facets` and +`provenance`. ymf poses it whole, with its bounds from its nodes when the +document holds them inline (from the solver, or `bounds`, when they are in an +HDF5 file), its boundary regions by its facet tags, and `mesh: {levels}` +only. Proteus reads it through `MeshYmfDomain` (cekees/proteus#186, the C7.1 +session's): that PR adds `ymf` as a file format; what it does not yet do is +turn an inline `discrete` block into a `MeshYmfDomain`, and the dispatch is +in `DomainYMF.py`, which #186 owns. Until it lands, Proteus refuses a +`discrete` block it cannot read, as before. + +The proposal, kept for the record: **keep them carried and refused, and make the refusal useful.** +The only consumer is `TriangulatedSurfaceDomain`, a surface. A volume mesh in +the document would have to be written out as `.node/.ele` (or `.msh`) and read +through `MeshTetgenDomain`/`Mesh2DMDomain`; that is a second path to a mesh +that `file` already gives for the price of one export step. Revisit if an +archive's own grids are wanted as a domain (the "replay on the mesh I got" +case), which is a closure question, not a Phase 5 one. The change now is the +message: say "write the mesh to a file and use `kind: file`". + +## 4. Periodic on a non-box domain + +Today's mechanism (`ADRProblem._periodic`) is a key function: a node on either +face of a periodic axis maps to its coordinates with that axis snapped to the +lower face, and Proteus gives nodes with equal keys one dof. That is purely +geometric, so it does not care whether the domain is a box. What it needs is a +mesh where **every node on one face has a partner on the other**. + +- `box` meshes (structured) have that for free. +- A `piecewise_linear` region is meshed by triangle/tetgen. If the boundary + segments on the two faces are translations of each other and the generator + adds no boundary points (triangle's `-Y`; tetgen's equivalent), the nodes + match. If it does add them, the faces get different nodes and periodicity + fails *silently* (unmatched nodes are simply not merged). + +**Proposal:** allow `periodic` on a `piecewise_linear` region when all hold, +and check them at validation, not by trusting the mesh: + +1. the periodic faces are bounding-box faces (the `_periodic_axes` test + already finds them, on the region's bounding box); +2. the boundary segments (2-D) or facets (3-D) on the lower face translate + onto those on the upper face, within tolerance; +3. the mesher is told not to refine the boundary. + +After meshing, assert that every keyed node has a partner, and fail with the +count if not; this check is cheap and is the one that catches (3) being wrong. +Non-axis-aligned or rotational periodicity is out of scope: it needs a shift +vector per pair rather than a snapped coordinate, and nobody has asked. +Periodic on a `file` domain is out of scope for the same reason plus Gmsh's +`Periodic` constraint, which Proteus's reader does not carry. + +**Decided and built.** ymf checks condition 2 at posing +(`periodic_faces_match`: the pieces on the lower face, moved by the box's +length, as point sets to 1e-8 of the box, with a count of those without a +partner). Proteus does 3 and the check after meshing: `Y` for triangle and +tetgen, and, since `Y` alone leaves a face with only its own vertices (2 +nodes on the notched channel's x face), every 2-D boundary segment cut into +equal pieces no longer than `h` before meshing, so translated faces get +translated nodes; `run` stops with the count of unpaired face nodes. On the +notched channel the old options left 0, 3 and 6 unpaired nodes at +h = 0.2, 0.1, 0.05; now 0, with P1 and P2 converging at rates 2 and 3 on a +solution whose x-derivative is nonzero on the periodic faces. 3-D faces are +not cut (tetgen's `Y` keeps the facets' own vertices); the check reports if +that is not enough. Periodic on `file` and `discrete` domains is refused. + +## 5. Order, and what to decide + +1. Non-box periodic (section 4): smallest, no schema change, one test with a + notched or L-shaped region that keeps two opposite faces. +2. `file` domains with tag selection (section 2): schema change plus + composition; the largest. +3. `discrete`: the message only. + +Decisions for cekees (answered 2026-11-04; the order was kept, periodic first, +in one round with `file` domains, D3): + +- **Tag selection in `boundary_regions`** (`tags:` beside `where`): yes/no. It is + the one schema change, and the whole value of `file` domains rests on it. + **D1: yes, exactly one of the two.** +- **`bounds` required on `file` domains** (refuse without), or read from the + mesh later. **D2: optional; without it, read from the mesh; posing stays + pure whenever bounds are given.** +- **Order**: periodic first, or `file` first. **D3: one round, periodic + first.** +- **D4 (new):** `discrete` domains solved through `MeshYmfDomain` and the + facet tags, bounds from the nodes; § 3 rewritten. + +## 6. Not verified (updated 2026-11-04) + +- ~~Whether Proteus passes `-Y`~~: it did not (`VApq30Dena`, + `VApq1.35q12feena`); it does now for a periodic region. +- ~~That `levels > 1` works from a `file` domain~~: it does (tetgen cube, + levels 1-3, rates 1.94 and 1.96). +- ~~Which classes `_FILE_CLASSES` covers~~: 2dm, 3dm, hex and tetgen, plus + `ymf` with cekees/proteus#186; Proteus reads bounds from tetgen, 2dm and 3dm. +- The cost of the partner check on a large 3-D mesh (it is two set + constructions over the face nodes; not timed). +- The `gmsh_geo` route with `h` and without `bounds` (Gmsh must make the mesh + before the bounds can be read; not exercised). +- Proteus solving a `discrete` domain: waits on cekees/proteus#186. diff --git a/docs/problem-specification.md b/docs/problem-specification.md index 231f1f0..7899ff6 100644 --- a/docs/problem-specification.md +++ b/docs/problem-specification.md @@ -14,6 +14,8 @@ that have been closed and why, [changes.md](changes.md). - [Writing it once: the notation](#writing-it-once-the-notation) - [Running a spec](#running-a-spec) - [Branching: one problem, several solution paths](#branching-one-problem-several-solution-paths) +- [Gauges](#gauges) +- [Splitting a discretization into sub-problems](#splitting-a-discretization-into-sub-problems) - [Document state](#document-state) - [Composing documents](#composing-documents) - [Provenance](#provenance) @@ -87,6 +89,27 @@ Optional: `initial_conditions`, `boundary_conditions`, `coefficients`. `dimensional_check` used to be accepted and are refused now; see [changes.md](changes.md).) +An equation is a string, and the strings are taken in the order of the +unknowns. Or it is named, and says which unknown it is `for`: + +```yaml +equations: + - name: continuity + for: p + equation: "∇·v = 0 in Ω" + - name: momentum + for: v + equation: "ρ ∂v/∂t + ∇·(ρ v⊗v) − ∇·(μ∇v) + ∇p = f in Ω" +``` + +With `for` on every equation the order written no longer matters (it is +`for` on all of them or on none). A name is unique, is what a file that +extends this one merges onto (replacing `momentum`'s text replaces that +equation and is recorded as an override), and is how a discretization will +say which equations each of its sub-problems solves. An abstract model may +leave an unknown without an equation for the files that extend it to +supply; a well-posed document may not. + ```yaml initial_conditions: - field: T @@ -112,7 +135,51 @@ as `value` or `formula` like any other right-hand side. All three are expressions, so they may vary along the boundary. It is one relation, so a vector unknown wants one condition per component. `beta: 0` is refused: that is a Dirichlet condition, and belongs where it can be imposed strongly. -`periodic` conditions need a box (see [The domain](#the-domain)). +`periodic` conditions need a box, or a region whose two faces are +translations of each other (see [The domain](#the-domain)). + +Four more kinds are for a [split discretization](#splitting-a-discretization-into-sub-problems)'s +solver; one ADR system (`adr_problem`) runs none of them and says so: + +| type | means | value | +|---|---|---| +| `slip` | on a vector: its normal component is the value (0 if none is given), its tangential traction zero -- a free-slip wall | optional | +| `open` | a declared open boundary: the natural flux applies, nothing is imposed | none | +| `inflow` | the value, imposed only where the flow enters -- the phase that comes in through an open top | required | +| `characteristic` | a state imposed through the incoming characteristics only -- a tailwater stage that a supercritical outflow ignores | required | + +A condition may also name a coefficient of the state rather than an +unknown: `{variable: χ, type: inflow, value: 1.0}` says the fluid entering +is air, whatever carries `χ` in a given discretization. Checked on load: +the region exists, the variable is an unknown or a coefficient, `slip` is +on a vector, and a kind gets a value exactly when it takes one. + +**Condition sets.** What a free-slip wall or an open top means is part of +the model: written once there, over its own unknowns and coefficients, and +applied to regions by each problem. + +```yaml +# the model +condition_sets: + free_slip: + - {variable: v, type: slip} + atmosphere: + - {variable: p, type: dirichlet, value: 0.0} + - {variable: v, type: open} + - {variable: χ, type: inflow, value: 1.0} # what enters is air +# a problem extending it +boundary_sets: + - {region: walls, set: free_slip} + - {region: top, set: atmosphere} +``` + +`ymf.boundary_conditions(strong_form)` is every condition with the sets +expanded, each member carrying its region and the set it came from. A +condition written out under `boundary_conditions` for the same region and +variable replaces the set's member, so a case can apply `atmosphere` and +give the pressure its own value; two sets that both decide one region's +variable are an error, as is a set nobody defines. A set's members are +checked where they are defined, as any condition is. A value is given as a number (`value`), an expression (`formula`), or a table, `series: [[t, value], ...]` with `t` ascending: a measured or design @@ -260,6 +327,8 @@ solution_paths: | `x² π³ x^n x**n` | powers | | `2π²κt`, `κ ΔT`, `sin(πx)` | implicit multiplication | | `sin cos tan exp log sqrt sinh cosh tanh abs` | functions, with parentheses | +| `heaviside(φ)`, `sign(e)`, `min(a, b, ...)`, `max(a, b, ...)` | sharp functions; a solver needs `heaviside` and `sign` of an unknown regularized, since their derivatives are Dirac deltas | +| `norm(v)`, `sym(∇v)` | `sqrt(v·v)`, and the symmetric part `(A + Aᵀ)/2` of a square tensor | | `grad div lap dt dx dy dz dot outer` | the operators spelled out, for ASCII-only text | | `lhs = rhs in Ω` | an equation; the `in ...` clause is ignored | | `y = 0 or x = 0 and y = H` | a boundary region; `∂Ω` is the whole boundary | @@ -277,6 +346,44 @@ unknown name 'kapa' (not a coefficient, an unknown, a coordinate, t, π or a function, nor a product of single-letter names) at position 2 in '2 kapa T' ``` +**A coefficient may be a function of the state.** Its formula may name the +unknowns as well: a conductivity that grows with temperature, a density +that follows a phase indicator, a friction law in the depth and discharge. + +```yaml +strong_form: + unknowns: [T] + equations: + - "∂T/∂t = ∇·(κ∇T) in Ω" + coefficients: + κ: {formula: "κ0 (1 + T²)", units: m2/s} + κ0: {value: 1.0e-3, units: m2/s} +``` + +It is substituted into the equations like any other coefficient, so the ADR +form carries the dependence: `κ` here makes the diffusion of `T` nonlinear +in `T`, with its derivative, exactly as if `κ0 (1 + T²)` were written in the +equation. Write a variable coefficient inside the divergence (`∇·(κ∇T)`, +`∇·(μ∇v)`); `κΔT` is the divergence form only while `κ` is constant, and is +refused otherwise. A formula may name the unknowns, not their derivatives. + +Such a coefficient is not a number, and three things that need one say so +rather than go quiet. Boundary data, initial conditions and exact solutions +are evaluated at points in `x`, `y`, `z`, `t`, so one that names it is an +error (the heat example's exact solution names `κ`, and stops being exact +the moment `κ` depends on `T`). A characteristic scale or a dimensionless +number built on it has no value, so `ymf.units` warns and asks for a +constant coefficient instead (`ρw`, not `ρ = ρw (1 − χ) + ρa χ`). And +`coefficient_numbers` leaves it out. A coefficient may not share its name +with an unknown. + +`adr_problem` substitutes such a coefficient, and also lists it by name +under `coefficients` -- its code, what it depends on, and its derivative +in each unknown (`sharp: true` and no derivative for `heaviside` or `sign` +of an unknown) -- so a code generator can emit `ρ` once as a function of +the state rather than its expansion in every term. A problem with no such +coefficient has no `coefficients` key and reads exactly as before. + Tensors follow the continuum-mechanics convention `grad(v)[i, j] = ∂v_j/∂x_i`, so `v·∇v` is (v·∇)v and `∇·(μ∇v)` is the vector Laplacian. @@ -415,10 +522,57 @@ two levels, as with strong and weak forms. `kind: box` and `kind: piecewise_linear` can be solved today: a box is meshed by `cells` per side, any other region by `h`, the element size of an -unstructured mesh (see the discretization's `mesh`). Boundary regions are -still named by `where` predicates, which hold on any mesh of the region; -periodic conditions need a box. `discrete` and `file` domains validate and -are carried, and `adr_problem` says so rather than failing obscurely. +unstructured mesh (see the discretization's `mesh`). A periodic condition +on a `piecewise_linear` region holds when the region's boundary pieces on +the two faces of each periodic axis are translations of each other: the +segments (2-D) or facets (3-D) on the lower face, moved by the length of +the bounding box along the axis, must be exactly those on the upper face. +That is checked at posing, with a count of the pieces that have no partner; +the mesher must then add no points to them, and the solver checks after +meshing that every node on one face has its partner on the other. + +A `file` domain is carried whole for the solver to read. A mesh file +(`2dm`, `3dm`, `hex`, `tetgen`, `ymf`) fixes the mesh, so its discretization +gives `mesh: {levels: ...}` only, the levels of uniform refinement of the +file's mesh, and a mesh size there is reported as what it lacks; a geometry +(`gmsh_geo`) is meshed by `h` like a region. `bounds: {lower, upper}` is +optional: it sets the scale of the boundary tolerance, and without it the +solver takes the bounds from the mesh, so for this kind alone posing +follows meshing; the dimension must then be given. A file's boundary +regions are naturally its boundary tags. A file domain cannot be periodic +(Gmsh's `Periodic` constraint is not carried). + +A `discrete` domain is a mesh as ymf stores one (`topology`, `geometry`, and +optionally tagged `facets` and `provenance`; see `ymf.mesh`): a mesh is a +type of domain. It is posed whole, like a file. Its bounds are its nodes' +when the geometry holds them inline; when they are in an HDF5 file the +solver takes them, or `bounds` gives them, and the dimension comes from the +geometry's type (`XY`, `XYZ`). Its boundary regions are its facets' tags, a +discretization gives `mesh: {levels}` only, and it cannot be periodic. The +design is [design/domains-beyond-box.md](design/domains-beyond-box.md). + +**Boundary regions by position or by tag.** A region is named by `where`, a +predicate on the point that holds on any mesh of the domain, or by `tags`, +the mesh's boundary tags: ints, or names the domain's `boundary_tags` +resolve. Exactly one of the two: + +```yaml +domain: + kind: piecewise_linear + vertices: [[0, 0], [2, 0], [2, 1], [1, 1], [1, 2], [0, 2]] + segments: [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]] + segment_flags: [1, 2, 3, 3, 4, 5] + boundary_tags: {bottom: 1, right: 2, notch: 3, top: 4, left: 5} +boundary_regions: + - {name: dry, tags: [bottom, left]} + - {name: outlet, where: "x = 2"} +``` + +A tag region reaches a solver as `flags: [1, 5]`, which it tests against the +boundary flag of the point, and carries no `where`, so a solver that does +not read `flags` fails on it rather than taking a missing predicate for the +whole boundary. A box carries no tags (name its faces by position), and a +periodic region is found by position, so it needs `where`. ## Branching: one problem, several solution paths @@ -494,6 +648,135 @@ Adding a branch adds entries; it never edits existing ones. That is what makes a document safe for several contributors, or an LLM and a person, to extend. +## Gauges + +What a run measures as it goes is part of the problem, beside its +conditions: `model.gauges`, one entry per gauge file. + +```yaml +gauges: + - name: heights # the level set along four verticals + kind: line + fields: [φ] + lines: [[[2.724, 0.0], [2.724, 1.8]], [[2.228, 0.0], [2.228, 1.8]]] + file: height1.csv + - name: pressure + kind: point + fields: [p] + points: [[3.22, 0.16], [3.22, 0.584]] + - kind: boundary_integral # forces and moments per boundary region + tags: [cylinder] + quantities: [force, moment] + moment_point: [0.2, 0.2, 0.0] + sample_rate: 0.01 +``` + +`point`, `line` and `line_integral` gauges sample `fields` (unknowns or +coefficients, `ρ` as well as `p`) at `points` or along `lines`, each line +its two ends. A `boundary_integral` gauge integrates `quantities` +(`wetted_area`, `force_pressure`, `force_viscous`, `force`, `moment`, the +moment about `moment_point`) over each boundary region in `tags`. Checked +on load: what each kind needs and may not have, that fields and tags +exist, and that a gauge's points share one dimension. A solver reads them; +Proteus's reader is C22's next step. + +## Splitting a discretization into sub-problems + +A model with several unknowns is often solved as a sequence of smaller +problems in each time step. Two-phase level-set flow is the standard case: +the model is momentum, continuity and the interface transport of a level +set `φ`, with `ρ` and `μ` coefficients of `χ = heaviside(φ)`; Proteus solves +it as flow, then a volume-of-fluid carrier of `χ`, then the level set, then +redistancing, then a mass correction. That splitting is a property of the +discretization, not of the model, so another discretization of the same +model can split it differently (CLSVOF and a projection scheme) or not at +all. [`tests/fixtures/two_phase`](../tests/fixtures/two_phase) holds the +model and the dam break split this way; [`tests/fixtures/shallow_water`](../tests/fixtures/shallow_water) +holds the shallow-water family the same way (one sub-problem, as Proteus's +SW2DCV solves the whole system at once), with a lake at rest and the +Torino flume as cases. + +```yaml +discretizations: + - name: level_set_five + provenance: proteus_derived + mesh: {h: 0.0322} + binds: # symbols the model leaves to the discretization + h: {meaning: element_size} + ε: {formula: "1.5 h"} + regularize: # sharp model functions made smooth + heaviside: {by: smoothed_heaviside, width: ε} + step_control: {type: min_cfl, cfl: 0.25} + solver: {type: nonlinear, nonlinear_solver: newton, linear_solver: petsc} + subproblems: # solved in this order every step + - name: flow + solves: [momentum, continuity] + lags: {φ: previous_step, χ: {via: vof, at: previous_step}} + finite_element: {velocity: {family: CG, order: 1}, pressure: {family: CG, order: 1}} + solver: {max_iterations: 100} # refines the entry's solver + - name: vof + introduces: {name: χh, approximates: χ} + equation: "∂χh/∂t + ∇·(v χh) = 0" + lags: {v: {via: flow, at: current_step}} + - name: ncls + solves: [interface] + lags: {v: {via: flow, at: current_step}} + - name: rdls + enforces: "norm(∇φ) = 1" + pseudo_time: steady + writes_back: {φ: replace} + - name: mcorr + enforces: "heaviside(φ) = χh" + introduces: {name: φc, correction_of: φ} + writes_back: {φ: add} +``` + +A sub-problem **solves** model equations by name (so the model's equations +must be [named](#document-structure)), or **enforces** a constraint the +others lose, or **introduces** an unknown of its own (an auxiliary carrier +of a model quantity, or a correction to a model unknown), or spins the run +up. It says what it **lags** -- reads from elsewhere -- and when: as the +previous step left it, or as an earlier sub-problem of this step left it +(`via`). What it **writes back** into another field it says too. Its +`solver` refines the entry's, which stays required as the default. +`binds` names the symbols the discretization supplies and the model may not +use (the element size, a smoothing width), and `regularize` the sharp +functions it replaces with smooth ones (`heaviside`, `sign`). + +Checked on load: sub-problem names are unique; every model equation is +solved by exactly one sub-problem; every name a sub-problem uses (an +equation, a lagged field, a `via`, what it introduces, approximates, +corrects or writes back) exists, and an introduced or bound name does not +shadow one that does. The split is part of what an output's key covers, so +changing a sub-problem's solver plans a new output. + +A sub-problem may give an `equation` of its own, `for` the unknown it +determines (by default the one it introduces); a lag with `of` reads a +copy of a field under a name of its own (redistancing reads `φs`, the +level set before it changes `φ`); and a sub-problem's own `regularize` +replaces the entry's (redistancing smooths its sign over `0.75 h`, the +rest their Heaviside over `1.5 h`). + +**Each sub-problem is classified on its own.** `adr_problem` reads a model +as one system, which this one is not (`v·∇φ` is a Hamiltonian in another +unknown). `ymf.symbolic.split.split_problem(doc, "level_set_five")` gives +each sub-problem's ADR form instead: its unknowns are those of the +equations it solves and of its own equation; every other field it reads is +a known field, declared in its `lags` (a read that is not is an error that +says so), which the plain data names (`v_0`) and whose gradient it names +`grad_v_0_0`; a coefficient lagged through the sub-problem that carries it +(`χ` via `vof`, which introduces `χh`) is that carrier. Bound symbols are +substituted, and the sharp functions are smoothed as Proteus smooths them: +`heaviside` becomes `gf.H` of `equivalent_polynomials.h`, whose derivative +is its `gf.D`, and `sign` becomes `2 H_ε − 1`. Each sub-problem also +carries what it reads and when, its named coefficients, the +regularizations it used, an introduced unknown's initial condition, its +effective solver, and its **boundary conditions**: those on the unknowns it +solves for, and those on a coefficient its introduced unknown carries +(atmosphere's `χ` inflow is VOF's `χh` inflow), with condition sets +expanded. A condition on a field a sub-problem only reads is another's. A solver maps a sub-problem onto its own modules (a +dispatcher); `proteus:` blocks are interim hints for it, carried as text. + ## Document state A YMF document is the record of a progression: from a model, through a @@ -845,13 +1128,14 @@ These are real limitations of the current schema and code, recorded here so nobody has to rediscover them. Gaps that have been closed, and how, are in [changes.md](changes.md). -- **Only `box` and `piecewise_linear` domains are solved.** A domain of any - kind validates and is carried (see [The domain](#the-domain)), but - `adr_problem` refuses `discrete` and `file`, saying which kind it was - given rather than failing obscurely. -- **Periodic conditions only on a box.** A periodic region must name two - opposite faces of the box; on a `piecewise_linear` domain the condition is - refused. +- **A solver may not read every domain kind yet.** ymf poses all four; + Proteus reads `discrete` domains once cekees/proteus#186 lands, and takes + the bounds of a bounds-free file only from tetgen, 2dm and 3dm files. +- **Periodic conditions only along the axes.** A periodic region names + two faces of the domain's bounding box, on a box or on a region whose two + faces match; a rotation, or a shift that is not along an axis, is not + expressible, and periodic conditions on a `file` domain are not + supported (Gmsh's `Periodic` constraint is not carried). - **The units of a flux condition are partly checked.** `neumann` and `robin` see the diffusive flux only, and `t` has no declared units anywhere, so an `α` or `β` that depends on time is never decided. diff --git a/docs/xdmf-model.yaml b/docs/xdmf-model.yaml index 4e5d767..98b0c34 100644 --- a/docs/xdmf-model.yaml +++ b/docs/xdmf-model.yaml @@ -416,9 +416,12 @@ XDMF_Elements: Set: status: not_supported notes: |- - Boundary and material regions are archived today as ordinary - attributes (Proteus writes nodeMaterialTypes and - elementMaterialTypes as Node- and Cell-centred Int fields). + Material regions are archived today as ordinary attributes (Proteus + writes nodeMaterialTypes and elementMaterialTypes as Node- and + Cell-centred Int fields). Tagged boundary faces are a mesh's facets + key (node tuples and one integer tag per facet, see + archive-format.md, "Meshes"), which is ymf's own and is not written + to the .xmf. Time: ymf: step.Time diff --git a/examples/README.md b/examples/README.md index b9b495d..db16c55 100644 --- a/examples/README.md +++ b/examples/README.md @@ -17,6 +17,7 @@ python examples/read_archive.py /tmp/ymf-demo/heat.ymf # read one back | [poisson.ymf](poisson.ymf) | The smallest complete spec: one weak form, one discretization, an exact solution. v0.1-style bare unknowns. | | [poisson_v02_enriched.ymf](poisson_v02_enriched.ymf) | The same problem with v0.2 additions: units, a standard name, characteristic scales, a typed boundary condition. | | [heat_equation.ymf](heat_equation.ymf) | A time-dependent problem: initial and boundary conditions with units, a time integrator, an analytical solution. The archive examples embed this spec. | +| [nonlinear_diffusion.ymf](nonlinear_diffusion.ymf) | A coefficient that is a function of the state: the conductivity `k = 1 + u²` names the unknown, so the problem is nonlinear and solved by Newton. Manufactured exact solution; P1 and P2 rates 2 and 3. | | [kovasznay_flow.ymf](kovasznay_flow.ymf) | Branching: Navier–Stokes with two weak forms (Taylor–Hood, stabilized equal-order), a discretization for each, a Reynolds number computed from the scales, and a verification plan against the exact solution. | ### Composed from several files: [navier_stokes/](navier_stokes/) diff --git a/examples/nonlinear_diffusion.ymf b/examples/nonlinear_diffusion.ymf new file mode 100644 index 0000000..838be65 --- /dev/null +++ b/examples/nonlinear_diffusion.ymf @@ -0,0 +1,59 @@ +# ============================================================ +# YMF Example: steady nonlinear diffusion, -∇·((1 + u²)∇u) = f +# +# Demonstrates: a coefficient that is a function of the state. k names +# the unknown u, so the diffusion is nonlinear in u and the solver needs +# Newton; the ADR form carries k's derivative in u. The source term is +# manufactured from u = sin(πx) sin(πy), so the run has an exact answer. +# +# With Proteus's ymf_run: P1 converges at rate 1.99, P2 at 3.00 (8, 16, +# 32 cells), and --check reproduces every output bitwise. +# ============================================================ + +model: + name: "Nonlinear diffusion, conductivity 1 + u², unit square" + physical_model: + provenance: human_specified + processes: [heat_conduction] + assumptions: [steady] + strong_form: + provenance: human_specified + unknowns: [u] + equation_formulation: "steady nonlinear diffusion" + equations: + - "-∇·(k∇u) = f in Ω" + domain: "Ω = [0, 1] × [0, 1]" + boundary_regions: + - {name: all_boundaries, where: "∂Ω"} + boundary_conditions: + - {region: all_boundaries, variable: u, type: dirichlet, value: 0.0} + coefficients: + k: {formula: "1 + u²"} + s: {formula: "sin(πx) sin(πy)"} + f: {formula: "2π² (3 s² − sin(πx)² − sin(πy)² + 1) s"} + weak_forms: + - label: galerkin + provenance: human_specified + derivation: "multiply by v in H1_0, integrate the diffusion by parts" + solution_spaces: {trial: H1, test: H1_0} + bilinear: "a(u; u, v) = ∫Ω (1 + u²) ∇u·∇v dx" + linear: "l(v) = ∫Ω f v dx" + stabilization_method: "none" +solution_paths: + analytical: + - name: manufactured + provenance: human_specified + method: manufactured + solution: + formula: "u = sin(πx) sin(πy)" + discretizations: + - name: P1 + provenance: human_specified + finite_element: {fields: {family: CG, order: 1}} + mesh: {cells: [8, 16, 32]} + solver: {type: nonlinear, nonlinear_solver: newton, linear_solver: superlu, tolerance: 1.0e-10} + - name: P2 + provenance: human_specified + finite_element: {fields: {family: CG, order: 2}} + mesh: {cells: [8, 16, 32]} + solver: {type: nonlinear, nonlinear_solver: newton, linear_solver: superlu, tolerance: 1.0e-10} diff --git a/src/ymf/__init__.py b/src/ymf/__init__.py index 7850b22..2896bea 100644 --- a/src/ymf/__init__.py +++ b/src/ymf/__init__.py @@ -41,7 +41,8 @@ validate_domain, write_ymf, ) -from ymf.normalize import normalize_unknowns +from ymf.normalize import (boundary_conditions, equations_in_order, normalize_equations, + normalize_unknowns) from ymf.xdmf import ( build_xdmf_tree, grid_element_time, @@ -104,6 +105,9 @@ "grid_element_time", # --- core: stdlib helpers --- "normalize_unknowns", + "normalize_equations", + "equations_in_order", + "boundary_conditions", # --- lazy: front-end validation (see _LAZY_ATTRS) --- *_LAZY_ATTRS, ] diff --git a/src/ymf/archive.py b/src/ymf/archive.py index ffe9026..5623827 100644 --- a/src/ymf/archive.py +++ b/src/ymf/archive.py @@ -151,12 +151,19 @@ "DATA_TYPES", "DATA_ITEM_FORMATS", "GEOMETRY_TYPES", + "FACET_VIEWS", + "QUADRATIC_EDGES", + "geometry_order", + "nodes_per_cell", "data_item", "data_item_for", "attribute", "topology", "geometry", "grid", + "facets", + "provenance", + "check_facet_tags", "new_approximation", "add_instant", "uniform_grid", @@ -286,6 +293,94 @@ class YmfArchiveError(ValueError): #: Legal ``Geometry`` types. GEOMETRY_TYPES = frozenset({"XYZ", "XY", "X_Y_Z", "X_Y", "VXVYVZ", "ORIGIN_DXDYDZ", "ORIGIN_DXDY"}) +#: Where a grid's facet list sits: ``"global"``, a mesh as a whole, whose +#: tags are 0 (untagged) or a condition (>0); or ``"subdomain"``, one rank's +#: piece, whose tags may also be -1 for a face on the interface between +#: ranks. -1 is a partitioner's product, so a global mesh never stores it. +FACET_VIEWS = frozenset({"global", "subdomain"}) + +#: Vertices of each unstructured cell type with a fixed node count. +_VERTICES = { + "Edge_3": 2, + "Triangle": 3, "Tri_6": 3, "Triangle_6": 3, + "Quadrilateral": 4, "Quad_8": 4, "Quadrilateral_8": 4, + "Tetrahedron": 4, "Tet_10": 4, "Tetrahedron_10": 4, + "Pyramid": 5, "Pyramid_13": 5, + "Wedge": 6, "Wedge_15": 6, + "Hexahedron": 8, "Hex_20": 8, "Hexahedron_20": 8, +} + +#: The node order of the quadratic cell types. The vertices come first, in +#: the linear type's order; then one node on each edge, node +#: ``n_vertices + k`` sitting on the edge between the two vertices +#: ``QUADRATIC_EDGES[type][k]``. This is VTK's order, which is what a +#: viewer draws, since ParaView's XDMF readers pass connectivity through +#: unpermuted. Mesh generators do not all agree with it -- one widely used +#: generator numbers the last two edge nodes of a 10-node tetrahedron the +#: other way round, (2,3) before (1,3) -- so a writer whose source orders +#: edge nodes differently permutes them first; see +#: :func:`ymf.mesh.edge_node_permutation`. +QUADRATIC_EDGES = { + "Edge_3": ((0, 1),), + "Tri_6": ((0, 1), (1, 2), (2, 0)), + "Quad_8": ((0, 1), (1, 2), (2, 3), (3, 0)), + "Tet_10": ((0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)), + "Pyramid_13": ((0, 1), (1, 2), (2, 3), (3, 0), (0, 4), (1, 4), (2, 4), (3, 4)), + "Wedge_15": ((0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (0, 3), (1, 4), (2, 5)), + "Hex_20": ((0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), + (0, 4), (1, 5), (2, 6), (3, 7)), +} +for _long, _short in (("Triangle_6", "Tri_6"), ("Quadrilateral_8", "Quad_8"), + ("Tetrahedron_10", "Tet_10"), ("Hexahedron_20", "Hex_20")): + QUADRATIC_EDGES[_long] = QUADRATIC_EDGES[_short] +del _long, _short + +#: The facet types a cell type's facet list may hold. A 1D cell's facets are +#: points, a 2D cell's are edges and a 3D cell's are faces, of the cell's +#: own geometry order; a wedge or pyramid has faces of two shapes, so its +#: facet list may be ``Mixed``. +_FACET_TYPES = { + "Polyline": {"Polyvertex"}, "Edge_3": {"Polyvertex"}, + "Triangle": {"Polyline"}, "Quadrilateral": {"Polyline"}, "Polygon": {"Polyline"}, + "Tri_6": {"Edge_3"}, "Triangle_6": {"Edge_3"}, + "Quad_8": {"Edge_3"}, "Quadrilateral_8": {"Edge_3"}, + "Tetrahedron": {"Triangle"}, "Hexahedron": {"Quadrilateral"}, + "Wedge": {"Triangle", "Quadrilateral", "Mixed"}, + "Pyramid": {"Triangle", "Quadrilateral", "Mixed"}, + "Tet_10": {"Tri_6", "Triangle_6"}, "Tetrahedron_10": {"Tri_6", "Triangle_6"}, + "Hex_20": {"Quad_8", "Quadrilateral_8"}, "Hexahedron_20": {"Quad_8", "Quadrilateral_8"}, + "Wedge_15": {"Tri_6", "Triangle_6", "Quad_8", "Quadrilateral_8", "Mixed"}, + "Pyramid_13": {"Tri_6", "Triangle_6", "Quad_8", "Quadrilateral_8", "Mixed"}, +} + +_SHA256 = re.compile(r"^[0-9a-f]{64}$") + + +def geometry_order(topology_type: str) -> Optional[int]: + """The geometry order a topology type carries: 1 linear, 2 quadratic. + + The order is the type's, never a separate field, so the two cannot + disagree: a ``Tet_10`` grid is a quadratic tetrahedral mesh. ``Mixed`` + has no single order and gives ``None``. + """ + if topology_type not in KNOWN_TOPOLOGY_TYPES: + raise YmfArchiveError( + f"{topology_type!r} is not a known topology type. " + f"Known types: {sorted(KNOWN_TOPOLOGY_TYPES)}") + if topology_type == "Mixed": + return None + return 2 if topology_type in QUADRATIC_EDGES else 1 + + +def nodes_per_cell(topology_type: str) -> Optional[int]: + """The node count a cell of this type has, or ``None`` where the type does + not fix it: ``Polyvertex``, ``Polyline`` and ``Polygon`` (which say so in + ``nodes_per_element``), ``Mixed`` and the structured types.""" + if topology_type not in _VERTICES: + return None + return _VERTICES[topology_type] + len(QUADRATIC_EDGES.get(topology_type, ())) + + # XDMF's own documented default for DataItem/@Precision. Faithful to the # format rather than to what a solver usually writes -- use data_item_for() # to infer the right value from an actual array instead of relying on this. @@ -550,18 +645,81 @@ def geometry(item: Dict[str, Any], *, geometry_type: str = "XYZ") -> Dict[str, A return {"type": geometry_type, "data_item": item} +def facets( + topo: Dict[str, Any], + tags: Dict[str, Any], + *, + view: str = "global", +) -> Dict[str, Any]: + """Describe a grid's tagged facets: boundary faces (or edges, in 2D) and + an integer tag on each. + + ``topo`` is a :func:`topology` of the facets, as tuples of the grid's + node numbers (a triangle's facets are ``Polyline`` with + ``nodes_per_element=2``, a ``Tet_10``'s are ``Tri_6``); ``tags`` is a + data item of one integer per facet: 0 untagged, >0 a condition. The + names of conditions belong to the specification, not to the mesh. + ``view`` is ``"global"`` for a whole mesh, or ``"subdomain"`` for one + rank's piece, where -1 marks a face on the interface between ranks + (:data:`FACET_VIEWS`). A face per facet, not a tag per node: a node at + a corner belongs to two faces that may carry different conditions. + """ + return {"topology": topo, "tags": {"data_item": tags}, "view": view} + + +def provenance( + tool: str, + *, + version: Optional[str] = None, + files: Iterable[Dict[str, str]] = (), + extra: Optional[Dict[str, Any]] = None, +) -> Dict[str, Any]: + """Describe where a grid came from: the tool that made it, its version, + and the files it was made from with their SHA-256 checksums. + + ``files`` is a list of ``{"name": ..., "sha256": ...}``; + :func:`ymf.mesh.file_provenance` computes the checksums from paths. The + tool is free text -- the archive names no tools -- and ``extra`` is a + mapping the writer may fill with anything else it wants kept. + """ + block: Dict[str, Any] = {"tool": str(tool)} + if version is not None: + block["version"] = str(version) + block["files"] = [dict(f) for f in files] + if extra: + block["extra"] = dict(extra) + return block + + +def _mesh_keys(facet_list: Optional[Dict[str, Any]], + source: Optional[Dict[str, Any]]) -> Dict[str, Any]: + keys: Dict[str, Any] = {} + if facet_list is not None: + keys["facets"] = facet_list + if source is not None: + keys["provenance"] = source + return keys + + def grid( topo: Dict[str, Any], geom: Dict[str, Any], attributes: Optional[Iterable[Dict[str, Any]]] = None, *, name: Optional[str] = None, + facets: Optional[Dict[str, Any]] = None, + provenance: Optional[Dict[str, Any]] = None, ) -> Dict[str, Any]: - """Assemble one uniform grid: connectivity, coordinates, and fields.""" + """Assemble one uniform grid: connectivity, coordinates, and fields. + + A mesh may also carry its tagged ``facets`` (:func:`facets`) and its + ``provenance`` (:func:`provenance`). + """ g: Dict[str, Any] = {"topology": topo, "geometry": geom} if name is not None: g["name"] = str(name) g["attributes"] = list(attributes) if attributes is not None else [] + g.update(_mesh_keys(facets, provenance)) return g @@ -583,10 +741,16 @@ def uniform_grid( topo: Dict[str, Any], geom: Dict[str, Any], attributes: Optional[Iterable[Dict[str, Any]]] = None, + *, + facets: Optional[Dict[str, Any]] = None, + provenance: Optional[Dict[str, Any]] = None, ) -> Dict[str, Any]: - """The value of a space at an instant, as one grid: ``{topology, geometry, attributes}``.""" - return {"topology": topo, "geometry": geom, - "attributes": list(attributes) if attributes is not None else []} + """The value of a space at an instant, as one grid: ``{topology, geometry, attributes}``, + and for a mesh optionally its ``facets`` and ``provenance`` (see :func:`grid`).""" + space = {"topology": topo, "geometry": geom, + "attributes": list(attributes) if attributes is not None else []} + space.update(_mesh_keys(facets, provenance)) + return space def spatial_grid(grids: Iterable[Dict[str, Any]]) -> Dict[str, Any]: @@ -894,6 +1058,55 @@ def _canonicalize_geometry(geom: Dict[str, Any], where: str) -> Dict[str, Any]: } +def _canonicalize_facets(f: Dict[str, Any], where: str) -> Dict[str, Any]: + for key in ("topology", "tags"): + if key not in f: + raise YmfArchiveError(f"{where}: facets has no {key!r}") + if not isinstance(f["tags"], dict) or "data_item" not in f["tags"]: + raise YmfArchiveError(f"{where}.tags: must be {{data_item: ...}}") + return { + "topology": _canonicalize_topology(f["topology"], f"{where}.topology"), + "tags": {"data_item": _canonicalize_data_item( + f["tags"]["data_item"], f"{where}.tags.data_item")}, + "view": str(f.get("view", "global")), + } + + +def _canonicalize_provenance(p: Dict[str, Any], where: str) -> Dict[str, Any]: + if not isinstance(p, dict) or "tool" not in p: + raise YmfArchiveError(f"{where}: provenance has no 'tool'") + unknown = sorted(set(p) - {"tool", "version", "files", "extra"}) + if unknown: + raise YmfArchiveError( + f"{where}: provenance takes tool, version, files and extra, not " + f"{', '.join(unknown)}; put anything else under extra") + canonical: Dict[str, Any] = {"tool": str(p["tool"])} + if p.get("version") is not None: + canonical["version"] = str(p["version"]) + files = [] + for i, entry in enumerate(p.get("files") or []): + if not isinstance(entry, dict) or set(entry) != {"name", "sha256"}: + raise YmfArchiveError( + f"{where}.files[{i}]: a file is {{name, sha256}}, got {entry!r}") + files.append({"name": str(entry["name"]), "sha256": str(entry["sha256"]).lower()}) + canonical["files"] = files + if p.get("extra"): + if not isinstance(p["extra"], dict): + raise YmfArchiveError(f"{where}.extra: must be a mapping") + canonical["extra"] = dict(p["extra"]) + return canonical + + +def _canonicalize_mesh_keys(src: Dict[str, Any], canonical: Dict[str, Any], + where: str) -> None: + """Carry a grid's optional ``facets`` and ``provenance`` into its canonical form.""" + if "facets" in src: + canonical["facets"] = _canonicalize_facets(src["facets"], f"{where}.facets") + if "provenance" in src: + canonical["provenance"] = _canonicalize_provenance( + src["provenance"], f"{where}.provenance") + + def _canonicalize_grid(g: Dict[str, Any], where: str) -> Dict[str, Any]: canonical: Dict[str, Any] = {} if "name" in g: @@ -904,6 +1117,7 @@ def _canonicalize_grid(g: Dict[str, Any], where: str) -> Dict[str, Any]: _canonicalize_attribute(a, f"{where}.attributes[{i}]") for i, a in enumerate(g.get("attributes", [])) ] + _canonicalize_mesh_keys(g, canonical, where) return canonical @@ -939,6 +1153,7 @@ def _canonicalize_step(step: Dict[str, Any], where: str) -> Dict[str, Any]: _canonicalize_attribute(a, f"{where}.attributes[{i}]") for i, a in enumerate(step.get("attributes", [])) ] + _canonicalize_mesh_keys(step, canonical, where) return canonical @@ -1132,6 +1347,72 @@ def _validate_grid_like(g: Dict[str, Any], where: str) -> None: seen[key] = i _validate_data_item(attr["data_item"], f"{at}.data_item") + if "facets" in g: + _validate_facets(g, f"{where}.facets") + if "provenance" in g: + for i, entry in enumerate(g["provenance"]["files"]): + if not _SHA256.match(entry["sha256"]): + raise YmfArchiveError( + f"{where}.provenance.files[{i}].sha256: {entry['sha256']!r} is " + "not a SHA-256 digest (64 hexadecimal digits)") + + +def _validate_facets(g: Dict[str, Any], where: str) -> None: + """The facet list's structure, and its tags where they are inline. + + Tags held in an HDF5 file are not read here (validation stays O(1) per + field); :func:`ymf.mesh.check_facets` reads them. + """ + f = g["facets"] + cell_type = g["topology"]["type"] + facet_topo = f["topology"] + allowed = _FACET_TYPES.get(cell_type) + if allowed is None and cell_type != "Mixed": + raise YmfArchiveError( + f"{where}: a {cell_type} grid cannot carry a facet list; facets are " + "for unstructured grids") + if allowed is not None and facet_topo["type"] not in allowed: + raise YmfArchiveError( + f"{where}.topology.type: the facets of {cell_type} cells are " + f"{' or '.join(sorted(allowed))}, not {facet_topo['type']!r}") + if facet_topo["type"] not in KNOWN_TOPOLOGY_TYPES: + raise YmfArchiveError( + f"{where}.topology.type: {facet_topo['type']!r} is not a known topology type") + _validate_data_item(facet_topo["data_item"], f"{where}.topology.data_item") + if facet_topo["data_item"]["data_type"] not in _INTEGER_DATA_TYPES: + raise YmfArchiveError(f"{where}.topology.data_item: node numbers must be integers") + if f["view"] not in FACET_VIEWS: + raise YmfArchiveError( + f"{where}.view: {f['view']!r} is not one of {sorted(FACET_VIEWS)}") + tags = f["tags"]["data_item"] + _validate_data_item(tags, f"{where}.tags.data_item") + if tags["data_type"] not in ("Int", "UInt"): + raise YmfArchiveError( + f"{where}.tags.data_item.data_type: tags are integers (Int), not " + f"{tags['data_type']!r}") + n = facet_topo["number_of_elements"] + dims = list(tags["dimensions"]) + if dims not in ([n], [n, 1]): + raise YmfArchiveError( + f"{where}.tags.data_item.dimensions: one tag per facet, {[n]}, got {dims}") + if "values" in tags: + check_facet_tags(tags["values"], f["view"], f"{where}.tags") + + +def check_facet_tags(tags: Iterable[int], view: str, where: str = "facets.tags") -> None: + """Refuse facet tags the view cannot hold: below -1 anywhere, and -1 in a + global mesh (a rank interface exists only once the mesh is partitioned).""" + lowest = min((int(t) for t in tags), default=0) + if lowest < -1: + raise YmfArchiveError( + f"{where}: tag {lowest} is not a tag; 0 is untagged, >0 a condition, " + "-1 a rank interface") + if lowest == -1 and view == "global": + raise YmfArchiveError( + f"{where}: tag -1 marks a rank interface, which a global mesh never " + "stores; the partitioner (or a per-rank reader) sets it, in a " + "subdomain view") + def validate_approximation(domain: Dict[str, Any]) -> Dict[str, Any]: """Check ``domain`` is structurally sound; return it canonicalized. @@ -1506,6 +1787,9 @@ def _data_items(domain: Dict[str, Any]): yield where + "/Geometry", g["geometry"] for a in g["attributes"]: yield "%s/%s/%s" % (where, a.get("center", "Node"), a["name"]), a + if "facets" in g: + yield where + "/Facets/Topology", g["facets"]["topology"] + yield where + "/Facets/Tags", g["facets"]["tags"] def step_indices(spec: Any, n: int) -> List[int]: diff --git a/src/ymf/closure.py b/src/ymf/closure.py index f4a5aac..dfb14ab 100644 --- a/src/ymf/closure.py +++ b/src/ymf/closure.py @@ -83,8 +83,9 @@ def _set(entry: Dict[str, Any], path: Tuple[str, ...], value: Any) -> None: entry[path[-1]] = value -def realizations(discretization: Dict[str, Any]) -> List[Dict[str, Any]]: - """The realizations a discretization entry asks for, in study order. +def realizations(discretization: Dict[str, Any], domain: Any = None) -> List[Dict[str, Any]]: + """The realizations a discretization entry asks for, in study order + (``domain``: the strong form's, which decides what a mesh needs). Each is a mapping such as ``{"cells": 8, "levels": 1}`` or ``{"h": 0.1, "levels": 1}`` (and ``"dt"`` when the entry steps in time). @@ -99,7 +100,7 @@ def realizations(discretization: Dict[str, Any]) -> List[Dict[str, Any]]: "cells is for a box, h for any other region" % (discretization.get("name"),)) method = method_of(discretization) - if method.missing(discretization): + if method.missing(discretization, domain): return [] # see ymf.discretizations.missing block = _get(discretization, (method.block,)) axes = [] @@ -219,10 +220,15 @@ def output_key(doc: Dict[str, Any], discretization: str, def planned_outputs(doc: Dict[str, Any]) -> Iterator[Tuple[str, str, Dict[str, Any]]]: """Every (key, discretization name, realization) the spec asks for.""" for entry in doc["solution_paths"]["discretizations"]: - for realization in realizations(entry): + for realization in realizations(entry, _domain_of(doc)): yield output_key(doc, entry["name"], realization), entry["name"], realization +def _domain_of(doc: Dict[str, Any]) -> Any: + """The strong form's domain, which decides what a mesh needs.""" + return ((doc.get("model") or {}).get("strong_form") or {}).get("domain") + + def _base_key(key: str) -> str: """A key without its ``~N`` suffix (a kept, differing rerun).""" return key.split("~", 1)[0] diff --git a/src/ymf/discretizations.py b/src/ymf/discretizations.py index 84e0ff2..b1592e3 100644 --- a/src/ymf/discretizations.py +++ b/src/ymf/discretizations.py @@ -29,7 +29,8 @@ from dataclasses import dataclass from typing import Any, Callable, Dict, List, Tuple -__all__ = ["Missing", "Method", "METHODS", "DEFAULT", "method_of", "missing"] +__all__ = ["Missing", "Method", "METHODS", "DEFAULT", "method_of", "missing", + "subproblem_solver"] @dataclass(frozen=True) @@ -62,16 +63,22 @@ class Method: studied: Tuple[Tuple[str, Tuple[str, ...]], ...] #: values an axis takes when the block is present but the axis is not defaults: Dict[str, Any] - #: what an entry of this method lacks (block-level checks are shared) - lacks: Callable[[Dict[str, Any]], List[Missing]] + #: what an entry of this method lacks (block-level checks are shared), + #: given the entry and the domain (or None) + lacks: Callable[..., List[Missing]] #: False: the fields validate, but no output is planned yet implemented: bool = True def axes(self) -> Tuple[Tuple[str, Tuple[str, ...]], ...]: return self.studied + (_DT,) - def missing(self, entry: Dict[str, Any]) -> List[Missing]: - """What ``entry`` still lacks before it can plan an output.""" + def missing(self, entry: Dict[str, Any], domain: Any = None) -> List[Missing]: + """What ``entry`` still lacks before it can plan an output. + + ``domain`` is the strong form's domain, where it decides what a + mesh needs: a domain that is already a mesh, or a file of one, fixes + it (see :func:`mesh_from_domain`). Without it, a mesh needs a size. + """ out = [] for other in METHODS.values(): if other is not self and other.block in entry: @@ -79,7 +86,7 @@ def missing(self, entry: Dict[str, Any]) -> List[Missing]: "a %s block that matches its method" % self.block, "%s: belongs to method %s, and this entry's method is %s" % (other.block, other.name, self.name))) - out += self.lacks(entry) + out += self.lacks(entry, domain) if not out and not self.implemented: out.append(Missing( "an implementation", @@ -89,7 +96,11 @@ def missing(self, entry: Dict[str, Any]) -> List[Missing]: def key_parts(self, realization: Dict[str, Any]) -> List[str]: """How a realization reads in an output key, before the digest.""" - if self.name == "mesh": + if self.name == "mesh" and not ("cells" in realization or "h" in realization): + # a mesh the domain fixes (a file of one, or a discrete domain): + # only its levels of refinement are the discretization's + parts = ["levels=%s" % (realization.get("levels", 1),)] + elif self.name == "mesh": # the key format archives were first written with: cells or h, # then levels only when it is not 1 size = "cells" if "cells" in realization else "h" @@ -104,8 +115,30 @@ def key_parts(self, realization: Dict[str, Any]) -> List[str]: return parts -def _mesh_lacks(entry: Dict[str, Any]) -> List[Missing]: +#: file formats whose file is a mesh; the others (gmsh_geo) are meshed by h +MESH_FILE_FORMATS = ("2dm", "3dm", "hex", "tetgen", "ymf") + + +def mesh_from_domain(domain: Any) -> bool: + """Whether the domain fixes the mesh: a ``discrete`` domain, or a + ``file`` domain whose file is a mesh rather than a geometry.""" + if not isinstance(domain, dict): + return False + return domain.get("kind") == "discrete" or ( + domain.get("kind") == "file" and domain.get("format") in MESH_FILE_FORMATS) + + +def _mesh_lacks(entry: Dict[str, Any], domain: Any = None) -> List[Missing]: mesh = entry.get("mesh") + if mesh_from_domain(domain): + if not isinstance(mesh, dict): + return [Missing("a mesh", "give mesh: {levels: 1}: the domain is the mesh, " + "and levels refines it uniformly")] + if "cells" in mesh or "h" in mesh: + return [Missing("a mesh the domain does not fix", + "this %s domain is already a mesh; give mesh: {levels: ...} " + "only, not cells or h" % (domain.get("kind"),))] + return [] if not isinstance(mesh, dict) or not ("cells" in mesh or "h" in mesh): return [Missing("a mesh", "give mesh: {h: ...} (or {cells: ...} on a box), " "or --h (--cells) on the command line")] @@ -116,7 +149,7 @@ def _mesh_lacks(entry: Dict[str, Any]) -> List[Missing]: def _fields_lacking(block: str, required: Tuple[str, ...], meaning: Dict[str, str]): - def lacks(entry: Dict[str, Any]) -> List[Missing]: + def lacks(entry: Dict[str, Any], domain: Any = None) -> List[Missing]: given = entry.get(block) given = given if isinstance(given, dict) else {} return [Missing("%s.%s" % (block, f), "give %s: {%s: ...}: %s" % (block, f, meaning[f])) @@ -165,6 +198,25 @@ def method_of(entry: Dict[str, Any]) -> Method: return METHODS[name] -def missing(entry: Dict[str, Any]) -> List[Missing]: - """What ``entry`` lacks before it can plan an output (empty: nothing).""" - return method_of(entry).missing(entry) +def missing(entry: Dict[str, Any], domain: Any = None) -> List[Missing]: + """What ``entry`` lacks before it can plan an output (empty: nothing); + ``domain`` as in :meth:`Method.missing`.""" + return method_of(entry).missing(entry, domain) + + +def subproblem_solver(entry: Dict[str, Any], subproblem: Dict[str, Any]) -> Dict[str, Any]: + """A sub-problem's solver: the entry's ``solver``, refined by its own. + + The entry's solver is the default every sub-problem starts from; a key a + sub-problem's ``solver:`` gives replaces the entry's for that sub-problem + alone (``max_iterations: 50`` for the transport solves, say), and + ``petsc_options`` merge key by key. + """ + base = dict(entry.get("solver") or {}) + own = dict(subproblem.get("solver") or {}) + options = dict(base.get("petsc_options") or {}) + options.update(own.pop("petsc_options", None) or {}) + base.update(own) + if options: + base["petsc_options"] = options + return base diff --git a/src/ymf/domain.py b/src/ymf/domain.py index e76863f..50a9d5c 100644 --- a/src/ymf/domain.py +++ b/src/ymf/domain.py @@ -66,11 +66,15 @@ class YmfDomainError(ValueError): }, "discrete": { "required": {"topology", "geometry"}, - "optional": set(), + # a mesh as ymf stores one (ymf.mesh, C7.1): its tagged facets, its + # provenance, its attributes; and bounds, for nodes held in a file + "optional": {"facets", "provenance", "attributes", "bounds"}, }, "file": { "required": {"format", "path"}, - "optional": set(), + # the bounding box, which sets the scale of the boundary tolerance; + # without it a solver takes it from the mesh (posing follows meshing) + "optional": {"bounds"}, }, } @@ -135,6 +139,21 @@ def validate_domain_block(domain: Any, where: str = "domain") -> None: raise YmfDomainError( "%s: dimension %r but %d components" % (where, dim, len(lower))) + bounds = domain.get("bounds") + if bounds is not None: + lower, upper = bounds.get("lower") or [], bounds.get("upper") or [] + if not lower or len(lower) != len(upper): + raise YmfDomainError("%s.bounds: lower and upper need the same, nonzero, number " + "of components" % (where,)) + for i, (a, b) in enumerate(zip(lower, upper)): + if not b > a: + raise YmfDomainError("%s.bounds: upper[%d] = %r is not above lower[%d] = %r" + % (where, i, b, i, a)) + dim = domain.get("dimension") + if dim is not None and dim != len(lower): + raise YmfDomainError("%s: dimension %r but bounds of %d components" + % (where, dim, len(lower))) + time = domain.get("time") if time is not None: if len(time) != 2 or not time[1] > time[0]: diff --git a/src/ymf/mesh.py b/src/ymf/mesh.py new file mode 100644 index 0000000..8ea51a2 --- /dev/null +++ b/src/ymf/mesh.py @@ -0,0 +1,201 @@ +"""Helpers for writing and checking a mesh stored as a ymf grid. + +A mesh is a grid (``ymf.archive``): ``topology``, ``geometry`` and +``attributes``, plus two optional keys that only a mesh needs -- its tagged +``facets`` and its ``provenance``. The schema, and the cheap structural +checks, live in ``ymf.archive``, which depends on pyyaml alone; the helpers +here read files and arrays, so they import numpy (and h5py, to read a +referenced array) when called. + +Nothing here names an external mesh format. A converter from some code's +input files to ymf states that code's conventions as data -- for instance +the vertex pairs its quadratic cells' edge nodes sit on -- and these +helpers do the rest. +""" + +from __future__ import annotations + +import hashlib +from pathlib import Path +from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple + +from .archive import ( + QUADRATIC_EDGES, + YmfArchiveError, + check_facet_tags, + load_array, + provenance, +) + +__all__ = [ + "file_sha256", + "file_provenance", + "edge_node_permutation", + "reorder_cells", + "edge_midpoint_offsets", + "check_facets", +] + +#: The faces of each linear 3D cell and the edges of each linear 2D cell, as +#: tuples of the cell's vertex numbers (VTK's order, as for the node order +#: in :data:`ymf.archive.QUADRATIC_EDGES`). A facet is matched against them +#: as a set of vertices, so the order within a face does not matter here. +_CELL_FACETS = { + "Polyline": ((0,), (1,)), + "Triangle": ((0, 1), (1, 2), (2, 0)), + "Quadrilateral": ((0, 1), (1, 2), (2, 3), (3, 0)), + "Tetrahedron": ((0, 1, 3), (1, 2, 3), (2, 0, 3), (0, 2, 1)), + "Hexahedron": ((0, 4, 7, 3), (1, 2, 6, 5), (0, 1, 5, 4), (3, 7, 6, 2), + (0, 3, 2, 1), (4, 5, 6, 7)), + "Wedge": ((0, 1, 2), (3, 5, 4), (0, 3, 4, 1), (1, 4, 5, 2), (2, 5, 3, 0)), + "Pyramid": ((0, 3, 2, 1), (0, 1, 4), (1, 2, 4), (2, 3, 4), (3, 0, 4)), +} +_LINEAR = { + "Edge_3": "Polyline", "Tri_6": "Triangle", "Triangle_6": "Triangle", + "Quad_8": "Quadrilateral", "Quadrilateral_8": "Quadrilateral", + "Tet_10": "Tetrahedron", "Tetrahedron_10": "Tetrahedron", + "Pyramid_13": "Pyramid", "Wedge_15": "Wedge", + "Hex_20": "Hexahedron", "Hexahedron_20": "Hexahedron", +} +#: Vertices of a facet of each type (the leading nodes of a quadratic facet). +_FACET_VERTICES = {"Polyvertex": 1, "Polyline": 2, "Edge_3": 2, "Triangle": 3, + "Tri_6": 3, "Triangle_6": 3, "Quadrilateral": 4, "Quad_8": 4, + "Quadrilateral_8": 4} + + +def file_sha256(path: str | Path) -> str: + """The SHA-256 digest of a file's bytes, as 64 lowercase hex digits.""" + digest = hashlib.sha256() + with open(path, "rb") as f: + for block in iter(lambda: f.read(1 << 20), b""): + digest.update(block) + return digest.hexdigest() + + +def file_provenance( + tool: str, + paths: Iterable[str | Path], + *, + version: Optional[str] = None, + extra: Optional[Dict[str, Any]] = None, +) -> Dict[str, Any]: + """A :func:`ymf.archive.provenance` block for a mesh made from ``paths``. + + Each file is recorded by its base name and the SHA-256 of its bytes, so + a reader can tell whether the mesh still matches the files beside it + without trusting their timestamps. + """ + files = [{"name": Path(p).name, "sha256": file_sha256(p)} for p in paths] + return provenance(tool, version=version, files=files, extra=extra) + + +def edge_node_permutation(topology_type: str, + source_edges: Sequence[Tuple[int, int]]) -> List[int]: + """The column order that takes a source's quadratic cells to ymf's. + + ``source_edges`` lists, for the source's edge nodes in its own order, + the pair of vertices each sits on. The result ``perm`` gives + ``ymf_cells = source_cells[:, perm]``: vertices unchanged, edge nodes + reordered to :data:`ymf.archive.QUADRATIC_EDGES`. An edge is the same + whichever way round its pair is written. + """ + if topology_type not in QUADRATIC_EDGES: + raise YmfArchiveError( + f"{topology_type!r} is not a quadratic cell type; quadratic types are " + f"{sorted(QUADRATIC_EDGES)}") + target = QUADRATIC_EDGES[topology_type] + n_vertices = max(max(e) for e in target) + 1 + source = [frozenset(e) for e in source_edges] + if len(source) != len(target) or len(set(source)) != len(source): + raise YmfArchiveError( + f"{topology_type} has {len(target)} distinct edges; the source lists " + f"{list(source_edges)}") + perm = list(range(n_vertices)) + for edge in target: + try: + perm.append(n_vertices + source.index(frozenset(edge))) + except ValueError: + raise YmfArchiveError( + f"{topology_type}: edge {edge} is not among the source's edges " + f"{list(source_edges)}") from None + return perm + + +def reorder_cells(cells, topology_type: str, source_edges: Sequence[Tuple[int, int]]): + """``cells`` (one row per cell) with the edge nodes put in ymf's order.""" + import numpy + return numpy.asarray(cells)[:, edge_node_permutation(topology_type, source_edges)] + + +def edge_midpoint_offsets(topology_type: str, nodes, cells): + """How far each cell's edge nodes sit from the midpoints of their edges. + + One number per cell: the largest distance from an edge node to the + midpoint of the edge :data:`ymf.archive.QUADRATIC_EDGES` puts it on, + relative to that edge's length. For a straight-sided mesh it is zero up + to rounding when the nodes are in ymf's order, and of order one when + two edge nodes are swapped -- which makes it the check a converter runs + on a straight-sided sample of its source. A curved mesh legitimately + gives nonzero values. + """ + import numpy + if topology_type not in QUADRATIC_EDGES: + raise YmfArchiveError(f"{topology_type!r} is not a quadratic cell type") + edges = numpy.array(QUADRATIC_EDGES[topology_type]) + nodes = numpy.asarray(nodes, dtype=float) + cells = numpy.asarray(cells) + n_vertices = int(edges.max()) + 1 + a = nodes[cells[:, edges[:, 0]]] + b = nodes[cells[:, edges[:, 1]]] + mid = nodes[cells[:, n_vertices:n_vertices + len(edges)]] + length = numpy.linalg.norm(b - a, axis=-1) + offset = numpy.linalg.norm(mid - 0.5 * (a + b), axis=-1) / numpy.where(length > 0, length, 1.0) + return offset.max(axis=1) + + +def check_facets(grid: Dict[str, Any], base_dir: str | Path = ".") -> Dict[str, int]: + """Read a grid's facet list and check what structure alone cannot. + + The tags are legal for the view (:func:`ymf.archive.check_facet_tags`); + every node number is a node of the grid; and every facet is a face (an + edge, in 2D) of some cell, matched on its vertices. Returns + ``{tag: count}``. ``grid`` is a canonical grid or uniform space, as + :func:`ymf.archive.validate_approximation` returns it. + """ + import numpy + if "facets" not in grid: + raise YmfArchiveError("grid has no facets") + f = grid["facets"] + where = "facets" + tags = numpy.asarray(load_array(f["tags"]["data_item"], base_dir)).reshape(-1) + check_facet_tags(tags.tolist(), f["view"], f"{where}.tags") + facet_type = f["topology"]["type"] + facet_nodes = numpy.asarray(load_array(f["topology"]["data_item"], base_dir)) + n_facets = f["topology"]["number_of_elements"] + facet_nodes = facet_nodes.reshape(n_facets, -1) + if len(tags) != n_facets: + raise YmfArchiveError(f"{where}: {n_facets} facets but {len(tags)} tags") + geometry = numpy.asarray(load_array(grid["geometry"]["data_item"], base_dir)) + n_nodes = geometry.shape[0] + if facet_nodes.size and (facet_nodes.min() < 0 or facet_nodes.max() >= n_nodes): + raise YmfArchiveError( + f"{where}.topology: node numbers must lie in [0, {n_nodes}); got " + f"{int(facet_nodes.min())}..{int(facet_nodes.max())}") + cell_type = grid["topology"]["type"] + linear = _LINEAR.get(cell_type, cell_type) + faces = _CELL_FACETS.get(linear) + if faces is not None and facet_type != "Mixed": + cells = numpy.asarray(load_array(grid["topology"]["data_item"], base_dir)) + cells = cells.reshape(grid["topology"]["number_of_elements"], -1) + known = set() + for face in faces: + for row in cells[:, list(face)].tolist(): + known.add(frozenset(row)) + k = _FACET_VERTICES[facet_type] + for i, row in enumerate(facet_nodes[:, :k].tolist()): + if frozenset(row) not in known: + raise YmfArchiveError( + f"{where}: facet {i} (nodes {row}) is not a face of any " + f"{cell_type} cell of the grid") + values, counts = numpy.unique(tags, return_counts=True) + return {int(v): int(c) for v, c in zip(values, counts)} diff --git a/src/ymf/normalize.py b/src/ymf/normalize.py index 9add8b0..3b83cda 100644 --- a/src/ymf/normalize.py +++ b/src/ymf/normalize.py @@ -51,3 +51,118 @@ def normalize_unknowns( f"Unsupported unknowns entry type: {type(entry)!r} (value={entry!r})" ) return result + + +EquationEntry = Union[str, Dict[str, Any]] + + +def normalize_equations(equations: List[EquationEntry]) -> List[Dict[str, Any]]: + """Normalize ``strong_form.equations`` to ``{name, for, equation}`` entries. + + An equation is written either as a bare string, ordered like the + unknowns, or as a mapping that names it and, with ``for``, the unknown + it determines:: + + - name: momentum + for: v + equation: "ρ ∂v/∂t + ∇·(ρ v⊗v) − ∇·(μ∇v) + ∇p = ρ g in Ω" + + A name is what a discretization's sub-problem refers to (``solves: + [momentum]``), and what an extending file merges onto. ``name`` and + ``for`` are ``None`` for a bare string. + + >>> normalize_equations(["-Δu = f in Ω"]) + [{'name': None, 'for': None, 'equation': '-Δu = f in Ω'}] + """ + result: List[Dict[str, Any]] = [] + for entry in equations: + if isinstance(entry, str): + result.append({"name": None, "for": None, "equation": entry}) + elif isinstance(entry, dict): + result.append({"name": entry.get("name"), "for": entry.get("for"), + "equation": entry.get("equation")}) + else: + raise TypeError( + f"Unsupported equations entry type: {type(entry)!r} (value={entry!r})" + ) + return result + + +def equations_in_order(strong_form: Dict[str, Any], complete: bool = True) -> List[str]: + """The equations' text, one per unknown, in the order of the unknowns. + + Bare strings are taken in the order written, as they always were. When + every equation says which unknown it is ``for``, the order written does + not matter: they are put in the unknowns' order. Raises ``ValueError`` + when some say and some do not, when one names something that is not an + unknown, or when two are for the same unknown -- each is a document + that does not say which equation determines which unknown. + ``complete=False`` allows unknowns no equation is for yet (an abstract + file the rest of the chain completes); the result then holds only the + equations given. + """ + entries = normalize_equations(strong_form.get("equations") or []) + stated = [e for e in entries if e["for"] is not None] + if not stated: + return [e["equation"] for e in entries] + if len(stated) != len(entries): + raise ValueError( + "equations: give 'for' on every equation or on none; %s %s but %s %s not" + % (_listing([e["name"] for e in stated]), "says" if len(stated) == 1 else "say", + _listing([e["name"] or e["equation"] for e in entries if e["for"] is None]), + "does" if len(entries) - len(stated) == 1 else "do")) + names = [u["name"] for u in normalize_unknowns(strong_form.get("unknowns") or [])] + by_unknown: Dict[str, Dict[str, Any]] = {} + for e in entries: + if e["for"] not in names: + raise ValueError("equations[%s] is for %r, which is not an unknown (have: %s)" + % (e["name"], e["for"], ", ".join(names))) + if e["for"] in by_unknown: + raise ValueError("equations %s and %s are both for %s" + % (by_unknown[e["for"]]["name"], e["name"], e["for"])) + by_unknown[e["for"]] = e + missing = [n for n in names if n not in by_unknown] + if missing and complete: + raise ValueError("equations: none is for %s" % ", ".join(missing)) + return [by_unknown[n]["equation"] for n in names if n in by_unknown] + + +def _listing(items: List[Any]) -> str: + return ", ".join(str(i) for i in items) + + +def boundary_conditions(strong_form: Dict[str, Any]) -> List[Dict[str, Any]]: + """Every boundary condition, with the condition sets expanded. + + A model names sets of conditions without a region (``condition_sets: + {free_slip: [{variable: v, type: slip}]}``) and a problem applies them + to regions (``boundary_sets: [{region: walls, set: free_slip}]``). The + result is ``boundary_conditions`` followed by one condition per member + of each applied set, with the region filled in and ``set`` naming where + it came from. A condition written out for the same region and variable + as a set's member replaces it: a case can apply ``atmosphere`` and still + give the pressure its own value. Raises ``ValueError`` for a set nobody + defines, and for two sets that both give one region's variable a + condition, which nothing would choose between. + """ + sets = strong_form.get("condition_sets") or {} + direct = [dict(bc) for bc in strong_form.get("boundary_conditions") or []] + written = {(bc.get("region"), bc.get("variable")) for bc in direct} + expanded: List[Dict[str, Any]] = [] + given: Dict[Any, str] = {} + for applied in strong_form.get("boundary_sets") or []: + name, region = applied.get("set"), applied.get("region") + if name not in sets: + raise ValueError("boundary_sets[region=%r]: no condition set named %r (have: %s)" + % (region, name, ", ".join(sorted(sets)) or "none")) + for member in sets[name]: + key = (region, member.get("variable")) + if key in written: + continue + if key in given and given[key] != name: + raise ValueError( + "boundary_sets: sets %s and %s both give %s a condition on %s" + % (given[key], name, member.get("variable"), region)) + given[key] = name + expanded.append(dict(member, region=region, set=name)) + return direct + expanded diff --git a/src/ymf/references.py b/src/ymf/references.py index e98bffd..3a37f9f 100644 --- a/src/ymf/references.py +++ b/src/ymf/references.py @@ -137,12 +137,289 @@ def _check(named: Any, candidates: List[str], path: str, what: str) -> None: "%s must name a %s; found %r" % (path, what, named)) if named not in candidates: raise YmfReferenceError( - "%s names %r, which is not a %s in this document (have: %s)" - % (path, named, what, ", ".join(candidates) or "none")) + "%s names %r, which is not %s %s in this document (have: %s)" + % (path, named, "an" if what[:1] in "aeiou" else "a", what, + ", ".join(candidates) or "none")) + + +def _check_equations(doc: Dict[str, Any]) -> None: + """Named equations: each name once, and ``for`` naming the unknowns.""" + from ymf.normalize import equations_in_order, normalize_equations + strong = _problem(doc).get("strong_form") or {} + names = [e["name"] for e in normalize_equations(strong.get("equations") or []) + if e["name"] is not None] + twice = sorted({n for n in names if names.count(n) > 1}) + if twice: + raise YmfReferenceError("model.strong_form.equations: %s named more than once" + % ", ".join(twice)) + if strong.get("equations") and strong.get("unknowns"): + try: + equations_in_order(strong, complete=doc.get("ymf_state") != "abstract") + except ValueError as exc: + raise YmfReferenceError("model.strong_form.%s" % exc) from exc + + +#: sharp functions of the notation a discretization may regularize +SHARP_FUNCTIONS = ("heaviside", "sign") + + +def _check_split(doc: Dict[str, Any], entry: Dict[str, Any]) -> None: + """A discretization's own names: its sub-problems, bindings and + regularizations, against the model and each other.""" + from ymf.normalize import normalize_equations, normalize_unknowns + name = entry.get("name", "?") + where = "solution_paths.discretizations[%s]" % (name,) + subproblems = [s for s in entry.get("subproblems") or [] if isinstance(s, dict)] + strong = _problem(doc).get("strong_form") or {} + unknowns = [u["name"] for u in normalize_unknowns(strong.get("unknowns") or [])] + coefficients = list(strong.get("coefficients") or {}) + model_names = set(unknowns) | set(coefficients) + for symbol in entry.get("binds") or {}: + if symbol in model_names: + raise YmfReferenceError( + "%s.binds.%s: %s is already a name in the model; a discretization " + "binds symbols the model leaves to it" % (where, symbol, symbol)) + for function in entry.get("regularize") or {}: + if function not in SHARP_FUNCTIONS: + raise YmfReferenceError( + "%s.regularize.%s: only the sharp functions %s are regularized" + % (where, function, ", ".join(SHARP_FUNCTIONS))) + if not subproblems: + return + names = [s["name"] for s in subproblems] + twice = sorted({n for n in names if names.count(n) > 1}) + if twice: + raise YmfReferenceError("%s.subproblems: %s named more than once" + % (where, ", ".join(twice))) + equations = [e["name"] for e in normalize_equations(strong.get("equations") or [])] + introduced: Dict[str, str] = {} + for s in subproblems: + new = (s.get("introduces") or {}).get("name") + if new is None: + continue + if new in model_names or new in introduced: + raise YmfReferenceError( + "%s.subproblems[%s].introduces: %s is already a name (%s)" + % (where, s["name"], new, "in the model" if new in model_names + else "introduced by %s" % introduced[new])) + introduced[new] = s["name"] + readable = model_names | set(introduced) + solved_by: Dict[str, str] = {} + for s in subproblems: + here = "%s.subproblems[%s]" % (where, s["name"]) + if not (s.get("solves") or s.get("enforces") or s.get("introduces") + or s.get("spin_up")): + raise YmfReferenceError( + "%s neither solves a model equation, enforces a constraint, " + "introduces an unknown nor spins up: give one" % (here,)) + for equation in s.get("solves") or []: + if None in equations: + raise YmfReferenceError( + "%s.solves names %r, but the model's equations are not named; " + "name them ({name, for, equation}) to split them" % (here, equation)) + _check(equation, equations, here + ".solves", "model equation") + if equation in solved_by: + raise YmfReferenceError( + "equation %s is solved by both %s and %s in %s" + % (equation, solved_by[equation], s["name"], where)) + solved_by[equation] = s["name"] + new = s.get("introduces") or {} + for key, candidates in (("approximates", model_names), ("correction_of", unknowns)): + if key in new: + _check(new[key], sorted(candidates), "%s.introduces.%s" % (here, key), + "model unknown" if key == "correction_of" else + "model unknown or coefficient") + for field, lag in (s.get("lags") or {}).items(): + copy_of = lag.get("of") if isinstance(lag, dict) else None + if copy_of is not None: + _check(copy_of, sorted(readable), "%s.lags.%s.of" % (here, field), + "model unknown, coefficient or introduced unknown") + if field in readable: + raise YmfReferenceError( + "%s.lags.%s: a copy (of: %s) takes a name of its own, and %s is " + "already one" % (here, field, copy_of, field)) + else: + _check(field, sorted(readable), "%s.lags" % (here,), + "model unknown, coefficient or introduced unknown") + if isinstance(lag, dict) and "via" in lag: + _check(lag["via"], names, "%s.lags.%s.via" % (here, field), "sub-problem") + if "for" in s: + _check(s["for"], sorted(set(unknowns) | set(introduced)), here + ".for", + "model or introduced unknown") + for function in s.get("regularize") or {}: + if function not in SHARP_FUNCTIONS: + raise YmfReferenceError( + "%s.regularize.%s: only the sharp functions %s are regularized" + % (here, function, ", ".join(SHARP_FUNCTIONS))) + for field in s.get("writes_back") or {}: + _check(field, sorted(set(unknowns) | set(introduced)), + "%s.writes_back" % (here,), "model or introduced unknown") + for symbol in s.get("binds") or {}: + if symbol in readable: + raise YmfReferenceError( + "%s.binds.%s: %s is already a name" % (here, symbol, symbol)) + if None not in equations: + unsolved = [e for e in equations if e not in solved_by] + if unsolved: + raise YmfReferenceError( + "%s: no sub-problem solves %s" % (where, ", ".join(unsolved))) + + +#: what each kind of condition needs: a value (value, formula or series), +#: none, or either +_NEEDS_VALUE = {"dirichlet": True, "neumann": True, "robin": True, "inflow": True, + "characteristic": True, "periodic": False, "open": False} + + +def _check_conditions(doc: Dict[str, Any]) -> None: + """Boundary conditions: their region, their variable, and what their kind needs.""" + strong = _problem(doc).get("strong_form") or {} + unknowns = {u["name"]: int(u.get("rank") or 0) + for u in _unknown_entries(strong.get("unknowns") or [])} + coefficients = set(strong.get("coefficients") or {}) + regions = {r.get("name") for r in strong.get("boundary_regions") or [] + if isinstance(r, dict)} + complete = doc.get("ymf_state") != "abstract" + from ymf.normalize import boundary_conditions + try: + conditions = boundary_conditions(strong) + except ValueError as exc: + raise YmfReferenceError("model.strong_form.%s" % exc) from exc + members = [dict(m, set=name) for name, entries in (strong.get("condition_sets") or {}).items() + for m in entries or [] if isinstance(m, dict)] + for bc in conditions + members: + if not isinstance(bc, dict): + continue + where = "model.strong_form.%s[%s]" % ( + ("condition_sets.%s" % bc["set"]) if "region" not in bc else "boundary_conditions", + ("variable=%r" % bc.get("variable")) if "region" not in bc else + "region=%r, variable=%r%s" % (bc.get("region"), bc.get("variable"), + ", set %s" % bc["set"] if bc.get("set") else "")) + if "region" in bc and complete and regions and bc.get("region") not in regions: + raise YmfReferenceError("%s: no boundary region named %r (have: %s)" + % (where, bc.get("region"), ", ".join(sorted(regions)))) + variable, kind = bc.get("variable"), bc.get("type") + if unknowns and variable not in unknowns and variable not in coefficients: + raise YmfReferenceError("%s: %s is neither an unknown nor a coefficient" + % (where, variable)) + given = any(bc.get(k) is not None for k in ("value", "formula", "series")) + if kind == "slip" and variable in unknowns and unknowns[variable] != 1: + raise YmfReferenceError("%s: slip gives the normal component of a vector, and " + "%s is a scalar" % (where, variable)) + article = "an" if str(kind)[:1] in "aeiou" else "a" + if _NEEDS_VALUE.get(kind) and not given: + raise YmfReferenceError("%s: %s %s condition needs a value, formula or series" + % (where, article, kind)) + if _NEEDS_VALUE.get(kind) is False and given: + raise YmfReferenceError("%s: %s %s condition takes no value" + % (where, article, kind)) + + +def _unknown_entries(unknowns): + from ymf.normalize import normalize_unknowns + ranks = {u["name"]: u.get("rank") for u in unknowns if isinstance(u, dict)} + return [dict(u, rank=ranks.get(u["name"])) for u in normalize_unknowns(unknowns)] + + +#: gauge kind -> (keys it needs, keys it may not have) +_GAUGE_KEYS = { + "point": (("fields", "points"), ("lines", "tags", "quantities", "moment_point")), + "line": (("fields", "lines"), ("points", "tags", "quantities", "moment_point")), + "line_integral": (("fields", "lines"), ("points", "tags", "quantities", "moment_point")), + "boundary_integral": (("tags", "quantities"), ("fields", "points", "lines")), +} + + +def _check_gauges(doc: Dict[str, Any]) -> None: + """What each gauge's kind needs, and that its names exist.""" + model = _problem(doc) + strong = model.get("strong_form") or {} + names = ({u["name"] for u in _unknown_entries(strong.get("unknowns") or [])} + | set(strong.get("coefficients") or {})) + regions = {r.get("name") for r in strong.get("boundary_regions") or [] if isinstance(r, dict)} + complete = doc.get("ymf_state") != "abstract" + for i, gauge in enumerate(model.get("gauges") or []): + if not isinstance(gauge, dict): + continue + where = "model.gauges[%s]" % (gauge.get("name") or i,) + kind = gauge.get("kind") + needs, refuses = _GAUGE_KEYS.get(kind, ((), ())) + for key in needs: + if not gauge.get(key): + raise YmfReferenceError("%s: a %s gauge needs %s" % (where, kind, key)) + for key in refuses: + if key in gauge: + raise YmfReferenceError("%s: a %s gauge takes no %s" % (where, kind, key)) + if complete and names: + for field in gauge.get("fields") or []: + _check(field, sorted(names), where + ".fields", "unknown or coefficient") + if complete and regions: + for tag in gauge.get("tags") or []: + _check(tag, sorted(regions), where + ".tags", "boundary region") + points = list(gauge.get("points") or []) + for line in gauge.get("lines") or []: + if len(line) != 2: + raise YmfReferenceError("%s: a line is its two ends, not %d points" + % (where, len(line))) + points.extend(line) + if gauge.get("moment_point"): + points.append(gauge["moment_point"]) + if len({len(p) for p in points}) > 1: + raise YmfReferenceError("%s: its points do not all have the same dimension" + % (where,)) + + +def _check_regions(doc: Dict[str, Any]) -> None: + """A boundary region is named by ``where`` or by ``tags``, exactly one; a + tag is an int or a name the domain's ``boundary_tags`` resolves; a box's + faces are named by position; a periodic region is found by position.""" + strong = _problem(doc).get("strong_form") or {} + domain = strong.get("domain") + is_box = isinstance(domain, str) or (isinstance(domain, dict) + and domain.get("kind") == "box") + names = (domain.get("boundary_tags") if isinstance(domain, dict) else None) or {} + periodic = set() + try: + from ymf.normalize import boundary_conditions + periodic = {bc.get("region") for bc in boundary_conditions(strong) + if bc.get("type") == "periodic"} + except ValueError: + pass # an undefined set: _check_conditions says so + for region in strong.get("boundary_regions") or []: + if not isinstance(region, dict): + continue + where = "model.strong_form.boundary_regions[%s]" % (region.get("name"),) + given = [k for k in ("where", "tags") if region.get(k) is not None] + if len(given) != 1: + raise YmfReferenceError( + "%s: give where (a predicate) or tags (the mesh's boundary tags), " + "exactly one; this gives %s" % (where, " and ".join(given) or "neither")) + if "tags" not in given: + continue + if is_box and domain is not None: + raise YmfReferenceError( + "%s: a box carries no boundary tags; name its faces with where " + "(\"x = 0\")" % (where,)) + for tag in region["tags"]: + if isinstance(tag, str) and not tag.lstrip("-").isdigit() and tag not in names: + raise YmfReferenceError( + "%s: tag %r is not one of the domain's boundary_tags (have: %s)" + % (where, tag, ", ".join(sorted(names)) or "none")) + if region.get("name") in periodic: + raise YmfReferenceError( + "%s: a periodic condition finds its faces by position, so the region " + "needs where:, not tags:" % (where,)) def resolve(doc: Dict[str, Any]) -> None: """Check every reference, raising :class:`YmfReferenceError` on the first bad one.""" + _check_equations(doc) + _check_regions(doc) + _check_conditions(doc) + _check_gauges(doc) + for entry in _paths(doc).get("discretizations") or []: + if isinstance(entry, dict): + _check_split(doc, entry) labels = weak_form_labels(doc) discretizations = discretization_names(doc) diff --git a/src/ymf/schema.py b/src/ymf/schema.py index 5d8e9c4..5c7e0bb 100644 --- a/src/ymf/schema.py +++ b/src/ymf/schema.py @@ -48,7 +48,14 @@ # discharge hydrograph is a tabulated neumann condition, and a stage # hydrograph a tabulated dirichlet one, which naming it as a type made # impossible to write. -BCTypeEnum = Enum(["dirichlet", "neumann", "robin", "periodic"]) +# slip, open, inflow and characteristic (C42): the normal component of a +# vector given and its tangential traction zero; a declared open boundary +# (the natural flux, C12); a value imposed only where the flow enters; a +# state imposed through the incoming characteristics (a tailwater stage). +# They are for a split discretization's solver: one ADR system runs the +# first four only (ymf.symbolic.adr_problem). +BCTypeEnum = Enum(["dirichlet", "neumann", "robin", "periodic", + "slip", "open", "inflow", "characteristic"]) ICTypeEnum = Enum(["constant", "function", "data_source"]) @@ -121,6 +128,66 @@ } ) +# --- splitting a discretization into sub-problems (C42) --------------------- +# A split discretization solves the model as a sequence of sub-problems in +# one time step: two-phase level-set flow as flow, VOF, level set, +# redistancing and mass correction. Each names the model equations it +# solves, or the constraint it enforces, the unknowns it introduces, what it +# reads from the others and when, and how it is solved. ymf checks the +# structure and the names (ymf.references); a solver's dispatcher maps a +# sub-problem onto a module. + +#: when a sub-problem reads another's field: as the previous system step +#: left it, or as an earlier sub-problem of this step left it +LagAtEnum = Enum(["previous_step", "current_step"]) + +SymbolBinding = Map({ + Optional("meaning"): Str(), # e.g. element_size + Optional("formula"): Str(), # e.g. "1.5 h" + Optional("value"): Float(), + Optional("proteus"): Any(), # an interim hint, see SubproblemDef +}) + +SubproblemDef = Map({ + "name": Str(), + # model equations (by name) this sub-problem solves + Optional("solves"): Seq(Str()), + # a constraint it restores ("|∇φ| = 1", "heaviside(φ) = χh"), in the notation + Optional("enforces"): Str(), + # an unknown it adds: an auxiliary carrier of a model quantity + # (approximates), or a correction to a model unknown (correction_of) + Optional("introduces"): Map({ + "name": Str(), + Optional("approximates"): Str(), + Optional("correction_of"): Str(), + Optional("units"): Str(), + }), + # the equation of an introduced unknown, or of a constraint's solve, and + # the unknown it determines (default: the one this sub-problem introduces) + Optional("equation"): Str(), + Optional("for"): Str(), + Optional("initial_condition"): Str(), + # what it reads from elsewhere: name -> when, or {via: sub-problem, at: + # when}; with `of`, the name is a copy of that field under a name of its + # own (redistancing reads φs, the level set before it changes φ) + Optional("lags"): MapPattern(Str(), LagAtEnum | Map({ + Optional("via"): Str(), "at": LagAtEnum, Optional("of"): Str()})), + # what it hands back: name -> replace | add | mass_quadrature | ... + Optional("writes_back"): MapPattern(Str(), Str()), + Optional("pseudo_time"): Enum(["steady"]), + Optional("spin_up"): Bool(), + Optional("binds"): MapPattern(Str(), SymbolBinding), + # replaces the entry's regularize: for this sub-problem's sharp functions + Optional("regularize"): MapPattern(Str(), Map({"by": Str(), Optional("width"): Str()})), + Optional("finite_element"): FiniteElementDef, + Optional("time_integration"): Str(), + Optional("stabilization"): Any(), + Optional("solver"): SolverConfig, + # options a dispatcher cannot infer yet, for one solver: an interim hint + # it checks against what it infers, and drops once it covers them + Optional("proteus"): Any(), +}) + DiscretizationEntry = Map( { "name": Str(), @@ -133,7 +200,22 @@ Optional("from_weak_form"): Str(), Optional("finite_element"): FiniteElementDef, Optional("discretization"): Any(), + # with sub-problems, the default each one's own solver: refines "solver": SolverConfig, + # split into sub-problems, solved in this order each step + Optional("subproblems"): Seq(SubproblemDef), + # symbols the discretization supplies and the model may not use: + # the element size h, a smoothing width ε = 1.5 h + Optional("binds"): MapPattern(Str(), SymbolBinding), + # sharp model functions replaced by smooth ones: heaviside -> a + # smoothed Heaviside of width ε + Optional("regularize"): MapPattern(Str(), Map({ + "by": Str(), Optional("width"): Str()})), + Optional("step_control"): Map({ + Optional("type"): Str(), + Optional("cfl"): Float(), + Optional("step_exact"): Bool(), + }), Optional("flux_method"): Str(), Optional("numerical_flux"): Str(), Optional("discretization_details"): Str(), @@ -236,37 +318,41 @@ } ) -BoundaryConditionDef = Map( - { - "region": Str(), - "variable": Str(), - "type": BCTypeEnum, - Optional("value"): Float(), - # A value that varies along the boundary or in time: an expression - # in x, y, z, t and the coefficients; "(a, b)" for a vector variable. - Optional("formula"): Str(), - # The value as tabulated data rather than a formula -- a measured or - # design record, which is the one thing that cannot be written as an - # expression. Pairs of [t, value], t ascending. Linear between - # samples; outside the range the end value is held, because - # extrapolating a record past its end is not something to do - # quietly. - # - # Available to any type: with neumann this is a discharge - # hydrograph, with dirichlet a stage hydrograph or a measured inlet - # temperature. - Optional("series"): Seq(Seq(Float())), - # A robin condition, alpha*u + beta*(q.n) = g: alpha and beta here, - # g in value or formula like any other right-hand side. Expressions, - # so they may vary along the boundary. - Optional("alpha"): Str(), - Optional("beta"): Str(), - Optional("units"): Str(), - Optional("source"): Str(), - Optional("source_url"): Str(), - Optional("temporal"): Str(), - } -) +_CONDITION_FIELDS = { + "variable": Str(), + "type": BCTypeEnum, + Optional("value"): Float(), + # A value that varies along the boundary or in time: an expression + # in x, y, z, t and the coefficients; "(a, b)" for a vector variable. + Optional("formula"): Str(), + # The value as tabulated data rather than a formula -- a measured or + # design record, which is the one thing that cannot be written as an + # expression. Pairs of [t, value], t ascending. Linear between + # samples; outside the range the end value is held, because + # extrapolating a record past its end is not something to do + # quietly. + # + # Available to any type: with neumann this is a discharge + # hydrograph, with dirichlet a stage hydrograph or a measured inlet + # temperature. + Optional("series"): Seq(Seq(Float())), + # A robin condition, alpha*u + beta*(q.n) = g: alpha and beta here, + # g in value or formula like any other right-hand side. Expressions, + # so they may vary along the boundary. + Optional("alpha"): Str(), + Optional("beta"): Str(), + Optional("units"): Str(), + Optional("source"): Str(), + Optional("source_url"): Str(), + Optional("temporal"): Str(), +} + +BoundaryConditionDef = Map({"region": Str(), **_CONDITION_FIELDS}) + +# A condition set: conditions without a region, defined once in a model +# (free_slip, atmosphere) and applied to regions by a problem's +# boundary_sets. ymf.normalize.boundary_conditions expands them. +ConditionDef = Map(dict(_CONDITION_FIELDS)) # A domain is the notation short form -- "Ω = [0, 1] × [0, 1], t ∈ (0, 50]", # which means a box and is parsed by ymf.symbolic -- or a mapping tagged by @@ -308,6 +394,7 @@ # kind: file -- someone else's language Optional("format"): Str(), Optional("path"): Str(), + Optional("bounds"): Map({"lower": Seq(Float()), "upper": Seq(Float())}), } ) @@ -405,14 +492,31 @@ # "where" is a membership predicate, not geometry: "∂Ω" (all of # it), "y = 0 or y = H", "x = 0 and y = 0". adr_problem already # calls this thing "where" in the conditions it returns. - "boundary_regions": Seq(Map({"name": Str(), "where": Str()})), + # a region is named by a predicate (where) or by the mesh's boundary + # tags (tags: ints, or names the domain's boundary_tags resolve), + # exactly one of the two (ymf.references) + "boundary_regions": Seq(Map({"name": Str(), Optional("where"): Str(), + Optional("tags"): Seq(Int() | Str())})), Optional("initial_conditions"): Seq(InitialConditionDef), Optional("boundary_conditions"): Seq(BoundaryConditionDef), + # named sets of conditions, and the regions they apply to: a model + # says what free_slip and atmosphere mean, a problem where they hold + Optional("condition_sets"): MapPattern(Str(), Seq(ConditionDef)), + Optional("boundary_sets"): Seq(Map({"region": Str(), "set": Str()})), Optional("coefficients"): MapPattern(Str(), CoefficientDef), # The equations, one string per (vector or scalar) equation, in the # notation of ymf.symbolic.notation -- "∂T/∂t = ∇·(κ∇T) in Ω", # "ρ ∂v/∂t + ∇·(ρ v⊗v) − μΔv + ∇p = f" -- ordered like the unknowns. - Optional("equations"): Seq(Str()), + # Or named, saying which unknown each determines: a sub-problem of + # a discretization names the equations it solves, and an extending + # file merges onto an equation by its name. + Optional("equations"): Seq(Str() | Map({ + "name": Str(), + Optional("for"): Str(), + "equation": Str(), + Optional("provenance"): ProvenanceEnum, + Optional("notes"): Str(), + })), } ) @@ -420,6 +524,30 @@ # model (updated: + characteristic_scales, units, mesh_generation) # --------------------------------------------------------------------------- +# --- gauges: what a run measures as it goes (C22) ----------------------------- +# The shape planned in FORCE_HISTORY_GAUGES_PLAN step 4: proteus's point, +# line and line-integral gauges, and the boundary-integral gauge (forces and +# moments per boundary region) that replaces RANS2P's and RANS3PF's history +# files. What each kind needs is checked on load (ymf.references). +GaugeKindEnum = Enum(["point", "line", "line_integral", "boundary_integral"]) +GaugeQuantityEnum = Enum(["wetted_area", "force_pressure", "force_viscous", "force", "moment"]) + +GaugeDef = Map({ + Optional("name"): Str(), + "kind": GaugeKindEnum, + # point, line, line_integral: the fields sampled (unknowns or coefficients) + Optional("fields"): Seq(Str()), + Optional("points"): Seq(Seq(Float())), + # each line is two points, its ends + Optional("lines"): Seq(Seq(Seq(Float()))), + # boundary_integral: boundary regions, and what is integrated over each + Optional("tags"): Seq(Str()), + Optional("quantities"): Seq(GaugeQuantityEnum), + Optional("moment_point"): Seq(Float()), + Optional("sample_rate"): Float(), + Optional("file"): Str(), +}) + ModelDef = Map( { "name": Str(), @@ -430,6 +558,7 @@ Optional("dimensionless_numbers"): MapPattern(Str(), DimensionlessNumberDef), Optional("units"): UnitsDef, Optional("mesh_generation"): MeshGenerationDef, + Optional("gauges"): Seq(GaugeDef), } ) diff --git a/src/ymf/state.py b/src/ymf/state.py index 9f0639d..f8c0426 100644 --- a/src/ymf/state.py +++ b/src/ymf/state.py @@ -65,7 +65,7 @@ def branches(doc: Dict[str, Any]) -> List[Branch]: for entry in _discretizations(doc): name = entry.get("name") try: - lacking = lacks(entry) + lacking = lacks(entry, closure._domain_of(doc)) except ValueError as exc: # an unknown method out.append(Branch(name, "well-posed", (str(exc),))) continue @@ -74,7 +74,7 @@ def branches(doc: Dict[str, Any]) -> List[Branch]: "%s is missing %s: %s" % (name, m.what, m.hint) for m in lacking))) continue planned = tuple(closure.output_key(doc, name, r) - for r in closure.realizations(entry)) + for r in closure.realizations(entry, closure._domain_of(doc))) owed = tuple(k for k in planned if k not in present) if owed: out.append(Branch(name, "discretized", (), planned, owed)) diff --git a/src/ymf/symbolic/__init__.py b/src/ymf/symbolic/__init__.py index be1f5f0..8626472 100644 --- a/src/ymf/symbolic/__init__.py +++ b/src/ymf/symbolic/__init__.py @@ -8,8 +8,10 @@ from ymf.symbolic.language import (D, Dt, Space, SymbolicError, parse_equation, parse_expression) from ymf.symbolic.printer import equation_notation, to_notation -from ymf.symbolic.problem import adr_problem, strong_form_equations +from ymf.symbolic.problem import adr_problem, state_coefficients, strong_form_equations +from ymf.symbolic.split import split_problem __all__ = ["D", "Dt", "Space", "SymbolicError", "parse_equation", "parse_expression", "classify", "adr_form", "grad_name", "to_code", "adr_problem", - "to_notation", "equation_notation", "strong_form_equations"] + "to_notation", "equation_notation", "strong_form_equations", + "state_coefficients", "split_problem"] diff --git a/src/ymf/symbolic/adr.py b/src/ymf/symbolic/adr.py index 5662b6b..a34ba3a 100644 --- a/src/ymf/symbolic/adr.py +++ b/src/ymf/symbolic/adr.py @@ -39,6 +39,7 @@ from typing import Any, Dict, List, Optional, Sequence, Tuple import sympy +from sympy.core.function import AppliedUndef from sympy.printing.numpy import NumPyPrinter from ymf.symbolic.language import D, Dt, Space, SymbolicError @@ -64,9 +65,50 @@ def _print_Float(self, expr): # repr is the shortest text that reads back to the same double return repr(float(expr)) + # A known field (another sub-problem's unknown, lagged) is a function of + # the coordinates and time to sympy, and data to the solver: its name, + # and its gradient by grad_name, which the solver supplies the same way. + def _print_Function(self, expr): + if isinstance(expr, AppliedUndef): + return expr.func.__name__ + return super()._print_Function(expr) + + def _print_Derivative(self, expr): + field = expr.expr + if isinstance(field, AppliedUndef) and len(expr.variable_count) == 1: + variable, count = expr.variable_count[0] + if count == 1 and variable.name in ("x", "y", "z"): + return grad_name(field.func.__name__, "xyz".index(variable.name)) + if count == 1 and variable.name == "t": + raise SymbolicError( + "the time derivative of the known field %s: a sub-problem reads " + "another's field at one time, not its rate" % (field.func.__name__,)) + raise SymbolicError("%s: only first space derivatives of a known field have a " + "name (grad__)" % (expr,)) + + # numpy alone, nested, rather than sympy's functools.reduce: the code is + # evaluated with numpy and nothing else in scope + def _print_Min(self, expr): + return self._nest("numpy.minimum", expr.args) + + def _print_Max(self, expr): + return self._nest("numpy.maximum", expr.args) + + def _nest(self, function, args): + text = self._print(args[0]) + for arg in args[1:]: + text = "%s(%s, %s)" % (function, text, self._print(arg)) + return text + def to_code(expr) -> str: - return _Printer().doprint(sympy.sympify(expr)) + expr = sympy.sympify(expr) + if expr.has(sympy.DiracDelta): + raise SymbolicError( + "%s has no value at a point: it is the derivative of heaviside or sign. " + "A sharp function of an unknown needs a discretization that regularizes " + "it before a solver can use its derivative" % (expr,)) + return _Printer().doprint(expr) class Equation: @@ -126,7 +168,8 @@ def classify(residual, space: Space, where: str = "equation", forms. """ eq = Equation(space.dim) - allowed = set(space.component_symbols) | set(space.x) | {space.t} + allowed = (set(space.component_symbols) | set(space.x) | {space.t} + | set(space.parameters)) stray = sympy.sympify(residual).free_symbols - allowed if stray: raise SymbolicError("%s: undefined names %s (not a coordinate, an unknown " diff --git a/src/ymf/symbolic/language.py b/src/ymf/symbolic/language.py index bac6d4e..7790adc 100644 --- a/src/ymf/symbolic/language.py +++ b/src/ymf/symbolic/language.py @@ -63,10 +63,26 @@ class Space: expression in the coordinates and time); they are substituted as the equations are parsed, so a derivative of a variable coefficient comes out right. + + ``parameters`` are symbols the solver supplies, which an expression may + keep free (a discretization's element size ``h``). + + ``functions`` replace functions of the notation by name as expressions + are read (``heaviside`` -> a smoothed Heaviside of width ε): what is + written is regularized, and what differentiation produces is not. + + ``fields`` (name -> rank) are known fields: the unknowns of another + sub-problem, read here as data (a velocity lagged into a transport + solve). Each component is an undefined function of the coordinates and + time, so it is differentiated as a field, not as a constant; the code + printer writes it by name and its gradient as ``grad__``. """ def __init__(self, dim: int, unknowns: Dict[str, int], - coefficients: Dict[str, sympy.Expr] | None = None): + coefficients: Dict[str, sympy.Expr] | None = None, + fields: Dict[str, int] | None = None, + parameters: Iterable[sympy.Symbol] = (), + functions: Dict[str, object] | None = None): if dim not in (1, 2, 3): raise SymbolicError("space dimension must be 1, 2 or 3, not %r" % (dim,)) self.dim = dim @@ -89,6 +105,33 @@ def __init__(self, dim: int, unknowns: Dict[str, int], % (name, rank)) self.component_symbols = [sympy.Symbol(n, real=True) for n in self.components] self.coefficients = dict(coefficients or {}) + #: functions of the notation replaced for this space: a split + #: discretization's regularized heaviside and sign + self.functions = dict(functions or {}) + #: symbols the solver supplies and an expression may keep free: a + #: split discretization's bound symbols (h, the element size) + self.parameters = tuple(parameters) + #: known fields: name -> a function of (x.., t), or a column of them + self.known: Dict[str, object] = {} + #: and their scalar components by name ("v_0" -> v_0(x, y, t)) + self.known_components: Dict[str, object] = {} + args = tuple(self.x) + (self.t,) + for name, rank in (fields or {}).items(): + if name in self.unknowns: + raise SymbolicError("%r is both an unknown and a known field" % (name,)) + if rank == 0: + f = sympy.Function(name, real=True)(*args) + self.known[name] = f + self.known_components[name] = f + elif rank == 1: + parts = [sympy.Function("%s_%d" % (name, i), real=True)(*args) + for i in range(dim)] + self.known[name] = sympy.Matrix(parts) + self.known_components.update( + {"%s_%d" % (name, i): f for i, f in enumerate(parts)}) + else: + raise SymbolicError("field %r has rank %r; only 0 and 1 are supported" + % (name, rank)) # -- operators --------------------------------------------------------- @@ -177,6 +220,7 @@ def namespace(self) -> Dict[str, object]: ns[name] = self.x[axis] ns["d" + name] = (lambda e, a=axis: self.d(e, a)) ns.update(self.coefficients) + ns.update(self.known) ns.update(self._fields) return ns diff --git a/src/ymf/symbolic/notation.py b/src/ymf/symbolic/notation.py index da3cfc7..41e6483 100644 --- a/src/ymf/symbolic/notation.py +++ b/src/ymf/symbolic/notation.py @@ -24,7 +24,9 @@ ``·`` (a dot product of vectors, a product otherwise) and ``⊗`` (outer product, binding tighter than ``·`` and ``*``); tuples ``(a, b)`` are vectors; functions sin cos tan exp log sqrt sinh cosh - tanh abs. + tanh abs, the sharp ones heaviside sign min max (min and max of any + number of arguments), and on vectors and tensors norm(v) = sqrt(v·v) + and sym(A) = (A + Aᵀ)/2. Differential operators, each applied to the operand that follows it ``∇u`` gradient, ``∇·F`` divergence, ``Δu`` Laplacian, ``∂u/∂t`` and ``∂u/∂x`` partial derivatives; also the spelled-out @@ -131,9 +133,30 @@ def tokens(text: str) -> List[Tuple[str, str, int]]: "sin": sympy.sin, "cos": sympy.cos, "tan": sympy.tan, "exp": sympy.exp, "log": sympy.log, "sqrt": sympy.sqrt, "sinh": sympy.sinh, "cosh": sympy.cosh, "tanh": sympy.tanh, "abs": sympy.Abs, + # sharp functions: their derivatives are Dirac deltas, so a solver needs + # them regularized (a discretization's business); min and max take any + # number of arguments + "heaviside": sympy.Heaviside, "sign": sympy.sign, + "min": sympy.Min, "max": sympy.Max, } +def _norm(v): + """The Euclidean norm of a vector, sqrt(v·v); of a scalar, its absolute value.""" + if isinstance(v, sympy.MatrixBase): + if v.shape[1] != 1: + raise ValueError("norm needs a vector, got a %dx%d tensor" % v.shape) + return sympy.sqrt(sum(e ** 2 for e in v)) + return sympy.Abs(v) + + +def _sym(a): + """The symmetric part of a square tensor, (A + Aᵀ)/2: sym(∇v) is the strain rate.""" + if not isinstance(a, sympy.MatrixBase) or a.shape[0] != a.shape[1] or a.shape[1] == 1: + raise ValueError("sym needs a square tensor, such as ∇v") + return (a + a.T) / 2 + + class Parser: """Parse notation against a :class:`~ymf.symbolic.language.Space`.""" @@ -147,6 +170,8 @@ def __init__(self, space, locals_: Optional[Dict[str, Any]] = None): # works on them -- the scalarized ADR form -- can be written back # (ymf.symbolic.printer); the vector itself is still just "v" self.names.update(zip(space.components, space.component_symbols)) + self.names.update(space.known_components) + self.names.update(space.known) self.names.update(space._fields) self.names.update(locals_ or {}) self.functions = dict(_FUNCTIONS) @@ -154,9 +179,11 @@ def __init__(self, space, locals_: Optional[Dict[str, Any]] = None): "div": space.div, "grad": space.grad, "lap": lambda e: space.div(space.grad(e)), "dt": space.dt, "dot": self._dot, "outer": space.outer, "vec": lambda *xs: sympy.Matrix([sympy.sympify(v) for v in xs]), + "norm": _norm, "sym": _sym, }) for axis, n in enumerate("xyz"[:space.dim]): self.functions["d" + n] = (lambda e, a=axis: space.d(e, a)) + self.functions.update(getattr(space, "functions", {})) # -- driving ------------------------------------------------------------- diff --git a/src/ymf/symbolic/printer.py b/src/ymf/symbolic/printer.py index 7da1b67..f2844ca 100644 --- a/src/ymf/symbolic/printer.py +++ b/src/ymf/symbolic/printer.py @@ -17,12 +17,13 @@ (``u``, ``v_0``), and any other symbol by its name; * ``+ - * /`` and powers, with ``² ³ ⁴`` for small integer exponents and ``^`` otherwise; ``sqrt`` for a square root; -* the functions sin cos tan exp log sqrt sinh cosh tanh abs; +* the functions sin cos tan exp log sqrt sinh cosh tanh abs heaviside + sign, and min and max of any number of arguments; * derivatives: ``D(e, axis)`` as ``∂e/∂x``, ``Dt(e)`` as ``∂e/∂t``; * a column vector as a tuple ``(a, b)``. Anything else (a tensor, an integral, a function the notation lacks, or -what sympy forms only from complex values -- ``I``, ``sign`` -- out of a +what sympy forms only from complex values -- ``I`` -- out of a root of a negative quantity) is refused with :class:`NotationError` rather than written in a form the parser would read differently. @@ -49,8 +50,9 @@ _FUNCTIONS = { sympy.sin: "sin", sympy.cos: "cos", sympy.tan: "tan", sympy.exp: "exp", sympy.log: "log", sympy.sinh: "sinh", sympy.cosh: "cosh", sympy.tanh: "tanh", - sympy.Abs: "abs", + sympy.Abs: "abs", sympy.Heaviside: "heaviside", sympy.sign: "sign", } +_VARIADIC = {sympy.Min: "min", sympy.Max: "max"} # precedence of what a printed form is, for deciding on parentheses _ATOM, _POWER, _PRODUCT, _SUM = 4, 3, 2, 1 @@ -115,6 +117,9 @@ def _text(e) -> tuple: for function, name in _FUNCTIONS.items(): if isinstance(e, function): return "%s(%s)" % (name, _text(e.args[0])[0]), _ATOM + for function, name in _VARIADIC.items(): + if isinstance(e, function): + return "%s(%s)" % (name, ", ".join(_text(a)[0] for a in e.args)), _ATOM raise NotationError("%s (%s) has no notation" % (e, type(e).__name__)) diff --git a/src/ymf/symbolic/problem.py b/src/ymf/symbolic/problem.py index 938e7ea..64675dd 100644 --- a/src/ymf/symbolic/problem.py +++ b/src/ymf/symbolic/problem.py @@ -24,15 +24,16 @@ from __future__ import annotations import math -from typing import Any, Dict, List, Optional, Tuple +from typing import Any, Dict, List, Optional, Sequence, Tuple import sympy from ymf.domain import box_bounds from ymf.symbolic.adr import adr_form, to_code -from ymf.symbolic.language import Space, SymbolicError, parse_equation, parse_expression +from ymf.symbolic.language import (D, Dt, Space, SymbolicError, parse_equation, + parse_expression) -__all__ = ["adr_problem", "strong_form_equations"] +__all__ = ["adr_problem", "strong_form_equations", "state_coefficients"] def _number(value): @@ -48,8 +49,12 @@ def _number(value): return None -def _coefficients(spec: Dict[str, Any], dim: int, - strict: bool = True) -> Dict[str, sympy.Expr]: +def _coefficients(spec: Dict[str, Any], dim: int, strict: bool = True, + unknowns: Optional[Dict[str, int]] = None, + fields: Optional[Dict[str, int]] = None, + preset: Optional[Dict[str, Any]] = None, + functions: Optional[Dict[str, Any]] = None, + parameters: Sequence[sympy.Symbol] = ()) -> Dict[str, sympy.Expr]: """Resolve coefficient values and formulas; formulas may use each other. ``strict`` decides what an unresolvable formula means. Building a problem @@ -57,10 +62,33 @@ def _coefficients(spec: Dict[str, Any], dim: int, defines cannot be run. Reading the same spec to report on it -- which is what :func:`coefficient_numbers` does for ``ymf.units`` -- it is not, and the entry is simply left out. + + ``unknowns`` (name -> rank) lets a formula name the unknowns too: a + density ``ρw (1 - χ) + ρa χ`` with ``χ`` an unknown, a conductivity + ``1 + u²``. Such a coefficient is a function of the state, substituted + into the equations like any other, so its dependence on the unknowns + shows up in the ADR form's ``depends_on`` and derivatives. Without + ``unknowns`` it is left unresolved, which is what a reader that wants + numbers needs. A formula may name the unknowns, not their derivatives. + + ``fields`` (name -> rank) are known fields a formula may name too (a + sub-problem's lagged reads, see ymf.symbolic.split), and ``preset`` gives + coefficients their value outright, in place of their formula (a + coefficient read through the sub-problem that carries it: χ -> χh). """ - values: Dict[str, sympy.Expr] = {} + unknowns = dict(unknowns or {}) + clash = sorted(set(unknowns) & set(spec)) + if clash: + raise SymbolicError("coefficient %s has the name of an unknown; a name is one " + "or the other" % ", ".join(clash)) + state = Space(dim, unknowns) + allowed = ({"x", "y", "z", "t"} | set(state.components) + | {p.name for p in parameters}) + values: Dict[str, sympy.Expr] = dict(preset or {}) formulas: Dict[str, str] = {} for name, entry in spec.items(): + if name in values: + continue # preset: read from elsewhere, not computed if isinstance(entry, dict): if entry.get("value") is not None: v = _number(entry["value"]) @@ -81,13 +109,19 @@ def _coefficients(spec: Dict[str, Any], dim: int, progressed = False for name, text in list(formulas.items()): try: - value = parse_expression(text, Space(dim, {}, values)) + value = parse_expression(text, Space(dim, unknowns, values, fields, + functions=functions)) except SymbolicError: continue - stray = (value.free_symbols if not isinstance(value, sympy.MatrixBase) - else set().union(*[e.free_symbols for e in value])) - if {s.name for s in stray} - {"x", "y", "z", "t"}: + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + stray = set().union(*[sympy.sympify(e).free_symbols for e in entries]) + if {s.name for s in stray} - allowed: continue # depends on a coefficient not resolved yet + if any(sympy.sympify(e).has(D, Dt) for e in entries): + raise SymbolicError( + "coefficient %s: %r takes a derivative of an unknown; a coefficient " + "may depend on the unknowns, not (yet) on their derivatives" + % (name, text)) values[name] = value del formulas[name] progressed = True @@ -127,8 +161,55 @@ def coefficient_numbers(spec: Dict[str, Any], dim: int = 3) -> Dict[str, float]: return numbers +def _given(value, space: Space, where: str): + """``value``, refused if it depends on the unknowns. + + Boundary data, initial conditions and exact solutions are evaluated by + the solver at points, in x, y, z and t alone. Since a coefficient may now + be a function of the state, naming one there could slip an unknown in. + """ + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + named = set().union(*[sympy.sympify(e).free_symbols for e in entries]) + state = sorted(s.name for s in named & set(space.component_symbols)) + if state: + raise SymbolicError("%s depends on the unknowns (%s), directly or through a " + "coefficient; give it in x, y, z and t" + % (where, ", ".join(state))) + return value + + +def state_coefficients(strong_form: Dict[str, Any], dim: int = 3) -> Dict[str, List[str]]: + """The coefficients that are functions of the state, and of what. + + ``{coefficient: [component, ...]}`` for every coefficient whose formula + names an unknown, directly or through another coefficient. Such a + coefficient is not a number, so a scale or a dimensionless number built + on it has no value; ``ymf.units`` says so rather than leave it out + silently. Coefficients that do not resolve are left out, as in + :func:`coefficient_numbers`. + """ + ranks = {} + for u in strong_form.get("unknowns") or []: + if isinstance(u, str): + ranks[u] = 0 + else: + ranks[u["name"]] = int(u.get("rank", 0)) + values = _coefficients(strong_form.get("coefficients") or {}, dim, strict=False, + unknowns=ranks) + components = {s.name: s for s in Space(dim, ranks).component_symbols} + out = {} + for name, value in values.items(): + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + named = set().union(*[sympy.sympify(e).free_symbols for e in entries]) + state = sorted(s.name for s in named if s.name in components) + if state: + out[name] = state + return out + + def _components(value, name: str, rank: int, space: Space, where: str) -> Dict[str, str]: """Per-component code strings for a scalar or vector value.""" + _given(value, space, where) if rank == 0: if isinstance(value, sympy.MatrixBase): raise SymbolicError("%s: %s is a scalar but got a vector" % (where, name)) @@ -147,13 +228,15 @@ def _robin_part(bc, key, space): "robin condition on %s: give alpha and beta, the coefficients of " "u and of the outward flux in alpha u + beta (q.n) = g" % (bc["region"],)) - return to_code(parse_expression(str(text), space)) + return to_code(_given(parse_expression(str(text), space), space, + "robin condition on %s: %s" % (bc["region"], key))) def _robin_rhs(bc, space): """g, from formula or value like any other right-hand side.""" if bc.get("formula") is not None: - return to_code(parse_expression(bc["formula"], space)) + return to_code(_given(parse_expression(bc["formula"], space), space, + "robin condition on %s: g" % (bc["region"],))) if bc.get("value") is not None: return to_code(_number(bc["value"])) raise SymbolicError( @@ -189,6 +272,37 @@ def _series(bc, name, rank): return [[float(t), float(v)] for t, v in series] +def region_tags(region: Dict[str, Any], domain: Any) -> List[int]: + """A tag region's boundary flags: its ints, and its names resolved + through the domain's ``boundary_tags``.""" + names = (domain.get("boundary_tags") if isinstance(domain, dict) else None) or {} + flags = [] + for tag in region["tags"]: + if isinstance(tag, int) or (isinstance(tag, str) and tag.lstrip("-").isdigit()): + flags.append(int(tag)) + elif tag in names: + flags.append(int(names[tag])) + else: + raise SymbolicError("boundary region %s: tag %r is not one of the domain's " + "boundary_tags (have: %s)" + % (region["name"], tag, ", ".join(sorted(names)) or "none")) + return flags + + +def _selector(region: Dict[str, Any], domain: Any, space: Space) -> Dict[str, Any]: + """How a solver finds a region's boundary points. + + ``{"where": code}`` for a predicate (None: the whole boundary), or + ``{"flags": [ints]}`` for a tag region, which the solver tests against the + boundary flag of the point. A tag region has no ``where`` key at all, so a + solver that does not know ``flags`` fails on it rather than reading a + missing predicate as "the whole boundary". + """ + if region.get("tags") is not None: + return {"flags": region_tags(region, domain)} + return {"where": _where(region["where"], space)} + + def _where(text: str, space: Space) -> Optional[str]: """A region's ``where`` predicate -> boolean code using ``tol`` (None: the whole boundary).""" from ymf.symbolic.notation import NotationError, Parser @@ -219,6 +333,68 @@ def _domain_dimension(text: str) -> int: return count +def periodic_faces_match(region: Dict[str, Any], axes: List[int], lower, upper, + name: str) -> None: + """Check that a region's periodic faces are translations of each other. + + A periodic axis pairs each boundary node on the lower face with the node + the translation along that axis puts on the upper face; a mesh has such + partners only if the boundary pieces on the two faces match. So the + segments (2-D) or facets (3-D) lying in the lower face, moved by the + box's length along the axis, must be exactly those lying in the upper + face, to a tolerance of 1e-8 of the box. (The mesher must then add no + points to them; that, and the partners themselves, are the solver's to + check after meshing.) + """ + vertices = [[float(c) for c in v] for v in region["vertices"]] + dim = len(vertices[0]) + scale = max(u - l for u, l in zip(upper, lower)) + tol = 1e-8 * scale + if dim == 2: + pieces = [list(s) for s in region.get("segments") or []] + if not pieces: + raise SymbolicError("periodic condition on %s: the region gives no segments to " + "match" % (name,)) + else: + pieces = [] + for facet in region.get("facets") or []: + # a PLC facet is a list of polygons; a bare polygon is allowed too + polygons = facet if facet and isinstance(facet[0], (list, tuple)) else [facet] + for polygon in polygons: + pieces.append(list(polygon)) + if not pieces: + raise SymbolicError("periodic condition on %s: the region gives no facets to " + "match" % (name,)) + + def key(points): + return frozenset(tuple(round(c / tol) for c in p) for p in points) + + for a in axes: + shift = upper[a] - lower[a] + on = {side: [] for side in ("lower", "upper")} + for piece in pieces: + points = [vertices[i] for i in piece] + for side, face in (("lower", lower[a]), ("upper", upper[a])): + if all(abs(p[a] - face) <= tol for p in points): + on[side].append(points) + moved = {key([[c + (shift if i == a else 0.0) for i, c in enumerate(p)] + for p in points]) for points in on["lower"]} + there = {key(points) for points in on["upper"]} + unmatched = len(moved ^ there) + if not on["lower"] or not on["upper"]: + raise SymbolicError( + "periodic condition on %s: no boundary %s lies in the %s face of axis %s" + % (name, "segment" if dim == 2 else "facet", + "lower" if not on["lower"] else "upper", "xyz"[a])) + if unmatched: + raise SymbolicError( + "periodic condition on %s: the faces of axis %s do not match: %d of %d " + "boundary %ss on them have no translated partner (the region's lower and " + "upper faces must be translations of each other)" + % (name, "xyz"[a], unmatched, len(on["lower"]) + len(on["upper"]), + "segment" if dim == 2 else "facet")) + + def _periodic_axes(where, lower, upper, space: Space, region: str) -> List[int]: """The axes whose two opposite faces both lie in a periodic region.""" code = _where(where, space) @@ -247,23 +423,77 @@ def _periodic_axes(where, lower, upper, space: Space, region: str) -> List[int]: return axes -def _posed(doc: Dict[str, Any]): - """What adr_problem and strong_form_equations share. +def _geometry_dimension(region: Dict[str, Any]) -> Optional[int]: + """A discrete domain's space dimension, from its geometry's type.""" + kind = ((region.get("geometry") or {}).get("type") or "").upper() + return {"XY": 2, "XYZ": 3}.get(kind) - Returns ``(strong, box, region, ranks, coefficients, space, residuals)``. + +def _inline_nodes(region: Dict[str, Any]) -> Optional[List[List[float]]]: + """A discrete domain's node coordinates, when its geometry holds them + inline (``values``); None when they are in a file (an HDF5 reference), + which the solver reads.""" + item = (region.get("geometry") or {}).get("data_item") or {} + values = item.get("values") + if not values: + return None + flat = [float(v) for v in (values if not isinstance(values[0], (list, tuple)) + else [c for row in values for c in row])] + dim = _geometry_dimension(region) or int((item.get("dimensions") or [0, 3])[-1]) + return [flat[i:i + dim] for i in range(0, len(flat), dim)] + + +def _frame(doc: Dict[str, Any]): + """The strong form's setting, before any equation is read. + + Returns ``(strong, box, region, ranks, dim)``: the strong form, the + domain's box (the bounding box of a region), the region itself (None + for a box), each unknown's rank, and the space dimension. Shared by + :func:`adr_problem` and the split reader, ymf.symbolic.split. """ problem = doc["model"] strong = problem["strong_form"] domain = strong["domain"] box = box_bounds(domain) region = None + if box is None and not isinstance(domain, str) and domain.get("kind") in ("file", + "discrete"): + # someone else's language, which the solver reads (or meshes, a + # geometry): carried whole. Its bounds set the boundary tolerance's + # scale; without them the solver takes them from the mesh, so that + # posing follows meshing for this kind alone (cekees, D2). + # A discrete domain is a mesh already (C7.1: a mesh is a type of + # domain), solved as it stands; its bounds are its nodes' (D4), when + # they are in the document, and the solver's to take otherwise. + from ymf.closure import plain + region = plain(domain) + bounds = region.get("bounds") + nodes = _inline_nodes(region) if region["kind"] == "discrete" else None + if bounds is None and nodes: + bounds = {"lower": [min(c) for c in zip(*nodes)], + "upper": [max(c) for c in zip(*nodes)]} + if bounds is not None: + box = [[float(c) for c in bounds["lower"]], [float(c) for c in bounds["upper"]]] + dim = len(box[0]) + elif region.get("dimension") or _geometry_dimension(region): + dim = int(region.get("dimension") or _geometry_dimension(region)) + else: + raise SymbolicError("strong_form.domain: a %s domain gives its bounds, or at " + "least its dimension, so the problem can be posed before " + "the mesh is read" % (region["kind"],)) + if not strong.get("equations"): + raise SymbolicError("strong_form.equations is needed to run a spec") + ranks = {} + for u in strong["unknowns"]: + ranks[u if isinstance(u, str) else u["name"]] = ( + 0 if isinstance(u, str) else int(u.get("rank", 0))) + return strong, box, region, ranks, dim if box is None and not isinstance(domain, str): if domain.get("kind") != "piecewise_linear": raise SymbolicError( - "strong_form.domain: adr_problem can run a box or a " - "piecewise_linear region; this one is kind %r. Other kinds " - "validate and are carried, but nothing solves them yet." - % (domain.get("kind"),)) + "strong_form.domain: adr_problem can run a box, a piecewise_linear " + "region or a file; this one is kind %r. Other kinds validate and are " + "carried, but nothing solves them yet." % (domain.get("kind"),)) # A boundary representation: the solver meshes it (its # discretization gives mesh.h). The bounding box stands in for the # box where one is needed: the boundary tolerance's scale. @@ -280,10 +510,24 @@ def _posed(doc: Dict[str, Any]): ranks[u] = 0 else: ranks[u["name"]] = int(u.get("rank", 0)) - coefficients = _coefficients(strong.get("coefficients") or {}, dim) + return strong, box, region, ranks, dim + + +def _posed(doc: Dict[str, Any]): + """What adr_problem and strong_form_equations share. + + Returns ``(strong, box, region, ranks, coefficients, space, residuals)``. + """ + strong, box, region, ranks, dim = _frame(doc) + coefficients = _coefficients(strong.get("coefficients") or {}, dim, unknowns=ranks) space = Space(dim, ranks, coefficients) residuals: List[sympy.Expr] = [] - for text in strong["equations"]: + from ymf.normalize import equations_in_order + try: + texts = equations_in_order(strong) + except ValueError as exc: + raise SymbolicError(str(exc)) from exc + for text in texts: residuals.extend(parse_equation(text, space)) return strong, box, region, ranks, coefficients, space, residuals @@ -314,6 +558,11 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] lower, upper = box interval = domain.get("time") time = tuple(interval) if interval else None + elif isinstance(domain, dict): + # a file domain without bounds: the solver takes them from the mesh + lower = upper = None + interval = domain.get("time") + time = tuple(interval) if interval else None else: try: lower, upper, time = Parser(Space(dim, {}, coefficients)).domain(domain) @@ -322,29 +571,47 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] adr = adr_form(residuals, space, hamiltonian_gradients=tuple(hamiltonian_gradients)) - regions = {r["name"]: r["where"] for r in strong.get("boundary_regions", [])} + regions = {r["name"]: r for r in strong.get("boundary_regions", [])} dirichlet = [] neumann = [] robin = [] periodic = [] - for bc in strong.get("boundary_conditions") or []: + from ymf.normalize import boundary_conditions + try: + conditions = boundary_conditions(strong) + except ValueError as exc: + raise SymbolicError(str(exc)) from exc + for bc in conditions: if bc["type"] not in ("dirichlet", "neumann", "periodic", "robin"): raise SymbolicError("boundary condition on %s: %s is not a kind of " - "boundary condition that runs" - % (bc["region"], bc["type"])) + "boundary condition that runs in one ADR system (%s " + "conditions are for a split discretization's solver)" + % (bc["region"], bc["type"], bc["type"])) name = bc["variable"] if name not in ranks: - raise SymbolicError("boundary condition on %s: %s is not an unknown" - % (bc["region"], name)) + raise SymbolicError("boundary condition on %s: %s is not an unknown%s" + % (bc["region"], name, + " (a condition on a coefficient is for a split " + "discretization's solver)" + if name in (strong.get("coefficients") or {}) else "")) if bc["region"] not in regions: raise SymbolicError("boundary condition: no boundary region named %r" % bc["region"]) - where = regions[bc["region"]] + selector = _selector(regions[bc["region"]], domain, space) + where = regions[bc["region"]].get("where") if bc["type"] == "periodic": - if region is not None: + if where is None: raise SymbolicError( - "periodic condition on %s: periodic faces are a box's; " - "this domain is piecewise_linear" % (bc["region"],)) + "periodic condition on %s: periodic faces are found by position, " + "so the region needs where:, not tags:" % (bc["region"],)) + if region is not None and region.get("kind") != "piecewise_linear": + raise SymbolicError( + "periodic condition on %s: a %s domain cannot be periodic (a mesh's " + "periodic pairing, or Gmsh's Periodic constraint, is not carried); a " + "box or a piecewise_linear region whose faces match can" + % (bc["region"], region.get("kind"))) axes = _periodic_axes(where, lower, upper, space, bc["region"]) + if region is not None: + periodic_faces_match(region, axes, lower, upper, bc["region"]) comps = [name] if ranks[name] == 0 else ["%s_%d" % (name, i) for i in range(dim)] periodic.extend({"component": c, "region": bc["region"], "axes": axes} for c in comps) continue @@ -354,8 +621,7 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] "robin condition on %s: %s is a vector. alpha u + beta q.n = g " "is one relation, so give one per component." % (bc["region"], name)) - robin.append({"component": name, "region": bc["region"], - "where": _where(where, space), + robin.append({"component": name, "region": bc["region"], **selector, "alpha": _robin_part(bc, "alpha", space), "beta": _robin_part(bc, "beta", space), "g": _robin_rhs(bc, space)}) @@ -366,8 +632,7 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] # neumann, a stage hydrograph or a measured inlet value under # dirichlet. One number per time, so a vector wants one # condition per component. - entries.append({"component": name, "region": bc["region"], - "where": _where(where, space), + entries.append({"component": name, "region": bc["region"], **selector, "series": _series(bc, name, ranks[name])}) continue if bc.get("formula") is not None: @@ -380,8 +645,7 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] raise SymbolicError("boundary condition on %s: give value or formula" % bc["region"]) for component, code in _components(value, name, ranks[name], space, "boundary condition on %s" % bc["region"]).items(): - entries.append({"component": component, "region": bc["region"], - "where": _where(where, space), + entries.append({"component": component, "region": bc["region"], **selector, "value": code}) initial = {} @@ -432,6 +696,58 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] # operator; errors in them are measured modulo constants. Absent # when there are none. out["up_to_constant"] = up_to_constant + named = _state_coefficient_data(coefficients, space) + if named: + # The coefficients that follow the state, by name. The ADR form above + # has them substituted, which is all a solver of one system needs; a + # code generator or a dispatcher wants the name too, to emit rho as + # one inline function of the state rather than its expansion in every + # slot. Absent when there are none, so a problem without them reads + # exactly as before. + out["coefficients"] = named + return out + + +def _state_coefficient_data(coefficients, space: Space) -> Dict[str, Any]: + """``{name: {code, depends_on, d}}`` for each coefficient of the state. + + ``code`` is the numpy code of the coefficient (a list for a vector), + ``depends_on`` maps each component it depends on to ``linear`` or + ``nonlinear``, and ``d`` holds its derivative in each. A sharp function + of an unknown (``heaviside(φ)``) has no derivative with a value at a + point; it is marked ``sharp: true`` and its ``d`` is left out, for the + discretization to regularize. + """ + from sympy.core.function import AppliedUndef + from ymf.symbolic.adr import _flag + unknowns = set(space.component_symbols) + out: Dict[str, Any] = {} + for name, value in coefficients.items(): + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + entries = [sympy.sympify(e) for e in entries] + fields = sorted({f.func.__name__ for e in entries for f in e.atoms(AppliedUndef)}) + if not (set().union(*[e.free_symbols for e in entries]) & unknowns or fields): + continue + flags: Dict[str, str] = {} + for e in entries: + for component, flag in _flag(e, space).items(): + if flags.get(component) != "nonlinear": + flags[component] = flag + vector = isinstance(value, sympy.MatrixBase) + record: Dict[str, Any] = { + "code": [to_code(e) for e in entries] if vector else to_code(entries[0]), + "depends_on": flags} + if fields: + # known fields it reads (a sub-problem's lagged data, ymf.symbolic.split) + record["reads"] = fields + try: + record["d"] = { + c: ([to_code(sympy.diff(e, s)) for e in entries] if vector + else to_code(sympy.diff(entries[0], s))) + for c, s in zip(space.components, space.component_symbols) if c in flags} + except SymbolicError: + record["sharp"] = True + out[name] = record return out diff --git a/src/ymf/symbolic/split.py b/src/ymf/symbolic/split.py new file mode 100644 index 0000000..1a04a09 --- /dev/null +++ b/src/ymf/symbolic/split.py @@ -0,0 +1,423 @@ +"""A split discretization -> one plain-data problem per sub-problem. + +A discretization with ``subproblems`` solves the model as a sequence of +smaller problems in each step (see the guide, "Splitting a discretization +into sub-problems"). The whole model need not be one ADR system -- two-phase +level-set flow is not: ``v·∇φ`` is a Hamiltonian in another unknown -- but +each sub-problem is, once what it reads from the others is data. So each one +is classified on its own: + +* its **unknowns** are the unknowns of the model equations it ``solves`` + (their ``for``), and the unknown of its own ``equation`` (its ``for``, by + default the unknown it ``introduces``); +* every other field it reads must be declared in ``lags``; it is a **known + field**, a function of the coordinates and time, which the plain data + names (``v_0``) and whose gradient it names ``grad__``; a + coefficient lagged ``via`` the sub-problem that carries it (``χ`` via the + one that introduces ``χh``, approximating ``χ``) is that carrier; +* the symbols the discretization ``binds`` (``h``, ``ε = 1.5 h``) are + substituted, and the sharp functions it ``regularize``\\ s become smooth + ones: Proteus's smoothed Heaviside of width ε, + ``H_ε(φ) = (1 + φ/ε + sin(πφ/ε)/π)/2`` on ``|φ| ≤ ε`` (``gf.H`` in + ``equivalent_polynomials.h``), whose derivative is its ``gf.D``, and the + smoothed sign ``2 H_ε − 1`` of its redistancing. + +The result is plain data, like :func:`ymf.symbolic.adr_problem`'s: what a +dispatcher maps onto a solver's modules, and what a kernel generator reads. +""" + +from __future__ import annotations + +from typing import Any, Dict, List, Optional, Tuple + +import sympy + +from ymf.symbolic.adr import adr_form, to_code +from ymf.symbolic.language import Space, SymbolicError, parse_equation, parse_expression + +__all__ = ["split_problem", "smoothed_heaviside", "smoothed_sign", "REGULARIZATIONS"] + + +def smoothed_heaviside(phi, eps): + """Proteus's ``gf.H(eps, phi)``: 0 below -ε, 1 above ε, smooth between.""" + return sympy.Piecewise( + (sympy.Integer(1), phi > eps), + (sympy.Integer(0), phi < -eps), + ((1 + phi / eps + sympy.sin(sympy.pi * phi / eps) / sympy.pi) / 2, True)) + + +def smoothed_sign(phi, eps): + """``2 H_ε(φ) − 1``, the sign Proteus's redistancing uses.""" + return 2 * smoothed_heaviside(phi, eps) - 1 + + +#: sharp function -> {regularization name: f(argument, width)} +REGULARIZATIONS = { + "heaviside": {"smoothed_heaviside": smoothed_heaviside}, + "sign": {"smoothed_sign": smoothed_sign}, +} + + +def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: + """The plain-data problem of each sub-problem of ``discretization``. + + ``doc`` is a validated, composed spec (its names were checked on load, + ymf.references). Returns:: + + {"discretization": name, "dim": d, + "domain": {"lower": [...], "upper": [...]}, "time": [t0, t1] | None, + "bound": {symbol: meaning or code}, + "initial": {component: code}, # the model unknowns' + "subproblems": [{"name", "unknowns", "adr", "reads", "coefficients", + "regularized", "initial", "boundary_conditions", + ...}, ...]} + + A sub-problem that solves nothing it can classify (one that only spins + the run up, or enforces a constraint it gives no equation for) has + ``"classified": False`` and the reason instead of ``adr``. + """ + from ymf.discretizations import subproblem_solver + from ymf.normalize import normalize_equations + from ymf.symbolic.problem import _coefficients, _frame, _state_coefficient_data + entries = (doc.get("solution_paths") or {}).get("discretizations") or [] + entry = next((d for d in entries if d.get("name") == discretization), None) + if entry is None: + raise SymbolicError("no discretization named %r (have: %s)" + % (discretization, ", ".join(d.get("name", "?") for d in entries))) + subproblems = entry.get("subproblems") or [] + if not subproblems: + raise SymbolicError("discretization %s is not split into sub-problems; " + "adr_problem reads it as one system" % (discretization,)) + strong, box, region, ranks, dim = _frame(doc) + equations = {e["name"]: e for e in normalize_equations(strong.get("equations") or [])} + spec = strong.get("coefficients") or {} + introduced = {s["introduces"]["name"]: s for s in subproblems if s.get("introduces")} + approximates = {s["introduces"]["approximates"]: s["introduces"]["name"] + for s in subproblems if (s.get("introduces") or {}).get("approximates")} + entry_bound = _bind(entry.get("binds") or {}, {}, dim, "the discretization's binds") + + numbers = _coefficients(spec, dim, strict=False) + domain = strong["domain"] + time = None + if box is not None: + lower, upper = box + if isinstance(domain, dict) and domain.get("time"): + time = list(domain["time"]) + else: + from ymf.symbolic.notation import Parser + lower, upper, interval = Parser(Space(dim, {}, numbers)).domain(domain) + time = list(interval) if interval else None + out_subs = [] + for sub in subproblems: + where = "sub-problem %s" % (sub["name"],) + bound = _bind(sub.get("binds") or {}, entry_bound, dim, where + " binds") + record: Dict[str, Any] = {"name": sub["name"]} + for key in ("solves", "enforces", "writes_back", "pseudo_time", "spin_up", + "finite_element", "time_integration", "stabilization", "proteus"): + if key in sub: + record[key] = sub[key] + record["solver"] = subproblem_solver(entry, sub) + + # its unknowns, and the text of every equation it classifies + own: Dict[str, int] = {} + texts: List[Tuple[str, str]] = [] + for name in sub.get("solves") or []: + equation = equations[name] + if equation["for"] is None: + raise SymbolicError("%s solves %s, which does not say which unknown it is " + "for; give the model's equations 'for'" % (where, name)) + own[equation["for"]] = ranks[equation["for"]] + texts.append((name, equation["equation"])) + if sub.get("equation"): + target = sub.get("for") or (sub.get("introduces") or {}).get("name") + if target is None: + raise SymbolicError("%s gives an equation but not the unknown it is for: " + "give for:" % (where,)) + own[target] = ranks.get(target, 0) + texts.append(("equation", sub["equation"])) + + # what it reads: known fields, and coefficients read through a carrier + fields: Dict[str, int] = {} + preset: Dict[str, Any] = {} + reads: Dict[str, Any] = {} + for name, lag in (sub.get("lags") or {}).items(): + lag = lag if isinstance(lag, dict) else {"at": lag} + target = lag.get("of", name) + read = {"at": lag["at"]} + if "via" in lag: + read["via"] = lag["via"] + if "of" in lag: + read["of"] = lag["of"] + if target in ranks or target in introduced: + fields[name] = ranks.get(target, 0) + elif target in spec and target in approximates: + carrier = approximates[target] + fields[carrier] = 0 + read["carried_by"] = carrier + reads[name] = read + known = Space(dim, own, {}, fields).known + for name, read in reads.items(): + field = read.get("carried_by", name) + if field in known: + if "carried_by" in read: + preset[name] = known[field] + read["components"] = ([field] if fields[field] == 0 else + ["%s_%d" % (field, i) for i in range(dim)]) + if not own: + record["classified"] = False + record["why"] = ("it spins the run up" if sub.get("spin_up") else + "it enforces %r and gives no equation for it" % sub.get("enforces")) + record["reads"] = reads + out_subs.append(record) + continue + + smooth = _regularizations(entry, sub, bound, dim, where) + functions = {fn: (lambda a, *rest, s=smoother, w=width: s(a, w)) + for fn, (_, width, smoother) in smooth.items()} + supplied = [s for s in bound.values() if isinstance(s, sympy.Symbol)] + coefficients = _coefficients(spec, dim, strict=False, unknowns=own, + fields=fields, preset=preset, functions=functions, + parameters=supplied) + declared = (set(spec) | set(ranks) | set(introduced) | set(bound) + | set(sub.get("lags") or {})) + used = _used([t for _, t in texts], spec, declared, carried=set(preset)) + values = dict(coefficients) + values.update(bound) + space = Space(dim, own, values, fields, supplied, functions) + residuals: List[sympy.Expr] = [] + for label, text in texts: + try: + residuals.extend(parse_equation(text, space)) + except SymbolicError as exc: + raise SymbolicError( + "%s, %s: %s (a field a sub-problem reads must be one of its unknowns " + "or declared in its lags)" % (where, label, exc)) from exc + written = [t for _, t in texts] + [ + (spec[c].get("formula") if isinstance(spec[c], dict) else spec[c]) or "" + for c in used if c not in preset] + [sub.get("initial_condition") or ""] + sharp = _written(written, smooth) + try: + record["adr"] = adr_form(residuals, space) + except SymbolicError as exc: + raise SymbolicError("%s: %s" % (where, exc)) from exc + record["unknowns"] = {name: ([name] if rank == 0 else + ["%s_%d" % (name, i) for i in range(dim)]) + for name, rank in own.items()} + record["reads"] = reads + named = _state_coefficient_data( + {k: v for k, v in values.items() if k in spec and k in used}, space) + if named: + record["coefficients"] = named + if sub.get("initial_condition"): + record["initial"] = _initial(sub, own, ranks, values, functions, dim, where) + if sharp: + record["regularized"] = {fn: {"by": by, "width": to_code(width)} + for fn, (by, width, _) in smooth.items() if fn in sharp} + if sub.get("enforces"): + try: + parse_equation(sub["enforces"], space) + except SymbolicError as exc: + raise SymbolicError("%s, enforces: %s" % (where, exc)) from exc + carried = {(s.get("introduces") or {}).get("approximates"): sub["introduces"]["name"] + for s in [sub] if (s.get("introduces") or {}).get("approximates")} + record["boundary_conditions"] = _conditions(strong, own, carried, space, dim, + lower, upper, where) + out_subs.append(record) + + initial = {} + for ic in strong.get("initial_conditions") or []: + name = ic["field"] + if ic.get("formula") is not None: + value = parse_expression(ic["formula"], Space(dim, {}, numbers)) + else: + value = sympy.Float(ic.get("value")) + if ranks.get(name, 0) == 1 and not isinstance(value, sympy.MatrixBase): + value = sympy.Matrix([value] * dim) + if isinstance(value, sympy.MatrixBase): + initial.update({"%s_%d" % (name, i): to_code(v) for i, v in enumerate(value)}) + else: + initial[name] = to_code(value) + return {"discretization": discretization, "dim": dim, + "domain": dict({"lower": lower, "upper": upper}, + **({"region": region} if region else {})), "time": time, + "bound": {k: (v if isinstance(v, str) else to_code(v)) + for k, v in _describe_bound(entry.get("binds") or {}, entry_bound).items()}, + "initial": initial, "subproblems": out_subs} + + +def _conditions(strong, own, carried, space, dim, lower, upper, where) -> List[Dict[str, Any]]: + """The boundary conditions one sub-problem imposes, as plain data. + + Those on its own unknowns, and those on a coefficient its introduced + unknown carries (atmosphere's χ inflow is VOF's χh inflow). Conditions on + what it only reads are another sub-problem's. Each entry has ``type``, + ``variable``, ``region``, ``where`` (boolean code using ``tol``, None for + the whole boundary), ``set`` if it came from one, the ``components`` it + acts on, and its data: ``value`` (code; per component for a vector, + except slip's normal component), ``series``, robin's ``alpha``/``beta``/ + ``g``, or periodic's ``axes``. + """ + from ymf.normalize import boundary_conditions + from ymf.symbolic.problem import (_components, _given, _number, _periodic_axes, + _robin_part, _robin_rhs, _series) + from ymf.symbolic.problem import _selector + regions = {r["name"]: r for r in strong.get("boundary_regions") or []} + out = [] + try: + conditions = boundary_conditions(strong) + except ValueError as exc: + raise SymbolicError("%s: %s" % (where, exc)) from exc + for bc in conditions: + variable = bc["variable"] + if variable in own: + target, rank = variable, own[variable] + elif variable in carried: + target, rank = carried[variable], 0 + else: + continue + here = "%s, boundary condition %s on %s" % (where, bc["type"], bc["region"]) + entry: Dict[str, Any] = {"type": bc["type"], "variable": variable, + "region": bc["region"], + **_selector(regions[bc["region"]], strong["domain"], space)} + if target != variable: + entry["carried_by"] = target + if bc.get("set"): + entry["set"] = bc["set"] + entry["components"] = ([target] if rank == 0 else + ["%s_%d" % (target, i) for i in range(dim)]) + kind = bc["type"] + if kind == "periodic": + if lower is None: + raise SymbolicError("%s: periodic faces are a box's" % (here,)) + if regions[bc["region"]].get("where") is None: + raise SymbolicError("%s: periodic faces are found by position, so the " + "region needs where:, not tags:" % (here,)) + entry["axes"] = _periodic_axes(regions[bc["region"]]["where"], lower, upper, + space, bc["region"]) + domain = strong["domain"] + if isinstance(domain, dict) and domain.get("kind") == "piecewise_linear": + from ymf.symbolic.problem import periodic_faces_match + periodic_faces_match(domain, entry["axes"], lower, upper, bc["region"]) + elif kind == "robin": + entry.update(alpha=_robin_part(bc, "alpha", space), + beta=_robin_part(bc, "beta", space), g=_robin_rhs(bc, space)) + elif bc.get("series") is not None: + entry["series"] = _series(bc, target, 0 if kind == "slip" else rank) + elif bc.get("formula") is not None or bc.get("value") is not None: + value = (parse_expression(bc["formula"], space) if bc.get("formula") is not None + else _number(bc["value"])) + if kind == "slip" or rank == 0: + entry["value"] = to_code(_given(value, space, here)) + else: + if not isinstance(value, sympy.MatrixBase): + value = sympy.Matrix([value] * dim) + entry["value"] = _components(value, target, rank, space, here) + elif kind == "slip": + entry["value"] = "0.0" + out.append(entry) + return out + + +def _bind(binds: Dict[str, Any], outer: Dict[str, Any], dim: int, where: str) -> Dict[str, Any]: + """Resolve ``binds``: a symbol with only a ``meaning`` stays a symbol the + solver supplies (h, the element size); one with a formula or value is + that, in the others.""" + values = dict(outer) + formulas = {} + for name, b in binds.items(): + if b.get("value") is not None: + values[name] = sympy.Float(b["value"]) + elif b.get("formula") is not None: + formulas[name] = b["formula"] + else: + values[name] = sympy.Symbol(name, positive=True) + for _ in range(len(formulas) + 1): + for name, text in list(formulas.items()): + try: + values[name] = parse_expression(text, Space(dim, {}, values)) + except SymbolicError: + continue + del formulas[name] + if formulas: + raise SymbolicError("%s: %s could not be resolved" % (where, ", ".join(sorted(formulas)))) + return values + + +def _describe_bound(binds, values): + return {name: (b.get("meaning") if b.get("formula") is None and b.get("value") is None + else values[name]) for name, b in binds.items()} + + +def _regularizations(entry, sub, bound, dim, where): + """{sharp function: (by, width, smoother)} for this sub-problem.""" + chosen = dict(entry.get("regularize") or {}) + chosen.update(sub.get("regularize") or {}) + out = {} + for function, rule in chosen.items(): + smoothers = REGULARIZATIONS.get(function, {}) + if rule["by"] not in smoothers: + raise SymbolicError("%s: %s is regularized by %r, which ymf does not know (have: %s)" + % (where, function, rule["by"], ", ".join(smoothers) or "none")) + if rule.get("width") is None: + raise SymbolicError("%s: regularizing %s needs a width" % (where, function)) + width = parse_expression(rule["width"], Space(dim, {}, bound)) + out[function] = (rule["by"], width, smoothers[rule["by"]]) + return out + + +def _written(texts, smooth) -> set: + """The regularized sharp functions these texts use.""" + from ymf.symbolic.notation import tokens + names = set() + for text in texts: + names |= {v for k, v, _ in tokens(str(text)) if k == "NAME"} + return {fn for fn in smooth if fn in names} + + +def _initial(sub, own, ranks, values, functions, dim, where): + """An introduced unknown's initial condition: an expression in the model + unknowns at the start, which the solver reads as known fields.""" + target = sub.get("for") or sub["introduces"]["name"] + model_fields = {n: r for n, r in ranks.items() if n not in own} + space = Space(dim, {}, {k: v for k, v in values.items() if not _names_unknown(v, own)}, + model_fields, functions=functions) + try: + value = parse_expression(sub["initial_condition"], space) + except SymbolicError as exc: + raise SymbolicError("%s, initial_condition: %s" % (where, exc)) from exc + return {target: to_code(value), + "reads": sorted({f.func.__name__ for f in sympy.sympify(value).atoms( + sympy.core.function.AppliedUndef)})} + + +def _used(texts: List[str], spec: Dict[str, Any], declared: set, carried=()) -> set: + """The coefficients the texts name, directly or through other coefficients' + formulas; not through a ``carried`` one, which is read rather than + computed. A run of letters that is not a declared name counts for each + letter, as the notation reads an undeclared ``ρU`` as ``ρ·U``.""" + from ymf.symbolic.notation import tokens + + def names(text): + out = set() + for kind, value, _ in tokens(str(text)): + if kind == "NAME": + out.update([value] if value in declared or len(value) == 1 else value) + return out & set(spec) + + used, todo = set(), set().union(*[names(t) for t in texts]) if texts else set() + while todo: + name = todo.pop() + used.add(name) + if name in carried: + continue + entry = spec[name] + formula = entry.get("formula") if isinstance(entry, dict) else entry + if formula is not None: + todo |= names(formula) - used + return used + + +def _names_unknown(value, own) -> bool: + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + names = set().union(*[sympy.sympify(e).free_symbols for e in entries]) + return any(s.name in own or s.name.rsplit("_", 1)[0] in own for s in names) diff --git a/src/ymf/units.py b/src/ymf/units.py index b63559f..e5efd9b 100644 --- a/src/ymf/units.py +++ b/src/ymf/units.py @@ -94,15 +94,21 @@ def normalize_unit_string(unit_str: Optional[str]) -> Optional[str]: def to_pint_unit(unit_str: Optional[str]) -> "pint.Unit": """Parse a YMF unit string into a :class:`pint.Unit`. - ``None``, ``""``, and ``"dimensionless"`` all map to pint's dimensionless - unit. Raises :class:`UnitParseError` on anything pint can't parse even + ``None``, ``""``, ``"dimensionless"`` and a bare number such as ``"1"`` + (the usual way to write "no units") all map to pint's dimensionless + unit; only dimensions are ever compared, so a number's scale is not + kept. Raises :class:`UnitParseError` on anything pint can't parse even after :func:`normalize_unit_string`. """ if not unit_str or unit_str.strip().lower() == "dimensionless": return _UREG.dimensionless normalized = normalize_unit_string(unit_str) try: - return _UREG.parse_expression(normalized).units + parsed = _UREG.parse_expression(normalized) + if not hasattr(parsed, "units"): + # older pint (before 0.25) returns a bare number for "1" + return _UREG.dimensionless + return parsed.units except Exception as exc: # pint raises several distinct exception types raise UnitParseError( f"could not parse unit string {unit_str!r} (normalized: {normalized!r}): {exc}" @@ -295,9 +301,64 @@ def _check_against_unknown(entries, var_key, location_fmt): f"are {dimensionality(scale_units)}; state it once with " f"{{derived: true, formula: \"{name}\"}}")) + _check_state_coefficients_in_scales(problem, strong_form, result) return result +def _check_state_coefficients_in_scales(problem, strong_form, result) -> None: + """A scale or a dimensionless number named after a coefficient that is a + function of the state. + + A coefficient may name the unknowns (a density ``ρw (1 - χ) + ρa χ``). + It then has no single value, so a ``density: {derived: true, formula: + "ρ"}`` scale, or an ``Re`` built on it, cannot be evaluated and would + drop out of ``non_dimensionalize`` without a word. Say which, and how to + fix it: derive it from a constant coefficient instead. + """ + try: + from ymf.symbolic.notation import tokens + from ymf.symbolic.problem import state_coefficients + except ImportError: # no sympy: no formulas are read at all + return + try: + dim = _space_dimension(strong_form.get("domain")) + except _Undetermined: + dim = 3 + try: + state = state_coefficients(strong_form, dim) + except Exception: # an unreadable spec is reported elsewhere + return + if not state: + return + formulas = [] + for name, spec in (problem.get("characteristic_scales") or {}).items(): + if isinstance(spec, dict) and spec.get("formula"): + formulas.append((f"characteristic_scales[{name}]", spec["formula"])) + for name, spec in (problem.get("dimensionless_numbers") or {}).items(): + if isinstance(spec, dict) and spec.get("formula"): + formulas.append((f"dimensionless_numbers[{name}]", spec["formula"])) + for location, formula in formulas: + try: + words = [value for kind, value, _ in tokens(str(formula)) if kind == "NAME"] + except Exception: + continue + named = [] + for word in words: + # "ρU" is ρ·U when it is not itself a name, as in the notation + parts = [word] if word in state or len(word) == 1 else list(word) + named.extend(p for p in parts if p in state and p not in named) + for coefficient in named: + result.issues.append(UnitIssue( + "warning", location, + f"coefficients[{coefficient}] depends on the unknowns " + f"({', '.join(state[coefficient])}), so this has no value; derive it " + f"from a constant coefficient instead")) + + +#: conditions whose g is a value of the variable itself, like dirichlet's +_VALUE_KINDS = ("slip", "inflow", "characteristic") + + #: What :func:`check` does with a boundary condition it cannot check for #: want of declared units: say nothing, warn, or call it an error. MISSING_UNITS = ("skip", "warn", "error") @@ -437,7 +498,8 @@ def _check_boundary_conditions(strong_form: Dict[str, Any], result: "UnitCheckRe ``units`` on a condition is the units of g, its right-hand side. - - dirichlet, u = g: ``[g] = [u]``. + - dirichlet, u = g: ``[g] = [u]``; so too slip (the normal component + of u), inflow and characteristic. open has no g. - neumann, q.n = g, where q = -a grad u is the diffusive flux of the variable's own equation: ``[g] = [a][u]/[L]``. - robin, alpha u + beta q.n = g: ``[alpha][u] = [beta][q.n] = [g]``. @@ -451,8 +513,13 @@ def _check_boundary_conditions(strong_form: Dict[str, Any], result: "UnitCheckRe ymf[symbolic] to read alpha, beta and the equations; without it they are undecided too. """ - bcs = [bc for bc in (strong_form.get("boundary_conditions") or []) - if bc.get("type") in ("dirichlet", "neumann", "robin")] + from ymf.normalize import boundary_conditions + try: + conditions = boundary_conditions(strong_form) + except ValueError: + return # an undefined set: reported on load + bcs = [bc for bc in conditions + if bc.get("type") in ("dirichlet", "neumann", "robin") + _VALUE_KINDS] if not bcs: return dims = _Dimensions() @@ -483,6 +550,8 @@ def report(location, verdict): for bc in bcs: var, kind = bc.get("variable"), bc["type"] location = f"boundary_conditions[region={bc.get('region')!r}, variable={var!r}]" + if bc.get("set"): + location += f" (set {bc['set']})" if var not in unknown_dims: continue try: @@ -491,7 +560,7 @@ def report(location, verdict): result.issues.append(UnitIssue("error", location, str(exc))) continue u = unknown_dims[var] - if kind == "dirichlet": + if kind == "dirichlet" or kind in _VALUE_KINDS: verdict = dims.compare(g, u) if verdict and verdict[0] == "error": declared = _unknown_units(unknowns, var) @@ -629,9 +698,10 @@ def _flux_dimension(strong_form, component, names, dims, space, cache): if "equations" not in cache: residuals = [] try: - for text in strong_form.get("equations") or []: + from ymf.normalize import equations_in_order + for text in equations_in_order(strong_form): residuals.extend(parse_equation(text, space)) - except SymbolicError as exc: + except (SymbolicError, ValueError) as exc: raise _Undetermined(["the equations do not read: %s" % exc]) cache["equations"] = residuals residuals = cache["equations"] diff --git a/tests/fixtures/shallow_water/lake_at_rest.ymf b/tests/fixtures/shallow_water/lake_at_rest.ymf new file mode 100644 index 0000000..cdfabd6 --- /dev/null +++ b/tests/fixtures/shallow_water/lake_at_rest.ymf @@ -0,0 +1,36 @@ +# Well-balancing (C35's first row): a lake at rest over a parabolic bump +# stays at rest, h + Z = 0.5 and q = 0, to round-off. Exact. +extends: proteus_numerics.ymf + +model: + name: "Lake at rest over a bump, 25 m × 1 m" + physical_model: {provenance: human_specified} + strong_form: + domain: "Ω = [0, 25] × [0, 1], t ∈ (0, 10]" + boundary_regions: + - {name: all, where: "∂Ω"} + boundary_sets: + - {region: all, set: wall} + initial_conditions: + - {field: h, type: function, formula: "0.5 − Z", units: m} + - {field: q, type: constant, value: 0.0} + coefficients: + g: {value: 9.81} + n: {value: 0.01} + θ: {value: 1.0} + Z: {formula: "max(0, 0.2 − 0.05 (x − 10)²)"} + α: {value: 0.0} + β: {value: 0.0} + +solution_paths: + analytical: + - name: at_rest + provenance: human_specified + method: closed_form + solution: + formula: |- + h = 0.5 − Z + q = (0, 0) + discretizations: + - name: convex_limited + mesh: {cells: [50, 100, 200]} diff --git a/tests/fixtures/shallow_water/proteus_numerics.ymf b/tests/fixtures/shallow_water/proteus_numerics.ymf new file mode 100644 index 0000000..e21c7f7 --- /dev/null +++ b/tests/fixtures/shallow_water/proteus_numerics.ymf @@ -0,0 +1,31 @@ +# How Proteus solves the shallow-water model: SW2DCV, the convex-limited +# graph-viscosity scheme with SSP RK3, as one sub-problem. Shared by the +# cases; each adds its mesh and may change the free-form scheme options. +extends: shallow_water.ymf +ymf_state: abstract + +solution_paths: + analytical: [] + discretizations: + - name: convex_limited + provenance: proteus_derived + binds: + he: {meaning: element_size} # h is the depth here + regularize: # a canopy's edge, over one element + heaviside: {by: smoothed_heaviside, width: he} + step_control: {type: min_cfl, cfl: 0.33} + solver: {type: time_marching, linear_solver: petsc} + # the scheme is algebraic (no weak form); its options free-form until + # a second solver needs them typed + discretization: + scheme: graph_viscosity_convex_limiting + time_integration: SSP_RK3 + entropy_viscosity_cE: 1.0 + lumped_mass: false + bounds_relaxation: true + subproblems: + - name: swe + solves: [mass, momentum] + finite_element: {fields: {family: CG, order: 1}} + time_integration: SSP_RK3 + proteus: {model: SW2DCV, INTEGRATE_BY_PARTS_WALLS: 1} diff --git a/tests/fixtures/shallow_water/shallow_water.ymf b/tests/fixtures/shallow_water/shallow_water.ymf new file mode 100644 index 0000000..badbd40 --- /dev/null +++ b/tests/fixtures/shallow_water/shallow_water.ymf @@ -0,0 +1,42 @@ +# The shallow-water equations with porosity, bed friction and canopy drag +# (C42's second family, C35's model): the model only. Unknowns are the depth +# h and the discharge q = h u. The pressure is written as a flux so the +# system is one ADR system; Proteus's SW2DCV discretizes the well-balanced +# non-conservative form θ g h ∇(h + Z), which is the same strong form. +# Friction and drag are coefficients of the state; a dry node (h = 0) makes +# q/h singular, which is the discretization's to regularize (SW2DCV's hEps). +ymf_state: abstract + +model: + physical_model: + provenance: human_specified + processes: [free_surface_flow, bed_friction, canopy_drag] + assumptions: [hydrostatic, depth_averaged, nondispersive] + units: {convention: SI, default: true} + strong_form: + provenance: human_specified + unknowns: + - {name: h, units: m} + - {name: q, units: m2/s, rank: 1} + equation_formulation: "shallow-water equations, conservative variables, porous medium" + equations: + - name: mass + for: h + equation: "∂(θ h)/∂t + ∇·(θ q) = 0 in Ω" + - name: momentum + for: q + equation: "∂(θ q)/∂t + ∇·(θ q⊗q/h) + ∇(θ g h²/2) − (g h²/2)∇θ + θ g h ∇Z + θ τb + Fd = 0 in Ω" + coefficients: + τb: {formula: "g n² norm(q) q / h^(7/3)"} # Manning bed friction + Fd: {formula: "(α + β norm(q)/h) q"} # Darcy-Forchheimer drag, per unit plan area + g: {units: m/s2} + n: {units: s/m^(1/3)} + θ: {units: "1"} # porosity: 1 - m π D²/4 in a canopy, 1 elsewhere + Z: {units: m} # the bed + α: {units: 1/s} + β: {units: 1/m} + condition_sets: + wall: # a no-flux wall (SW2DCV's integrated-by-parts walls) + - {variable: q, type: slip} + outflow: # the natural outflow (an open end) + - {variable: q, type: open} diff --git a/tests/fixtures/shallow_water/torino.ymf b/tests/fixtures/shallow_water/torino.ymf new file mode 100644 index 0000000..8e6cecc --- /dev/null +++ b/tests/fixtures/shallow_water/torino.ymf @@ -0,0 +1,59 @@ +# The Torino flume (torino_flume/dswe/torino.py, sw_model 0): a dam break +# down a 3 % slope into an emergent rod canopy, the flume open at its far +# end. The deck's options carried over: dam_height 0.15 at x 5.173 along the +# bed, the canopy from 5.2255 − 0.5·0.0525 to 5.2255 + 75.5·0.0525 along the +# bed, 194 rods/m² of 6 mm, C_D under the Froude law, cE = 5, refinement 5 +# (he = 0.0725 m). Not carried: the 'age' drag law (it reads each node's +# time since wetting, a history no coefficient holds), the gate snapped to a +# node with half the depth there (a mesh-dependent initial condition), the +# structured 161 × 8 mesh (a mesh gives one size, not one per axis), and +# the load-cell forces (a derived-quantity gauge, C34). +extends: proteus_numerics.ymf + +model: + name: "Torino flume, S = 0.03, H = 0.15 m, canopy" + physical_model: {provenance: human_specified} + strong_form: + domain: "Ω = [0, Lx] × [0, 0.51], t ∈ (0, 15]" + boundary_regions: + - {name: walls, where: "x = 0 or y = 0 or y = 0.51"} + - {name: open_end, where: "x = Lx"} + boundary_sets: + - {region: walls, set: wall} + - {region: open_end, set: outflow} + initial_conditions: + - field: h + type: function + # the reservoir upstream of the gate, its level clipped at the bed + formula: "max(Z0 + H − Z, 0) heaviside(xg − x)" + units: m + - {field: q, type: constant, value: 0.0} + coefficients: + cS: {formula: "1/sqrt(1 + S²)"} # along-bed to horizontal + Lx: {formula: "11.6 cS"} + S: {value: 0.03} + Z: {formula: "−S x"} + H: {value: 0.15, units: m} + xg: {formula: "5.173 cS"} # the gate + Z0: {formula: "−S xg"} + xa: {formula: "(5.2255 − 0.5·0.0525) cS"} # the canopy + xb: {formula: "(5.2255 + 75.5·0.0525) cS"} + canopy: {formula: "heaviside(x − xa) heaviside(xb − x)"} + m: {value: 194.0, units: 1/m2} + D: {value: 0.006, units: m} + θ: {formula: "1 − canopy m π D²/4"} + CD: {value: 1.0} + frc: {value: 0.9} + frp: {value: 0.8} + Fr: {formula: "norm(q) / (h sqrt(g h))"} + β: {formula: "canopy CD min(1, frc Fr^(−frp)) m D / 2"} + α: {value: 0.0} + g: {value: 9.81} + n: {value: 0.0} + +solution_paths: + discretizations: + - name: convex_limited + mesh: {h: 0.0725} + discretization: + entropy_viscosity_cE: 5.0 diff --git a/tests/fixtures/two_phase/dambreak.ymf b/tests/fixtures/two_phase/dambreak.ymf new file mode 100644 index 0000000..25596d9 --- /dev/null +++ b/tests/fixtures/two_phase/dambreak.ymf @@ -0,0 +1,102 @@ +# The 2-D dam break of test/TwoPhaseFlow/damBreakLevelSet.py, split into +# proteus's five sub-problems (RANS2P, VOF, NCLS, RDLS, MCorr). The stored +# mesh (meshDambreak) waits on C7.1; this gives an element size instead. +extends: two_phase_level_set.ymf + +model: + name: "Dam break, 3.22 m × 1.8 m tank, 1.2 m × 0.6 m column (2-D)" + physical_model: {provenance: human_specified} + strong_form: + domain: "Ω = [0, 3.22] × [0, 1.8], t ∈ (0, 3]" + boundary_regions: + - {name: walls, where: "x = 0 or x = 3.22 or y = 0"} + - {name: top, where: "y = 1.8"} + boundary_sets: + - {region: walls, set: free_slip} + - {region: top, set: atmosphere} + initial_conditions: + - {field: v, type: constant, value: 0.0} + - {field: p, type: constant, value: 0.0} + - field: φ + type: function + formula: "sqrt(max(x − 1.2, 0)² + max(y − 0.6, 0)²) + min(max(x − 1.2, y − 0.6), 0)" + units: m + coefficients: + ρw: {value: 998.2} + νw: {value: 1.004e-6} + ρa: {value: 1.205} + νa: {value: 1.5e-5} + g: {formula: "(0, -9.8)"} + + # the deck's gauges: the level set along four verticals and one + # horizontal (height1.csv, height2.csv), the pressure at three points on + # the right wall (pressure.csv) + gauges: + - name: heights + kind: line + fields: [φ] + lines: [[[2.724, 0.0], [2.724, 1.8]], [[2.228, 0.0], [2.228, 1.8]], + [[1.732, 0.0], [1.732, 1.8]], [[0.582, 0.0], [0.582, 1.8]]] + file: height1.csv + - name: bottom + kind: line + fields: [φ] + lines: [[[0.0, 0.0], [3.22, 0.0]]] + file: height2.csv + - name: pressure + kind: point + fields: [p] + points: [[3.22, 0.16], [3.22, 0.584], [3.22, 0.12]] + file: pressure.csv + +solution_paths: + analytical: [] + discretizations: + - name: level_set_five + provenance: proteus_derived + mesh: {h: 0.0322} + binds: + h: {meaning: element_size} + ε: {formula: "1.5 h"} + regularize: + heaviside: {by: smoothed_heaviside, width: ε} + step_control: {type: min_cfl, cfl: 0.25, step_exact: true} + solver: {type: nonlinear, nonlinear_solver: newton, linear_solver: petsc, max_iterations: 50} + subproblems: + - name: flow + solves: [momentum, continuity] + lags: {φ: previous_step, χ: {via: vof, at: previous_step}} + finite_element: {velocity: {family: CG, order: 1}, pressure: {family: CG, order: 1}} + time_integration: backward_euler + solver: {max_iterations: 100} + proteus: {model: RANS2P, useVF: 1.0} + - name: vof + introduces: {name: χh, approximates: χ, units: "1"} + equation: "∂χh/∂t + ∇·(v χh) = 0" + initial_condition: "heaviside(φ)" + lags: {v: {via: flow, at: current_step}} + finite_element: {fields: {family: CG, order: 1}} + time_integration: backward_euler + - name: ncls + solves: [interface] + lags: {v: {via: flow, at: current_step}} + finite_element: {fields: {family: CG, order: 1}} + time_integration: backward_euler + - name: rdls + enforces: "norm(∇φ) = 1" + for: φ + equation: "sign(φs) (norm(∇φ) − 1) = 0" + lags: {φs: {of: φ, via: ncls, at: current_step}} + regularize: {sign: {by: smoothed_sign, width: εr}} + binds: {εr: {formula: "0.75 h"}} + pseudo_time: steady + writes_back: {φ: replace} + finite_element: {fields: {family: CG, order: 1}} + - name: mcorr + enforces: "heaviside(φ) = χh" + introduces: {name: φc, correction_of: φ, units: m} + equation: "heaviside(φ + φc) − χh − ∇·(εd ∇φc) = 0" + binds: {εd: {formula: "10 h"}} + lags: {φ: {via: rdls, at: current_step}, χh: {via: vof, at: current_step}} + writes_back: {φ: add} + finite_element: {fields: {family: CG, order: 1}} diff --git a/tests/fixtures/two_phase/two_phase_level_set.ymf b/tests/fixtures/two_phase/two_phase_level_set.ymf new file mode 100644 index 0000000..1036d7b --- /dev/null +++ b/tests/fixtures/two_phase/two_phase_level_set.ymf @@ -0,0 +1,48 @@ +# Incompressible two-phase flow with a sharp interface carried by a level +# set (C42): the model only. Phase 0 is water (φ < 0), phase 1 air (φ > 0). +# The density and viscosity are coefficients of the state, through +# χ = heaviside(φ). The whole model is not one ADR system (v·∇φ is a +# Hamiltonian in another unknown); a split discretization classifies each +# sub-problem with the others' fields lagged. +ymf_state: abstract + +model: + physical_model: + provenance: human_specified + processes: [incompressible_flow, two_phase_flow, interface_transport, gravity] + assumptions: [newtonian, immiscible, isothermal, no_phase_change, no_surface_tension] + strong_form: + provenance: human_specified + unknowns: + - {name: v, units: m/s, rank: 1} + - {name: p, units: Pa} + - {name: φ, units: m} + equation_formulation: "variable-density incompressible Navier-Stokes, interface by a level set" + equations: + - name: momentum + for: v + equation: "∂(ρ v)/∂t + ∇·(ρ v⊗v) − ∇·(2μ sym(∇v)) + ∇p = ρ g in Ω" + - name: continuity + for: p + equation: "∇·v = 0 in Ω" + - name: interface + for: φ + equation: "∂φ/∂t + v·∇φ = 0 in Ω" + coefficients: + χ: {formula: "heaviside(φ)", units: "1"} + ρ: {formula: "ρw (1 − χ) + ρa χ", units: kg/m3} + μ: {formula: "ρw νw (1 − χ) + ρa νa χ", units: Pa*s} + ρw: {units: kg/m3} + νw: {units: m2/s} + ρa: {units: kg/m3} + νa: {units: m2/s} + g: {units: m/s2} + # what the cases' walls and tops mean, written once (Proteus's + # BC_RANS.setFreeSlip and setAtmosphere, for this model's unknowns) + condition_sets: + free_slip: + - {variable: v, type: slip} + atmosphere: + - {variable: p, type: dirichlet, value: 0.0} + - {variable: v, type: open} + - {variable: χ, type: inflow, value: 1.0} # what enters is air diff --git a/tests/test_boundary_kinds.py b/tests/test_boundary_kinds.py new file mode 100644 index 0000000..f92752f --- /dev/null +++ b/tests/test_boundary_kinds.py @@ -0,0 +1,101 @@ +"""slip, open, inflow and characteristic; conditions on a coefficient. + +The kinds a split discretization's solver needs (C42): a free-slip wall, +an open boundary to the air, the phase entering through it, a tailwater +stage. They load and are checked; one ADR system runs none of them, and +adr_problem says which solver does. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + + +@pytest.fixture +def heat(): + return load_ymf(EXAMPLES / "heat_equation.ymf") + + +@pytest.fixture +def poiseuille(): + return load_ymf(EXAMPLES / "navier_stokes" / "plane_poiseuille.ymf") + + +def _with(doc, *conditions): + doc = copy.deepcopy(doc) + doc["model"]["strong_form"]["boundary_conditions"] = list(conditions) + return doc + + +def test_the_four_kinds_load(tmp_path): + text = (EXAMPLES / "navier_stokes" / "plane_poiseuille.ymf").read_text(encoding="utf-8") + start = text.index(" boundary_conditions:") + end = text.index(" coefficients:", start) + (tmp_path / "navier_stokes_model.ymf").write_text( + (EXAMPLES / "navier_stokes" / "navier_stokes_model.ymf").read_text(encoding="utf-8"), + encoding="utf-8") + (tmp_path / "p.ymf").write_text(text[:start] + """ boundary_conditions: + - {region: walls, variable: v, type: slip} + - {region: ends, variable: v, type: open} + - {region: ends, variable: p, type: inflow, value: 0.0} + - {region: pressure_datum, variable: p, type: characteristic, series: [[0, 0], [1, 1]]} +""" + text[end:], encoding="utf-8") + kinds = [bc["type"] for bc in load_ymf(tmp_path / "p.ymf")["model"]["strong_form"][ + "boundary_conditions"]] + assert kinds == ["slip", "open", "inflow", "characteristic"] + + +@pytest.mark.parametrize("condition, message", [ + ({"region": "walls", "variable": "T", "type": "slip"}, + "slip gives the normal component of a vector, and T is a scalar"), + ({"region": "walls", "variable": "T", "type": "inflow"}, + "an inflow condition needs a value, formula or series"), + ({"region": "walls", "variable": "T", "type": "characteristic"}, + "a characteristic condition needs a value"), + ({"region": "walls", "variable": "T", "type": "open", "value": 1.0}, + "an open condition takes no value"), + ({"region": "walls", "variable": "q", "type": "dirichlet", "value": 1.0}, + "q is neither an unknown nor a coefficient"), + ({"region": "lid", "variable": "T", "type": "dirichlet", "value": 1.0}, + "no boundary region named 'lid'"), +]) +def test_what_a_condition_needs_is_checked(heat, condition, message): + with pytest.raises(YmfReferenceError, match=message): + resolve(_with(heat, condition)) + + +def test_a_condition_may_name_a_coefficient(heat): + resolve(_with(heat, {"region": "walls", "variable": "κ", "type": "inflow", "value": 1.0})) + + +def test_slip_on_a_vector_needs_no_value(poiseuille): + resolve(_with(poiseuille, {"region": "walls", "variable": "v", "type": "slip"})) + + +def test_one_adr_system_runs_none_of_them(heat, poiseuille): + pytest.importorskip("sympy") + from ymf.symbolic import SymbolicError, adr_problem + with pytest.raises(SymbolicError, match="slip conditions are for a split discretization"): + adr_problem(_with(poiseuille, {"region": "walls", "variable": "v", "type": "slip"})) + with pytest.raises(SymbolicError, match="a condition on a coefficient is for a split"): + adr_problem(_with(heat, {"region": "walls", "variable": "κ", "type": "dirichlet", + "value": 1.0})) + + +def test_the_units_of_an_inflow_value_are_those_of_its_variable(heat): + pytest.importorskip("pint") + from ymf import check_units + wrong = _with(heat, {"region": "walls", "variable": "T", "type": "inflow", + "value": 300.0, "units": "m"}) + assert not check_units(wrong) + right = _with(heat, {"region": "walls", "variable": "T", "type": "inflow", + "value": 300.0, "units": "K"}) + assert check_units(right) diff --git a/tests/test_condition_sets.py b/tests/test_condition_sets.py new file mode 100644 index 0000000..7c63110 --- /dev/null +++ b/tests/test_condition_sets.py @@ -0,0 +1,97 @@ +"""Condition sets: written once in a model, applied to regions by a problem. + +The 25 air-water-vv two-phase cases are written in proteus's boundary- +condition setters (setFreeSlip, setAtmosphere, ...), each a bundle of +conditions on every model's variables. A ymf model names such bundles over +its own unknowns and coefficients; a case says where each holds. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.normalize import boundary_conditions # noqa: E402 +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 + +TWO_PHASE = Path(__file__).parent / "fixtures" / "two_phase" +EXAMPLES = Path(__file__).parent.parent / "examples" + + +@pytest.fixture +def dambreak(): + return load_ymf(TWO_PHASE / "dambreak.ymf") + + +def _strong(doc): + return doc["model"]["strong_form"] + + +def test_the_dam_break_applies_free_slip_and_atmosphere(dambreak): + assert [(bc["region"], bc["variable"], bc["type"], bc["set"]) + for bc in boundary_conditions(_strong(dambreak))] == [ + ("walls", "v", "slip", "free_slip"), + ("top", "p", "dirichlet", "atmosphere"), + ("top", "v", "open", "atmosphere"), + ("top", "χ", "inflow", "atmosphere"), + ] + + +def test_a_condition_written_out_replaces_the_sets_member(dambreak): + strong = _strong(dambreak) + strong["boundary_conditions"] = [ + {"region": "top", "variable": "p", "type": "dirichlet", "formula": "9.8 t"}] + conditions = boundary_conditions(strong) + pressure = [bc for bc in conditions if bc["region"] == "top" and bc["variable"] == "p"] + assert pressure == [{"region": "top", "variable": "p", "type": "dirichlet", + "formula": "9.8 t"}] + assert len(conditions) == 4 + + +def test_an_undefined_set_is_refused(dambreak): + _strong(dambreak)["boundary_sets"].append({"region": "walls", "set": "no_slip"}) + with pytest.raises(YmfReferenceError, match="no condition set named 'no_slip'"): + resolve(dambreak) + + +def test_two_sets_may_not_both_decide_one_regions_variable(dambreak): + strong = _strong(dambreak) + strong["condition_sets"]["no_slip"] = [{"variable": "v", "type": "dirichlet", "value": 0.0}] + strong["boundary_sets"].append({"region": "walls", "set": "no_slip"}) + with pytest.raises(YmfReferenceError, match="free_slip and no_slip both give v"): + resolve(dambreak) + + +def test_a_sets_members_are_checked_where_they_are_defined(dambreak): + strong = _strong(dambreak) + strong["condition_sets"]["broken"] = [{"variable": "p", "type": "slip"}] + with pytest.raises(YmfReferenceError, + match=r"condition_sets\.broken\[variable='p'\]: slip .* p is a scalar"): + resolve(dambreak) + + +def test_a_set_is_what_its_conditions_would_be(tmp_path): + """The heat equation's wall condition as a set gives the same problem.""" + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + plain = load_ymf(EXAMPLES / "heat_equation.ymf") + with_set = copy.deepcopy(plain) + strong = _strong(with_set) + (wall,) = strong.pop("boundary_conditions") + strong["condition_sets"] = {"held": [{k: v for k, v in wall.items() if k != "region"}]} + strong["boundary_sets"] = [{"region": wall["region"], "set": "held"}] + a, b = adr_problem(plain), adr_problem(with_set) + assert a["dirichlet"] == b["dirichlet"] and a["adr"] == b["adr"] + + +def test_the_units_check_says_which_set_a_condition_came_from(dambreak): + pytest.importorskip("pint") + from ymf import check_units + strong = _strong(dambreak) + strong["condition_sets"]["atmosphere"][0] = {"variable": "p", "type": "dirichlet", + "value": 0.0, "units": "m"} + issues = [i for i in check_units(dambreak).issues if i.severity == "error"] + assert any("(set atmosphere)" in i.location for i in issues) diff --git a/tests/test_discrete_domains.py b/tests/test_discrete_domains.py new file mode 100644 index 0000000..4023dd4 --- /dev/null +++ b/tests/test_discrete_domains.py @@ -0,0 +1,102 @@ +"""`discrete` domains: a mesh as ymf stores one, solved as it stands (C27 +Phase 5, D4; C7.1: a mesh is a type of domain). + +Its bounds are its nodes' when the document holds them, and the solver's to +take when they are in a file; its boundary regions select the tags of its +facets (#6); a discretization gives only its levels; it is not periodic. +The solver reads it through Proteus's MeshYmfDomain (cekees/proteus#186). +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") +pytest.importorskip("sympy") + +from ymf.domain import validate_domain_block # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 +from ymf.state import reached # noqa: E402 +from ymf.symbolic import SymbolicError, adr_problem # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + + +def _item(values, dims, data_type="Float"): + return {"data_type": data_type, "precision": 8 if data_type == "Float" else 4, + "dimensions": dims, "values": values} + + +#: the square [0, 2] x [0, 1] as two triangles, its four edges tagged +SQUARE = { + "kind": "discrete", + "topology": {"type": "Triangle", "number_of_elements": 2, + "data_item": _item([0, 1, 2, 0, 2, 3], [2, 3], "Int")}, + "geometry": {"type": "XY", "data_item": _item([0, 0, 2, 0, 2, 1, 0, 1], [4, 2])}, + "facets": {"topology": {"type": "Polyline", "number_of_elements": 4, + "data_item": _item([0, 1, 1, 2, 2, 3, 3, 0], [4, 2], "Int")}, + "tags": {"data_item": _item([1, 2, 3, 4], [4], "Int")}, + "view": "global"}, + "boundary_tags": {"bottom": 1, "right": 2, "top": 3, "left": 4}, +} + + +def _on(domain, mesh): + doc = load_ymf(EXAMPLES / "poisson.ymf") + strong = doc["model"]["strong_form"] + strong["domain"] = copy.deepcopy(domain) + strong["boundary_regions"] = [{"name": "fixed", "tags": ["bottom", "left"]}, + {"name": "rest", "tags": [2, 3]}] + strong["boundary_conditions"] = [ + {"region": "fixed", "variable": "u", "type": "dirichlet", "value": 0.0}, + {"region": "rest", "variable": "u", "type": "neumann", "value": 0.0}] + doc["solution_paths"]["analytical"] = [] + doc.pop("vvuq", None) + for d in doc["solution_paths"]["discretizations"]: + d["mesh"] = copy.deepcopy(mesh) + return doc + + +def test_a_mesh_is_a_domain(): + validate_domain_block(SQUARE) + + +def test_its_bounds_are_its_nodes(): + problem = adr_problem(_on(SQUARE, {"levels": 1})) + assert problem["domain"]["lower"] == [0.0, 0.0] + assert problem["domain"]["upper"] == [2.0, 1.0] + assert problem["domain"]["region"]["kind"] == "discrete" + assert problem["dirichlet"][0]["flags"] == [1, 4] + assert problem["neumann"][0]["flags"] == [2, 3] + + +def test_nodes_in_a_file_leave_the_bounds_to_the_solver(): + square = copy.deepcopy(SQUARE) + square["geometry"]["data_item"] = {"format": "HDF", "data_type": "Float", + "precision": 8, "dimensions": [4, 2], + "data": "square.h5:/nodes"} + problem = adr_problem(_on(square, {"levels": 1})) + assert problem["domain"]["lower"] is None + assert problem["adr"]["dim"] == 2 # from the geometry's type + square["bounds"] = {"lower": [0, 0], "upper": [2, 1]} + assert adr_problem(_on(square, {"levels": 1}))["domain"]["upper"] == [2.0, 1.0] + + +def test_a_discretization_refines_it_by_levels(): + from ymf.closure import planned_outputs + doc = _on(SQUARE, {"levels": [1, 2]}) + assert reached(doc)[0] == "discretized" + assert sorted(k.rsplit("/", 1)[0] for k, _, _ in planned_outputs(doc))[:2] == [ + "P1_linear/levels=1", "P1_linear/levels=2"] + assert "already a mesh" in reached(_on(SQUARE, {"cells": 4}))[1] + + +def test_a_discrete_domain_is_not_periodic(): + doc = _on(SQUARE, {"levels": 1}) + strong = doc["model"]["strong_form"] + strong["boundary_regions"].append({"name": "ends", "where": "x = 0 or x = 2"}) + strong["boundary_conditions"].append({"region": "ends", "variable": "u", + "type": "periodic"}) + with pytest.raises(SymbolicError, match="a discrete domain cannot be periodic"): + adr_problem(doc) diff --git a/tests/test_domain.py b/tests/test_domain.py index e7ddc9a..45f25e5 100644 --- a/tests/test_domain.py +++ b/tests/test_domain.py @@ -162,16 +162,25 @@ def test_the_two_box_forms_give_the_same_problem(): assert adr_problem(doc) == adr_problem(mapped) -def test_a_kind_nothing_solves_yet_says_so(tmp_path): - """Validating and carrying a domain is not the same as running it.""" +def test_a_file_domain_is_posed_from_its_bounds_or_its_dimension(tmp_path): + """A file domain is carried whole for the solver to read (C27 Phase 5, + D2). Its bounds set the boundary tolerance's scale; without them the + solver takes them from the mesh, but the dimension must still be known + to pose the problem.""" pytest.importorskip("sympy") from ymf.symbolic.language import SymbolicError from ymf.symbolic.problem import adr_problem - doc = load(tmp_path, '{kind: file, format: "stl", path: "hull.stl"}') + doc = load(tmp_path, '{kind: file, format: "tetgen", path: "cube"}') doc["model"]["strong_form"]["equations"] = ["−Δu = 0"] - with pytest.raises(SymbolicError, match="adr_problem can run a box"): + with pytest.raises(SymbolicError, match="gives its bounds, or at least its dimension"): adr_problem(doc) + doc["model"]["strong_form"]["domain"]["dimension"] = 2 + problem = adr_problem(doc) + assert problem["domain"]["lower"] is None and problem["domain"]["upper"] is None + assert problem["domain"]["region"]["format"] == "tetgen" + doc["model"]["strong_form"]["domain"]["bounds"] = {"lower": [0, 0], "upper": [2, 1]} + assert adr_problem(doc)["domain"]["upper"] == [2.0, 1.0] # -------------------------------------------------------------------------- diff --git a/tests/test_file_domains.py b/tests/test_file_domains.py new file mode 100644 index 0000000..6b74f9d --- /dev/null +++ b/tests/test_file_domains.py @@ -0,0 +1,121 @@ +"""`file` domains (C27 Phase 5, D2 and D3). + +A file domain is carried whole for the solver: a mesh file (2dm, 3dm, hex, +tetgen, ymf) fixes the mesh, so a discretization gives only its `levels` of +uniform refinement; a geometry (gmsh_geo) is meshed by `h`. Its `bounds` are +optional: without them the solver takes them from the mesh. Boundary +regions select its boundary tags; periodic conditions are not carried. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.domain import YmfDomainError, validate_domain_block # noqa: E402 +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 +from ymf.state import reached # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + +CUBE = {"kind": "file", "format": "tetgen", "path": "cube", "dimension": 3, + "boundary_tags": {"bottom": 1, "top": 2, "front": 3, "back": 4, "left": 5, + "right": 6}} + + +def _on_file(domain, mesh): + doc = load_ymf(EXAMPLES / "poisson.ymf") + strong = doc["model"]["strong_form"] + strong["domain"] = copy.deepcopy(domain) + strong["boundary_regions"] = [{"name": "fixed", "tags": ["bottom", "top"]}, + {"name": "sides", "tags": [3, 4, 5, 6]}] + strong["boundary_conditions"] = [ + {"region": "fixed", "variable": "u", "type": "dirichlet", "value": 0.0}, + {"region": "sides", "variable": "u", "type": "neumann", "value": 0.0}] + doc["solution_paths"]["analytical"] = [] + doc.pop("vvuq", None) + for d in doc["solution_paths"]["discretizations"]: + d["mesh"] = copy.deepcopy(mesh) + return doc + + +def test_a_mesh_file_plans_its_levels(): + from ymf.closure import planned_outputs + doc = _on_file(CUBE, {"levels": [1, 2]}) + assert reached(doc)[0] == "discretized" + keys = sorted(k for k, _, _ in planned_outputs(doc)) + assert [k.rsplit("/", 1)[0] for k in keys] == [ + "P1_linear/levels=1", "P1_linear/levels=2", "P2_quadratic/levels=1", + "P2_quadratic/levels=2"] + + +def test_a_mesh_file_takes_no_mesh_size(): + state, why = reached(_on_file(CUBE, {"h": 0.1})) + assert state == "well-posed" + assert "already a mesh; give mesh: {levels: ...} only" in why + + +def test_a_geometry_file_is_meshed_by_h(): + from ymf.closure import planned_outputs + geo = dict(CUBE, format="gmsh_geo", path="cube.geo") + keys = [k for k, _, _ in planned_outputs(_on_file(geo, {"h": [0.2, 0.1]}))] + assert keys[0].startswith("P1_linear/h=0.2/") + + +def test_the_tags_resolve_through_the_files_boundary_tags(): + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + problem = adr_problem(_on_file(CUBE, {"levels": 1})) + assert problem["domain"] == {"lower": None, "upper": None, + "region": problem["domain"]["region"]} + assert problem["dirichlet"][0]["flags"] == [1, 2] + assert problem["neumann"][0]["flags"] == [3, 4, 5, 6] + + +def test_bounds_set_the_box(): + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + cube = dict(CUBE, bounds={"lower": [0, 0, 0], "upper": [1, 2, 3]}) + domain = adr_problem(_on_file(cube, {"levels": 1}))["domain"] + assert domain["lower"] == [0.0, 0.0, 0.0] and domain["upper"] == [1.0, 2.0, 3.0] + + +@pytest.mark.parametrize("bounds, message", [ + ({"lower": [0, 0], "upper": [1, 1, 1]}, "same, nonzero, number"), + ({"lower": [0, 0, 1], "upper": [1, 1, 1]}, r"upper\[2\] = 1 is not above"), + ({"lower": [0, 0], "upper": [1, 1]}, "dimension 3 but bounds of 2"), +]) +def test_bounds_are_checked(bounds, message): + with pytest.raises(YmfDomainError, match=message): + validate_domain_block(dict(CUBE, bounds=bounds)) + + +def test_a_file_domain_is_not_periodic(): + pytest.importorskip("sympy") + from ymf.symbolic import SymbolicError, adr_problem + doc = _on_file(dict(CUBE, bounds={"lower": [0, 0, 0], "upper": [1, 1, 1]}), {"levels": 1}) + strong = doc["model"]["strong_form"] + strong["boundary_regions"].append({"name": "ends", "where": "x = 0 or x = 1"}) + strong["boundary_conditions"].append({"region": "ends", "variable": "u", + "type": "periodic"}) + resolve(doc) + with pytest.raises(SymbolicError, match="a file domain cannot be periodic"): + adr_problem(doc) + + +def test_a_file_domain_loads_with_bounds(tmp_path): + text = (EXAMPLES / "poisson.ymf").read_text(encoding="utf-8").replace( + ''' domain: "Ω = [0, 1] × [0, 1]"''', ''' domain: + kind: file + format: tetgen + path: cube + bounds: {lower: [0, 0, 0], upper: [1, 1, 1]} + boundary_tags: {walls: 1}''').replace(''' where: "∂Ω"''', ''' tags: [walls]''').replace( + "mesh: {cells: [4, 8, 16]}", "mesh: {levels: [1, 2]}") + (tmp_path / "p.ymf").write_text(text, encoding="utf-8") + doc = load_ymf(tmp_path / "p.ymf") + assert doc["model"]["strong_form"]["domain"]["bounds"]["upper"] == [1.0, 1.0, 1.0] + assert reached(doc)[0] == "discretized" diff --git a/tests/test_gauges.py b/tests/test_gauges.py new file mode 100644 index 0000000..baed51b --- /dev/null +++ b/tests/test_gauges.py @@ -0,0 +1,69 @@ +"""Gauges: what a run measures as it goes. + +The shape FORCE_HISTORY_GAUGES_PLAN step 4 planned (C22): point, line and +line-integral gauges of named fields, and the boundary-integral gauge of +forces and moments per boundary region, which replaces RANS2P's and +RANS3PF's history files. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 + +TWO_PHASE = Path(__file__).parent / "fixtures" / "two_phase" + + +@pytest.fixture +def dambreak(): + return load_ymf(TWO_PHASE / "dambreak.ymf") + + +def _with(doc, gauge): + doc = copy.deepcopy(doc) + doc["model"]["gauges"] = [gauge] + return doc + + +def test_the_dam_break_carries_its_decks_gauges(dambreak): + gauges = {g["name"]: g for g in dambreak["model"]["gauges"]} + assert len(gauges["heights"]["lines"]) == 4 + assert gauges["pressure"]["points"][2] == [3.22, 0.12] + assert gauges["pressure"]["file"] == "pressure.csv" + + +def test_a_boundary_integral_gauge(dambreak): + resolve(_with(dambreak, {"kind": "boundary_integral", "tags": ["walls"], + "quantities": ["force", "moment"], + "moment_point": [1.61, 0.9], "sample_rate": 0.01, + "file": "wall_forces.csv"})) + + +@pytest.mark.parametrize("gauge, message", [ + ({"kind": "point", "fields": ["p"]}, "a point gauge needs points"), + ({"kind": "line", "points": [[0.0, 0.0]], "fields": ["p"], "lines": [[[0, 0], [1, 1]]]}, + "a line gauge takes no points"), + ({"kind": "boundary_integral", "tags": ["walls"]}, "needs quantities"), + ({"kind": "boundary_integral", "tags": ["walls"], "quantities": ["force"], + "fields": ["p"]}, "takes no fields"), + ({"kind": "point", "fields": ["q"], "points": [[0.0, 0.0]]}, + "names 'q', which is not an unknown or coefficient"), + ({"kind": "boundary_integral", "tags": ["lid"], "quantities": ["force"]}, + "names 'lid', which is not a boundary region"), + ({"kind": "line_integral", "fields": ["φ"], "lines": [[[0, 0], [1, 1], [2, 2]]]}, + "a line is its two ends, not 3 points"), + ({"kind": "point", "fields": ["p"], "points": [[0.0, 0.0], [1.0, 1.0, 1.0]]}, + "do not all have the same dimension"), +]) +def test_what_a_gauge_needs_is_checked(dambreak, gauge, message): + with pytest.raises(YmfReferenceError, match=message): + resolve(_with(dambreak, gauge)) + + +def test_a_gauge_of_a_coefficient_of_the_state(dambreak): + resolve(_with(dambreak, {"kind": "point", "fields": ["ρ"], "points": [[1.0, 0.5]]})) diff --git a/tests/test_mesh.py b/tests/test_mesh.py new file mode 100644 index 0000000..e2e6f3b --- /dev/null +++ b/tests/test_mesh.py @@ -0,0 +1,301 @@ +"""A mesh as a ymf grid: tagged facets, provenance, and the order carried by the cell type.""" + +import pytest + +from ymf.archive import ( + KNOWN_TOPOLOGY_TYPES, + QUADRATIC_EDGES, + YmfArchiveError, + add_instant, + attribute, + data_item_for, + facets, + geometry, + geometry_order, + inline_approximation, + new_approximation, + nodes_per_cell, + provenance, + read_ymf, + select_approximation, + topology, + uniform_grid, + validate_approximation, + write_ymf, +) + +np = pytest.importorskip("numpy") + +from ymf.mesh import ( # noqa: E402 + check_facets, + edge_midpoint_offsets, + edge_node_permutation, + file_provenance, + reorder_cells, +) + + +# -------------------------------------------------------------------------- +# fixtures: a unit square in two triangles, its four sides tagged +# -------------------------------------------------------------------------- + +NODES = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]]) +CELLS = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int32) +SIDES = np.array([[0, 1], [1, 2], [2, 3], [3, 0]], dtype=np.int32) +SIDE_TAGS = np.array([1, 2, 0, 4], dtype=np.int32) # top side untagged + + +def square_space(tags=SIDE_TAGS, view="global", source=None, sides=SIDES): + return uniform_grid( + topology("Triangle", len(CELLS), data_item_for(CELLS, inline=True)), + geometry(data_item_for(NODES, inline=True)), + [attribute("elementMaterialTypes", + data_item_for(np.ones(2, dtype=np.int32), inline=True), center="Cell")], + facets=facets(topology("Polyline", len(sides), data_item_for(sides, inline=True), + nodes_per_element=2), + data_item_for(np.asarray(tags, dtype=np.int32), inline=True), view=view), + provenance=source) + + +def mesh_domain(space): + domain = new_approximation() + add_instant(domain, 0.0, {"Mesh Spatial_Domain": space}) + return domain + + +SOURCE = provenance("a mesh converter", version="0.1", + files=[{"name": "square.msh", "sha256": "ab" * 32}], + extra={"units": "m"}) + + +# -------------------------------------------------------------------------- +# facets and provenance survive the archive's round trips +# -------------------------------------------------------------------------- + + +def test_a_mesh_with_facets_and_provenance_round_trips(tmp_path): + domain = validate_approximation(mesh_domain(square_space(source=SOURCE))) + write_ymf(domain, tmp_path / "square.ymf") + back, _ = read_ymf(tmp_path / "square.ymf") + assert validate_approximation(back) == domain + space = back["times"][0]["spaces"]["Mesh Spatial_Domain"] + assert space["provenance"] == {"tool": "a mesh converter", "version": "0.1", + "files": [{"name": "square.msh", "sha256": "ab" * 32}], + "extra": {"units": "m"}} + assert space["facets"]["view"] == "global" + assert check_facets(space) == {0: 1, 1: 1, 2: 1, 4: 1} + + +def test_facets_held_in_hdf5_are_read_and_checked(tmp_path): + h5py = pytest.importorskip("h5py") + with h5py.File(tmp_path / "square.h5", "w") as f: + for name, array in (("nodes", NODES), ("elements", CELLS), + ("facets", SIDES), ("facet_tags", SIDE_TAGS)): + f[name] = array + ref = lambda name, a: data_item_for(a, "square.h5:/" + name) # noqa: E731 + space = uniform_grid( + topology("Triangle", 2, ref("elements", CELLS)), geometry(ref("nodes", NODES)), + facets=facets(topology("Polyline", 4, ref("facets", SIDES), nodes_per_element=2), + ref("facet_tags", SIDE_TAGS))) + domain = validate_approximation(mesh_domain(space)) + write_ymf(domain, tmp_path / "square.ymf") + back, _ = read_ymf(tmp_path / "square.ymf") + space = validate_approximation(back)["times"][0]["spaces"]["Mesh Spatial_Domain"] + assert check_facets(space, tmp_path) == {0: 1, 1: 1, 2: 1, 4: 1} + # inlining carries the facet arrays along with the rest + inlined = inline_approximation(back, tmp_path) + facet_list = inlined["times"][0]["spaces"]["Mesh Spatial_Domain"]["facets"] + assert facet_list["tags"]["data_item"]["values"] == SIDE_TAGS.tolist() + assert facet_list["topology"]["data_item"]["values"] == SIDES.reshape(-1).tolist() + + +def test_selecting_fields_keeps_the_facets(): + domain = mesh_domain(square_space()) + picked = select_approximation(domain, fields=["elementMaterialTypes"]) + assert "facets" in picked["times"][0]["spaces"]["Mesh Spatial_Domain"] + + +def test_ymf2xmf_writes_a_mesh_with_facets(tmp_path): + from ymf.xdmf import write_xdmf + write_xdmf(mesh_domain(square_space(source=SOURCE)), tmp_path / "square.xmf") + assert "Triangle" in (tmp_path / "square.xmf").read_text() + + +def test_file_provenance_records_names_and_checksums(tmp_path): + (tmp_path / "a.node").write_bytes(b"4 2 0 0\n") + block = file_provenance("a generator", [tmp_path / "a.node"], version="1.6") + import hashlib + assert block == {"tool": "a generator", "version": "1.6", + "files": [{"name": "a.node", + "sha256": hashlib.sha256(b"4 2 0 0\n").hexdigest()}]} + validate_approximation(mesh_domain(square_space(source=block))) + + +# -------------------------------------------------------------------------- +# what the facet list refuses +# -------------------------------------------------------------------------- + + +def test_a_global_mesh_never_stores_a_rank_interface(): + with pytest.raises(YmfArchiveError, match="-1 marks a rank interface"): + validate_approximation(mesh_domain(square_space(tags=[1, -1, 0, 4]))) + + +def test_a_subdomain_view_may_mark_rank_interfaces(): + space = validate_approximation(mesh_domain( + square_space(tags=[1, -1, 0, 4], view="subdomain")))["times"][0]["spaces"] + assert check_facets(space["Mesh Spatial_Domain"]) == {-1: 1, 0: 1, 1: 1, 4: 1} + + +@pytest.mark.parametrize("tags, match", [ + ([1, -2, 0, 4], "tag -2 is not a tag"), + ([1, 2, 0], r"one tag per facet, \[4\]"), +]) +def test_bad_tags_are_refused(tags, match): + with pytest.raises(YmfArchiveError, match=match): + validate_approximation(mesh_domain(square_space(tags=tags))) + + +def test_tags_are_integers(): + space = square_space() + space["facets"]["tags"]["data_item"] = data_item_for(np.array([1.0, 2.0, 0.0, 4.0]), + inline=True) + with pytest.raises(YmfArchiveError, match="tags are integers"): + validate_approximation(mesh_domain(space)) + + +def test_the_facet_type_follows_the_cell_type(): + space = square_space() + space["facets"]["topology"]["type"] = "Triangle" + with pytest.raises(YmfArchiveError, match="facets of Triangle cells are Polyline"): + validate_approximation(mesh_domain(space)) + + +def test_an_unknown_view_is_refused(): + with pytest.raises(YmfArchiveError, match="'partitioned' is not one of"): + validate_approximation(mesh_domain(square_space(view="partitioned"))) + + +def test_a_facet_must_be_a_face_of_a_cell(): + # (1, 3) joins two corners of the square but is an edge of neither triangle + sides = np.array([[0, 1], [1, 3]], dtype=np.int32) + space = validate_approximation(mesh_domain(square_space(sides=sides, tags=[1, 2]))) + with pytest.raises(YmfArchiveError, match=r"facet 1 \(nodes \[1, 3\]\) is not a face"): + check_facets(space["times"][0]["spaces"]["Mesh Spatial_Domain"]) + + +def test_a_facet_node_must_be_a_node_of_the_grid(): + space = validate_approximation(mesh_domain(square_space( + sides=np.array([[0, 1], [3, 4]], dtype=np.int32), tags=[1, 2]))) + with pytest.raises(YmfArchiveError, match=r"must lie in \[0, 4\)"): + check_facets(space["times"][0]["spaces"]["Mesh Spatial_Domain"]) + + +@pytest.mark.parametrize("block, match", [ + ({"version": "1"}, "provenance has no 'tool'"), + ({"tool": "t", "codes": ["x"]}, "put anything else under extra"), + ({"tool": "t", "files": [{"name": "a"}]}, r"a file is \{name, sha256\}"), + ({"tool": "t", "files": [{"name": "a", "sha256": "1234"}]}, "not a SHA-256 digest"), +]) +def test_bad_provenance_is_refused(block, match): + with pytest.raises(YmfArchiveError, match=match): + validate_approximation(mesh_domain(square_space(source=block))) + + +# -------------------------------------------------------------------------- +# geometry order is the cell type's +# -------------------------------------------------------------------------- + + +def test_every_known_type_has_an_order(): + for t in KNOWN_TOPOLOGY_TYPES: + order = geometry_order(t) + assert order in (1, 2) or (t == "Mixed" and order is None), t + assert geometry_order("Tri_6") == geometry_order("Tet_10") == 2 + assert geometry_order("Tetrahedron") == 1 + with pytest.raises(YmfArchiveError, match="not a known topology type"): + geometry_order("Tet_11") + + +@pytest.mark.parametrize("cell_type, count", [ + ("Triangle", 3), ("Tri_6", 6), ("Quad_8", 8), ("Tetrahedron", 4), ("Tet_10", 10), + ("Tetrahedron_10", 10), ("Pyramid_13", 13), ("Wedge_15", 15), ("Hex_20", 20), + ("Polyline", None), ("Mixed", None), +]) +def test_nodes_per_cell(cell_type, count): + assert nodes_per_cell(cell_type) == count + + +def test_quadratic_edges_name_each_edge_once(): + for t, edges in QUADRATIC_EDGES.items(): + assert len({frozenset(e) for e in edges}) == len(edges), t + + +# -------------------------------------------------------------------------- +# the Tet_10 order: a source that numbers edge nodes differently is permuted +# -------------------------------------------------------------------------- + +#: Two edge-node conventions met in practice, as test data. The first is a +#: mesh generator's: the last two edge nodes the other way round from ymf's. +#: The second is how Proteus numbers its quadratic tetrahedron's nodes +#: (FemTools: increments of 1, then 2, then 3 between vertex numbers). +GENERATOR_TET10 = ((0, 1), (1, 2), (2, 0), (0, 3), (2, 3), (1, 3)) +INCREMENT_TET10 = ((0, 1), (1, 2), (2, 3), (0, 2), (1, 3), (0, 3)) + + +def _straight_tet10(edges): + """One straight-sided Tet_10 with its edge nodes at the midpoints of ``edges``.""" + vertices = np.array([[0.0, 0, 0], [1.0, 0, 0], [0.0, 1, 0], [0.0, 0, 1]]) + nodes = np.vstack([vertices] + [0.5 * (vertices[a] + vertices[b]) for a, b in edges]) + return nodes, np.arange(10)[None, :] + + +def test_ymfs_own_order_is_the_identity(): + assert edge_node_permutation("Tet_10", QUADRATIC_EDGES["Tet_10"]) == list(range(10)) + nodes, cells = _straight_tet10(QUADRATIC_EDGES["Tet_10"]) + assert edge_midpoint_offsets("Tet_10", nodes, cells).max() < 1e-15 + + +@pytest.mark.parametrize("source, perm", [ + (GENERATOR_TET10, [0, 1, 2, 3, 4, 5, 6, 7, 9, 8]), + (INCREMENT_TET10, [0, 1, 2, 3, 4, 5, 7, 9, 8, 6]), +]) +def test_a_source_tet10_is_put_in_ymf_order(source, perm): + assert edge_node_permutation("Tet_10", source) == perm + nodes, cells = _straight_tet10(source) + # read as ymf's order the edge nodes are off their edges... + assert edge_midpoint_offsets("Tet_10", nodes, cells).max() > 0.3 + # ...and after the permutation every one is at its edge's midpoint + fixed = reorder_cells(cells, "Tet_10", source) + assert edge_midpoint_offsets("Tet_10", nodes, fixed).max() < 1e-15 + + +def test_a_pair_written_either_way_round_is_the_same_edge(): + flipped = tuple((b, a) for a, b in QUADRATIC_EDGES["Tri_6"]) + assert edge_node_permutation("Tri_6", flipped) == list(range(6)) + + +@pytest.mark.parametrize("cell_type, source, match", [ + ("Tetrahedron", QUADRATIC_EDGES["Tet_10"], "not a quadratic cell type"), + ("Tet_10", QUADRATIC_EDGES["Tet_10"][:5], "has 6 distinct edges"), + ("Tri_6", ((0, 1), (1, 2), (1, 3)), r"edge \(2, 0\) is not among"), +]) +def test_a_wrong_source_convention_is_refused(cell_type, source, match): + with pytest.raises(YmfArchiveError, match=match): + edge_node_permutation(cell_type, source) + + +def test_a_tet10_mesh_carries_tri6_facets(tmp_path): + nodes, cells = _straight_tet10(QUADRATIC_EDGES["Tet_10"]) + # the face (0, 1, 2) of the tetrahedron, with its edge nodes 4, 5, 6 + face = np.array([[0, 1, 2, 4, 5, 6]], dtype=np.int32) + space = uniform_grid( + topology("Tet_10", 1, data_item_for(cells.astype(np.int32), inline=True)), + geometry(data_item_for(nodes, inline=True)), + facets=facets(topology("Tri_6", 1, data_item_for(face, inline=True)), + data_item_for(np.array([3], dtype=np.int32), inline=True))) + space = validate_approximation(mesh_domain(space))["times"][0]["spaces"]["Mesh Spatial_Domain"] + assert check_facets(space) == {3: 1} + space["facets"]["topology"]["type"] = "Triangle" + with pytest.raises(YmfArchiveError, match="facets of Tet_10 cells are Tri_6 or Triangle_6"): + validate_approximation(mesh_domain(space)) diff --git a/tests/test_named_equations.py b/tests/test_named_equations.py new file mode 100644 index 0000000..1802a67 --- /dev/null +++ b/tests/test_named_equations.py @@ -0,0 +1,153 @@ +"""Named equations: ``{name, for, equation}`` beside the bare strings. + +A name is what a discretization's sub-problem refers to and what an +extending file merges onto; ``for`` says which unknown an equation +determines, so the order they are written in stops mattering. +""" + +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.normalize import equations_in_order, normalize_equations # noqa: E402 +from ymf.references import YmfReferenceError # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" +NS = EXAMPLES / "navier_stokes" + +MOMENTUM = "ρ ∂v/∂t + ∇·(ρ v⊗v) − μΔv + ∇p = f in Ω" +CONTINUITY = "∇·v = 0 in Ω" + + +def _named(**extra): + return [{"name": "continuity", "for": "p", "equation": CONTINUITY}, + {"name": "momentum", "for": "v", "equation": MOMENTUM, **extra}] + + +# --- the helpers, no schema ------------------------------------------------- + + +def test_bare_strings_keep_their_order(): + strong = {"unknowns": ["v", "p"], "equations": [MOMENTUM, CONTINUITY]} + assert equations_in_order(strong) == [MOMENTUM, CONTINUITY] + + +def test_for_puts_them_in_the_unknowns_order(): + strong = {"unknowns": ["v", "p"], "equations": _named()} + assert equations_in_order(strong) == [MOMENTUM, CONTINUITY] + + +def test_a_named_equation_without_for_keeps_its_place(): + strong = {"unknowns": ["u"], "equations": [{"name": "heat", "equation": "-Δu = f"}]} + assert normalize_equations(strong["equations"]) == [ + {"name": "heat", "for": None, "equation": "-Δu = f"}] + assert equations_in_order(strong) == ["-Δu = f"] + + +@pytest.mark.parametrize("equations, message", [ + ([{"name": "momentum", "for": "v", "equation": MOMENTUM}, CONTINUITY], + "give 'for' on every equation or on none"), + ([{"name": "momentum", "for": "w", "equation": MOMENTUM}, + {"name": "continuity", "for": "p", "equation": CONTINUITY}], + "is for 'w', which is not an unknown"), + ([{"name": "a", "for": "v", "equation": MOMENTUM}, + {"name": "b", "for": "v", "equation": CONTINUITY}], + "a and b are both for v"), + ([{"name": "momentum", "for": "v", "equation": MOMENTUM}], + "none is for p"), +]) +def test_what_does_not_say_which_equation_is_for_which_unknown(equations, message): + with pytest.raises(ValueError, match=message): + equations_in_order({"unknowns": ["v", "p"], "equations": equations}) + + +def test_an_incomplete_set_is_fine_where_the_chain_completes_it(): + strong = {"unknowns": ["v", "p"], + "equations": [{"name": "momentum", "for": "v", "equation": MOMENTUM}]} + assert equations_in_order(strong, complete=False) == [MOMENTUM] + + +# --- through the loader ----------------------------------------------------- + + +def _poiseuille_with(tmp_path, equations_yaml, model_extra=""): + """plane_poiseuille, extending a copy of the model whose equations are given.""" + model = (NS / "navier_stokes_model.ymf").read_text(encoding="utf-8") + start = model.index(" equations:") + end = model.index(" # The coefficients the equations use") + (tmp_path / "navier_stokes_model.ymf").write_text( + model[:start] + equations_yaml + model[end:] + model_extra, encoding="utf-8") + (tmp_path / "plane_poiseuille.ymf").write_text( + (NS / "plane_poiseuille.ymf").read_text(encoding="utf-8"), encoding="utf-8") + return tmp_path / "plane_poiseuille.ymf" + + +NAMED_YAML = ( + " equations:\n" + " - name: continuity\n" + " for: p\n" + ' equation: "∇·v = 0 in Ω"\n' + " - name: momentum\n" + " for: v\n" + ' equation: "ρ ∂v/∂t + ∇·(ρ v⊗v) − μΔv + ∇p = f in Ω"\n') + + +def test_named_and_reordered_equations_give_the_same_problem(tmp_path): + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + named = load_ymf(_poiseuille_with(tmp_path, NAMED_YAML)) + plain = load_ymf(NS / "plane_poiseuille.ymf") + assert adr_problem(named)["adr"] == adr_problem(plain)["adr"] + + +def test_an_extending_file_replaces_an_equation_by_its_name(tmp_path): + path = _poiseuille_with(tmp_path, NAMED_YAML) + (tmp_path / "variant.ymf").write_text( + "extends: plane_poiseuille.ymf\n" + "model:\n" + " strong_form:\n" + " equations:\n" + " - name: momentum\n" + ' equation: "ρ ∂v/∂t + ∇·(ρ v⊗v) − ∇·(μ∇v) + ∇p = f in Ω"\n', + encoding="utf-8") + assert path.exists() + doc = load_ymf(tmp_path / "variant.ymf") + equations = doc["model"]["strong_form"]["equations"] + assert [e["name"] for e in equations] == ["continuity", "momentum"] + assert equations[1]["for"] == "v" # kept from the model + assert "∇·(μ∇v)" in equations[1]["equation"] + assert any(o["path"] == "model.strong_form.equations[momentum].equation" + for o in doc["composition"]["overrides"]) + + +def test_an_abstract_model_may_leave_an_unknown_to_the_files_that_extend_it(tmp_path): + one = (" equations:\n" + " - name: momentum\n" + " for: v\n" + ' equation: "ρ ∂v/∂t + ∇·(ρ v⊗v) − μΔv + ∇p = f in Ω"\n') + path = _poiseuille_with(tmp_path, one) + load_ymf(tmp_path / "navier_stokes_model.ymf") # abstract: fine + with pytest.raises(YmfReferenceError, match="none is for p"): + load_ymf(path) # well-posed: not + + +@pytest.mark.parametrize("equations_yaml, message", [ + (NAMED_YAML.replace("for: p", "for: q"), "is for 'q', which is not an unknown"), + (NAMED_YAML.replace("name: continuity", "name: momentum"), "momentum named more than once"), +]) +def test_bad_names_are_refused_on_load(tmp_path, equations_yaml, message): + with pytest.raises(YmfReferenceError, match=message): + load_ymf(_poiseuille_with(tmp_path, equations_yaml)) + + +def test_the_flux_units_check_reads_named_equations(tmp_path): + pytest.importorskip("sympy") + pytest.importorskip("pint") + from ymf import check_units + named = load_ymf(_poiseuille_with(tmp_path, NAMED_YAML)) + plain = load_ymf(NS / "plane_poiseuille.ymf") + assert ([str(i) for i in check_units(named, missing_units="warn").issues] + == [str(i) for i in check_units(plain, missing_units="warn").issues]) diff --git a/tests/test_notation_functions.py b/tests/test_notation_functions.py new file mode 100644 index 0000000..1169f26 --- /dev/null +++ b/tests/test_notation_functions.py @@ -0,0 +1,101 @@ +"""heaviside, sign, min, max, norm and sym in the notation. + +The two-phase level-set model writes χ = heaviside(φ); its initial signed +distance is a min/max/sqrt formula; its momentum equation uses the strain +rate sym(∇v); shallow-water friction uses norm(q). heaviside and sign are +sharp: their derivatives are Dirac deltas, which a solver cannot evaluate, +so a sharp function of an unknown is refused where a derivative is needed +until a discretization regularizes it. +""" + +import math + +import pytest + +sympy = pytest.importorskip("sympy") + +from ymf.symbolic import Space, SymbolicError, adr_form, parse_equation # noqa: E402 +from ymf.symbolic import parse_expression, to_code, to_notation # noqa: E402 + + +def _space(**unknowns): + return Space(2, unknowns or {"u": 0}) + + +def _evaluate(code, **names): + numpy = pytest.importorskip("numpy") # not in ymf[symbolic] alone + return eval(code, {"__builtins__": {}, "numpy": numpy}, names) # noqa: S307 + + +@pytest.mark.parametrize("text", [ + "heaviside(u)", "sign(u - 1)", "min(x, y)", "max(x, y, 0)", + "min(max(x - 1.2, y - 0.6), 0)", +]) +def test_the_new_functions_round_trip(text): + space = _space() + value = parse_expression(text, space) + assert parse_expression(to_notation(value), space) == value + + +def test_min_and_max_are_numpy_alone(): + code = to_code(parse_expression("max(x, y, 0) + min(x, 2)", _space())) + assert "functools" not in code + assert _evaluate(code, x=-1.0, y=-2.0) == 0.0 + -1.0 + + +def test_the_box_signed_distance_of_the_dam_break_column(): + """The note's φ0 equals damBreakLevelSet.py's phi_init in all three branches.""" + text = "sqrt(max(x − 1.2, 0)² + max(y − 0.6, 0)²) + min(max(x − 1.2, y − 0.6), 0)" + code = to_code(parse_expression(text, _space())) + + def phi_init(x, y): # test/TwoPhaseFlow/damBreakLevelSet.py + dx, dy = x - 1.2, y - 0.6 + if dx < 0 and dy < 0: + return max(dx, dy) + if dx > 0 and dy > 0: + return math.hypot(dx, dy) + return max(dx, dy) + + for x, y in [(0.3, 0.2), (1.1, 0.55), (2.0, 1.4), (1.5, 0.1), (0.4, 1.0), (3.2, 1.8)]: + assert math.isclose(_evaluate(code, x=x, y=y), phi_init(x, y), abs_tol=1e-15) + + +def test_norm_of_a_vector_and_of_a_scalar(): + space = _space(q=1) + q0, q1 = space.component_symbols + assert parse_expression("norm(q)", space) == sympy.sqrt(q0 ** 2 + q1 ** 2) + assert parse_expression("norm(x - 1)", space) == sympy.Abs(space.x[0] - 1) + + +def test_sym_is_the_symmetric_part(): + space = _space(v=1) + strain = parse_expression("sym(∇v)", space) + assert strain == strain.T + with pytest.raises(SymbolicError, match="square tensor"): + parse_expression("sym(v)", space) + + +def test_the_symmetric_stress_is_a_diffusion_with_cross_terms(): + """∇·(2μ sym(∇v)) = ∇·(μ(∇v + ∇vᵀ)): v_0's equation diffuses v_1 too.""" + space = Space(2, {"v": 1, "p": 0}, {"μ": sympy.Float(0.5)}) + residuals = [] + for text in ["−∇·(2μ sym(∇v)) + ∇p = 0", "∇·v = 0"]: + residuals.extend(parse_equation(text, space)) + momentum = adr_form(residuals, space)["equations"][0] + assert set(momentum["diffusion"]) == {"v_0", "v_1"} + + +def test_a_sharp_function_of_an_unknown_has_no_derivative_to_give(): + space = _space() + residuals = parse_equation("∂(heaviside(u))/∂t = 0", space) + with pytest.raises(SymbolicError, match="regularizes it"): + adr_form(residuals, space) + + +def test_a_sharp_function_of_the_coordinates_is_data(): + """A zone indicator in x, y: a reaction coefficient, never differentiated.""" + space = _space() + residuals = parse_equation("∂u/∂t = ∇·∇u − heaviside(x − 0.5) u", space) + reaction = adr_form(residuals, space)["equations"][0]["reaction"] + assert reaction["depends_on"] == {"u": "linear"} + assert "numpy.select" in reaction["r"] diff --git a/tests/test_periodic_regions.py b/tests/test_periodic_regions.py new file mode 100644 index 0000000..9e07e6b --- /dev/null +++ b/tests/test_periodic_regions.py @@ -0,0 +1,91 @@ +"""Periodic conditions on a piecewise_linear region (C27 Phase 5, D3). + +Proteus pairs a node on a periodic face with the node the translation puts +on the opposite face, so a mesh has partners only if the boundary pieces on +the two faces match. ymf checks that at posing: the segments (2-D) or +facets (3-D) on the lower face, moved along the axis, must be exactly those +on the upper face. The mesher adding no points to them, and the partners +themselves, are checked by the solver after meshing. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") +sympy = pytest.importorskip("sympy") + +from ymf.schema import load_ymf # noqa: E402 +from ymf.symbolic import SymbolicError, adr_problem # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + +#: a channel with a notch in its top, away from the periodic x faces +NOTCHED = { + "kind": "piecewise_linear", "dimension": 2, + "vertices": [[0.0, 0.0], [4.0, 0.0], [4.0, 1.0], [3.0, 1.0], + [3.0, 0.5], [1.0, 0.5], [1.0, 1.0], [0.0, 1.0]], + "segments": [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 0]], +} + +#: the unit cube as a PLC, one quadrilateral facet per face +CUBE = { + "kind": "piecewise_linear", "dimension": 3, + "vertices": [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], + [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]], + "facets": [[[0, 1, 2, 3]], [[4, 5, 6, 7]], [[0, 1, 5, 4]], + [[2, 3, 7, 6]], [[1, 2, 6, 5]], [[0, 3, 7, 4]]], +} + + +def _periodic(domain, ends="x = 0 or x = 4"): + doc = load_ymf(EXAMPLES / "poisson.ymf") + strong = doc["model"]["strong_form"] + strong["domain"] = copy.deepcopy(domain) + strong["boundary_regions"] = [{"name": "ends", "where": ends}, + {"name": "rest", "where": "∂Ω"}] + strong["boundary_conditions"] = [ + {"region": "ends", "variable": "u", "type": "periodic"}, + {"region": "rest", "variable": "u", "type": "dirichlet", "value": 0.0}] + doc["solution_paths"]["analytical"] = [] + doc.pop("vvuq", None) + return doc + + +def test_a_notched_channel_is_periodic_along_x(): + problem = adr_problem(_periodic(NOTCHED)) + assert problem["periodic"] == [{"component": "u", "region": "ends", "axes": [0]}] + assert problem["domain"]["region"]["segments"] == NOTCHED["segments"] + + +def test_a_face_cut_differently_on_the_two_sides_is_refused(): + split = copy.deepcopy(NOTCHED) + split["vertices"].append([4.0, 0.5]) # the x = 4 face in two + split["segments"][1] = [1, 8] + split["segments"].append([8, 2]) + with pytest.raises(SymbolicError, match="faces of axis x do not match: 3 of 3"): + adr_problem(_periodic(split)) + + +def test_a_cube_is_periodic_along_x_and_y(): + doc = _periodic(CUBE, ends="x = 0 or x = 1 or y = 0 or y = 1") + assert adr_problem(doc)["periodic"][0]["axes"] == [0, 1] + + +def test_a_cube_face_with_its_facet_rotated_still_matches(): + """Facets match as point sets, not by vertex order.""" + cube = copy.deepcopy(CUBE) + cube["facets"][4] = [[2, 6, 5, 1]] + assert adr_problem(_periodic(cube, ends="x = 0 or x = 1"))["periodic"][0]["axes"] == [0] + + +def test_split_problem_checks_the_faces_too(): + from ymf.symbolic import split_problem + doc = _periodic(NOTCHED) + doc["model"]["strong_form"]["equations"] = [ + {"name": "poisson", "for": "u", "equation": "-Δu = f in Ω"}] + entry = doc["solution_paths"]["discretizations"][0] + entry["subproblems"] = [{"name": "one", "solves": ["poisson"]}] + (sub,) = split_problem(doc, entry["name"])["subproblems"] + assert [c["axes"] for c in sub["boundary_conditions"] if c["type"] == "periodic"] == [[0]] diff --git a/tests/test_printer.py b/tests/test_printer.py index 069cca7..569231c 100644 --- a/tests/test_printer.py +++ b/tests/test_printer.py @@ -64,7 +64,8 @@ def _space_of(doc): ranks = {(u if isinstance(u, str) else u["name"]): (0 if isinstance(u, str) else int(u.get("rank", 0))) for u in strong["unknowns"]} - return Space(dim, ranks, _coefficients(strong.get("coefficients") or {}, dim)) + return Space(dim, ranks, _coefficients(strong.get("coefficients") or {}, dim, + unknowns=ranks)) def _loadable(path): @@ -91,7 +92,8 @@ def test_every_example_spec_round_trips(path): assert strong.get("equations"), "an example spec with no equations" space = _space_of(doc) count = 0 - for text in strong["equations"]: + from ymf.normalize import equations_in_order + for text in equations_in_order(strong): for residual in parse_equation(text, space): round_trip(residual, space) count += 1 diff --git a/tests/test_references.py b/tests/test_references.py index b3cd657..fd0eadc 100644 --- a/tests/test_references.py +++ b/tests/test_references.py @@ -217,7 +217,7 @@ def test_the_old_branch_ref_spelling_says_what_it_became(tmp_path): def test_a_reference_to_a_missing_analytical_branch_is_caught(tmp_path): text = spec(weak_form("only"), discretization("d"), analytical("exact"), VVUQ % ' reference:\n analytical: "nope"\n') - with pytest.raises(YmfReferenceError, match="not a analytical solution"): + with pytest.raises(YmfReferenceError, match="not an analytical solution"): load_ymf(write(tmp_path, text)) diff --git a/tests/test_shallow_water.py b/tests/test_shallow_water.py new file mode 100644 index 0000000..e7031dc --- /dev/null +++ b/tests/test_shallow_water.py @@ -0,0 +1,120 @@ +"""The shallow-water family (C42's second, C35's model) as ymf documents. + +tests/fixtures/shallow_water: the model (depth h, discharge q; porosity, +Manning friction, Darcy-Forchheimer drag as coefficients of the state), the +numerics Proteus's SW2DCV uses (one sub-problem), a lake at rest over a bump +(well-balancing, exact), and the Torino flume (a dam break down a 3 % slope +into a rod canopy, from torino_flume/dswe/torino.py). +""" + +import math +from pathlib import Path + +import pytest + +sympy = pytest.importorskip("sympy") +pytest.importorskip("strictyaml") + +from ymf.schema import load_ymf # noqa: E402 +from ymf.state import reached # noqa: E402 +from ymf.symbolic import D, Dt, split_problem, strong_form_equations # noqa: E402 + +SW = Path(__file__).parent / "fixtures" / "shallow_water" + + +@pytest.fixture(scope="module") +def lake(): + return load_ymf(SW / "lake_at_rest.ymf") + + +@pytest.fixture(scope="module") +def torino(): + return load_ymf(SW / "torino.ymf") + + +def _swe(doc): + (sub,) = split_problem(doc, "convex_limited")["subproblems"] + return sub + + +def _equation(sub, component): + return next(e for e in sub["adr"]["equations"] if e["component"] == component) + + +def _evaluate(code, **names): + numpy = pytest.importorskip("numpy") # not in ymf[symbolic] alone + return eval(code, {"__builtins__": {}, "numpy": numpy}, names) # noqa: S307 + + +def test_both_cases_are_discretized(lake, torino): + assert reached(lake)[0] == "discretized" and reached(torino)[0] == "discretized" + # a case changes a free-form scheme option of the shared numerics + assert torino["solution_paths"]["discretizations"][0]["discretization"][ + "entropy_viscosity_cE"] == "5.0" + + +def test_shallow_water_is_one_adr_system(lake): + swe = _swe(lake) + assert swe["unknowns"] == {"h": ["h"], "q": ["q_0", "q_1"]} + momentum = _equation(swe, "q_0") + assert momentum["advection"]["f"][0] == "4.905*h**2 + 1.0*q_0**2/h" + assert "q_0*numpy.sqrt(q_0**2 + q_1**2)/h**(7/3)" in momentum["reaction"]["r"] + assert swe["coefficients"]["τb"]["depends_on"] == { + "h": "nonlinear", "q_0": "nonlinear", "q_1": "nonlinear"} + assert [(c["type"], c["variable"]) for c in swe["boundary_conditions"]] == [("slip", "q")] + + +def test_the_lake_at_rest_solves_the_model_exactly(lake): + """Well-balancing of the strong form: h = 0.5 - Z, q = 0 makes every + residual vanish identically, the bed's kink included.""" + space, residuals = strong_form_equations(lake) + bed = space.coefficients["Z"] # the spec's own bed + state = {sympy.Symbol("h", real=True): sympy.Float(0.5) - bed, + sympy.Symbol("q_0", real=True): 0, sympy.Symbol("q_1", real=True): 0} + + def substitute(r): + r = r.replace(D, lambda e, axis: sympy.diff(sympy.sympify(e).subs(state), + space.x[int(axis)])) + r = r.replace(Dt, lambda e: sympy.diff(sympy.sympify(e).subs(state), space.t)) + return sympy.simplify(r.subs(state).doit()) + + assert [substitute(r) for r in residuals] == [0, 0, 0] + + +def test_the_bed_kink_is_not_regularized(lake): + """Only sharp functions as written are smoothed; the bump's max gives a + bounded step in ∇Z, which is left as it is.""" + swe = _swe(lake) + assert "regularized" not in swe + assert "[0,1/2,1]" in _equation(swe, "q_0")["reaction"]["r"] + + +def test_the_torino_canopy(torino): + swe = _swe(torino) + assert swe["regularized"] == {"heaviside": {"by": "smoothed_heaviside", "width": "he"}} + porosity = 1 - 194.0 * math.pi * 0.006 ** 2 / 4 + mass = _equation(swe, "h")["mass"]["m"] + assert math.isclose(_evaluate(mass, h=1.0, x=7.0, he=0.0725), porosity, rel_tol=1e-12) + assert _evaluate(mass, h=1.0, x=2.0, he=0.0725) == 1.0 + assert [(c["type"], c["region"]) for c in swe["boundary_conditions"]] == [ + ("slip", "walls"), ("open", "open_end")] + + +def test_the_torino_drag_is_the_decks_froude_law(torino): + """Fd = (α + β|q|/h) q with β = C_D min(1, 0.9 Fr^-0.8) m D / 2 in the + canopy, Fr = |q|/(h sqrt(g h)): SW2DCV's -(α + β|u|) h u.""" + swe = _swe(torino) + h, q0, x = 0.1, 0.05, 7.0 + fr = (q0 / h) / math.sqrt(9.81 * h) + beta = 1.0 * min(1.0, 0.9 * fr ** -0.8) * 194.0 * 0.006 / 2 + expected = beta * abs(q0) / h * q0 + fd = swe["coefficients"]["Fd"]["code"][0] + assert math.isclose(_evaluate(fd, h=h, q_0=q0, q_1=0.0, x=x, he=0.0725), expected, + rel_tol=1e-12) + assert _evaluate(fd, h=h, q_0=q0, q_1=0.0, x=2.0, he=0.0725) == 0.0 # outside + + +def test_both_cases_pass_the_unit_check(lake, torino): + pytest.importorskip("pint") + from ymf import check_units + assert check_units(lake) and check_units(torino) diff --git a/tests/test_split.py b/tests/test_split.py new file mode 100644 index 0000000..26b705f --- /dev/null +++ b/tests/test_split.py @@ -0,0 +1,263 @@ +"""ymf.symbolic.split: each sub-problem of a split discretization, classified. + +The two-phase dam break (tests/fixtures/two_phase) is not one ADR system, +but each of its five sub-problems is once what it reads from the others is +data: flow is variable-density Navier-Stokes with ρ and μ read from VOF's +carrier, VOF and the level set are transport in a lagged velocity, +redistancing is an eikonal Hamiltonian, mass correction a regularized +reaction-diffusion. +""" + +import copy +import math +from pathlib import Path + +import pytest + +sympy = pytest.importorskip("sympy") +pytest.importorskip("strictyaml") + +from ymf.schema import load_ymf # noqa: E402 +from ymf.symbolic import Space, SymbolicError, adr_form, parse_equation # noqa: E402 +from ymf.symbolic.split import smoothed_heaviside, split_problem # noqa: E402 + +TWO_PHASE = Path(__file__).parent / "fixtures" / "two_phase" + + +@pytest.fixture(scope="module") +def dambreak(): + return load_ymf(TWO_PHASE / "dambreak.ymf") + + +@pytest.fixture(scope="module") +def split(dambreak): + return split_problem(dambreak, "level_set_five") + + +def _sub(split, name): + return next(s for s in split["subproblems"] if s["name"] == name) + + +def _equation(sub, component): + return next(e for e in sub["adr"]["equations"] if e["component"] == component) + + +def _evaluate(code, **names): + numpy = pytest.importorskip("numpy") # not in ymf[symbolic] alone + return eval(code, {"__builtins__": {}, "numpy": numpy}, names) # noqa: S307 + + +# --- the setting --------------------------------------------------------------- + + +def test_the_setting(split): + assert split["domain"] == {"lower": [0.0, 0.0], "upper": [3.22, 1.8]} + assert split["time"] == [0.0, 3.0] + assert split["bound"] == {"h": "element_size", "ε": "1.5*h"} + assert set(split["initial"]) == {"v_0", "v_1", "p", "φ"} + + +def test_each_sub_problem_has_its_own_unknowns(split): + assert [(s["name"], s["unknowns"]) for s in split["subproblems"]] == [ + ("flow", {"v": ["v_0", "v_1"], "p": ["p"]}), + ("vof", {"χh": ["χh"]}), + ("ncls", {"φ": ["φ"]}), + ("rdls", {"φ": ["φ"]}), + ("mcorr", {"φc": ["φc"]}), + ] + + +# --- the five sub-problems ------------------------------------------------------- + + +def test_flow_reads_density_and_viscosity_through_the_vof_carrier(split): + flow = _sub(split, "flow") + assert flow["reads"]["χ"] == {"at": "previous_step", "via": "vof", "carried_by": "χh", + "components": ["χh"]} + assert flow["coefficients"]["ρ"]["code"] == "998.2 - 996.995*χh" + assert flow["coefficients"]["ρ"]["reads"] == ["χh"] + momentum = _equation(flow, "v_0") + assert momentum["mass"]["m"] == "v_0*(998.2 - 996.995*χh)" + assert set(momentum["diffusion"]) == {"v_0", "v_1"} # 2μ sym(∇v) + assert flow["solver"]["max_iterations"] == 100 # refined + assert flow["proteus"] == {"model": "RANS2P", "useVF": "1.0"} # a hint, carried as text + assert "regularized" not in flow # χ is read through χh, not evaluated + + +def test_the_level_set_is_a_hamiltonian_in_the_lagged_velocity(split): + ncls = _sub(split, "ncls") + assert ncls["reads"] == {"v": {"at": "current_step", "via": "flow", + "components": ["v_0", "v_1"]}} + h = _equation(ncls, "φ")["hamiltonian"] + assert h["H"] == "grad_φ_0*v_0 + grad_φ_1*v_1" + assert h["dH"] == {"φ": ["v_0", "v_1"]} + + +def test_vof_is_conservative_transport(split): + vof = _equation(_sub(split, "vof"), "χh") + assert vof["advection"]["f"] == ["χh*v_0", "χh*v_1"] + + +def test_vof_starts_from_the_smoothed_level_set(split): + vof = _sub(split, "vof") + assert vof["initial"]["reads"] == ["φ"] + assert vof["regularized"] == {"heaviside": {"by": "smoothed_heaviside", "width": "1.5*h"}} + assert _evaluate(vof["initial"]["χh"], φ=0.0, h=0.1) == 0.5 + assert _evaluate(vof["initial"]["χh"], φ=1.0, h=0.1) == 1.0 + + +def test_redistancing_reads_a_copy_of_the_level_set_and_smooths_its_sign(split): + rdls = _sub(split, "rdls") + assert rdls["reads"]["φs"]["of"] == "φ" + assert rdls["regularized"] == {"sign": {"by": "smoothed_sign", "width": "0.75*h"}} + eq = _equation(rdls, "φ") + assert eq["hamiltonian"]["depends_on"] == {"φ": "nonlinear"} # |∇φ| + # the reaction is -S(φs): 1 - 2 H_ε(φs) + assert _evaluate(eq["reaction"]["r"], φs=1.0, h=0.1) == -1.0 + assert _evaluate(eq["reaction"]["r"], φs=-1.0, h=0.1) == 1.0 + + +def test_mass_correction_is_a_regularized_reaction_diffusion(split): + mcorr = _sub(split, "mcorr") + eq = _equation(mcorr, "φc") + assert eq["reaction"]["depends_on"] == {"φc": "nonlinear"} + assert eq["diffusion"]["φc"]["a"] == ["10*h", "10*h"] + # at φ + φc = 0 the reaction is H_ε(0) - χh = 1/2 - χh, its derivative 1/ε + assert _evaluate(eq["reaction"]["r"], φ=0.0, φc=0.0, χh=0.25, h=0.1) == 0.25 + dr = eq["reaction"]["dr"]["φc"] + assert math.isclose(_evaluate(dr, φ=0.0, φc=0.0, χh=0.25, h=0.1), 1 / 0.15) + + +# --- the regularization is Proteus's ---------------------------------------------- + + +def _proteus_H(eps, phi): # equivalent_polynomials.h, gf.H + if phi > eps: + return 1.0 + if phi < -eps: + return 0.0 + if phi == 0.0: + return 0.5 + return 0.5 * (1.0 + phi / eps + math.sin(math.pi * phi / eps) / math.pi) + + +def _proteus_D(eps, phi): # gf.D + if phi > eps or phi < -eps: + return 0.0 + return 0.5 * (1.0 + math.cos(math.pi * phi / eps)) / eps + + +def test_the_smoothed_heaviside_and_its_derivative_are_proteus_gf_H_and_gf_D(): + from ymf.symbolic.adr import to_code + phi, eps = sympy.symbols("phi eps", real=True) + H = to_code(smoothed_heaviside(phi, eps)) + D = to_code(sympy.diff(smoothed_heaviside(phi, eps), phi)) + for e in (0.01, 0.0483, 0.3): + for k in range(-40, 41): + p = 2 * e * k / 40 + assert math.isclose(_evaluate(H, phi=p, eps=e), _proteus_H(e, p), abs_tol=1e-15) + assert math.isclose(_evaluate(D, phi=p, eps=e), _proteus_D(e, p), + rel_tol=1e-13, abs_tol=1e-13) + + +# --- what is refused ----------------------------------------------------------------- + + +def _without(dambreak, sub, key): + doc = copy.deepcopy(dambreak) + entry = doc["solution_paths"]["discretizations"][0] + next(s for s in entry["subproblems"] if s["name"] == sub).pop(key) + return doc + + +def test_a_field_read_without_saying_when_is_refused(dambreak): + with pytest.raises(SymbolicError, match="sub-problem ncls.*declared in its lags"): + split_problem(_without(dambreak, "ncls", "lags"), "level_set_five") + + +def test_an_unregularized_sharp_function_of_an_unknown_is_refused(dambreak): + doc = copy.deepcopy(dambreak) + doc["solution_paths"]["discretizations"][0].pop("regularize") + with pytest.raises(SymbolicError, match="sub-problem mcorr.*regularizes it"): + split_problem(doc, "level_set_five") + + +def test_a_regularization_ymf_does_not_know_is_refused(dambreak): + doc = copy.deepcopy(dambreak) + doc["solution_paths"]["discretizations"][0]["regularize"]["heaviside"]["by"] = "tanh" + with pytest.raises(SymbolicError, match="regularized by 'tanh'"): + split_problem(doc, "level_set_five") + + +def test_only_split_discretizations_and_real_names(dambreak): + with pytest.raises(SymbolicError, match="no discretization named 'nope'"): + split_problem(dambreak, "nope") + doc = copy.deepcopy(dambreak) + doc["solution_paths"]["discretizations"][0].pop("subproblems") + with pytest.raises(SymbolicError, match="not split into sub-problems"): + split_problem(doc, "level_set_five") + + +def test_a_sub_problem_that_only_spins_up_is_carried_unclassified(dambreak): + doc = copy.deepcopy(dambreak) + doc["solution_paths"]["discretizations"][0]["subproblems"].append( + {"name": "warmup", "spin_up": True}) + warmup = _sub(split_problem(doc, "level_set_five"), "warmup") + assert warmup["classified"] is False and warmup["why"] == "it spins the run up" + + +# --- known fields in the language ---------------------------------------------------- + + +def test_a_known_field_is_differentiated_as_a_field(): + space = Space(2, {"c": 0}, {}, fields={"v": 1}) + reaction = adr_form(parse_equation("∂c/∂t = ∇·v", space), space)["equations"][0]["reaction"] + assert reaction["r"] == "-grad_v_0_0 - grad_v_1_1" + + +def test_the_rate_of_a_known_field_has_no_name(): + space = Space(2, {"c": 0}, {}, fields={"v": 1}) + with pytest.raises(SymbolicError, match="time derivative of the known field v_0"): + adr_form(parse_equation("∂c/∂t = ∂v_0/∂t", space), space) + + +def test_a_name_is_an_unknown_or_a_known_field_not_both(): + with pytest.raises(SymbolicError, match="both an unknown and a known field"): + Space(2, {"v": 1}, {}, fields={"v": 1}) + + +# --- boundary conditions, per sub-problem ----------------------------------------------- + + +def test_each_sub_problem_imposes_the_conditions_on_what_it_solves(split): + flow = [(c["type"], c["variable"], c["region"], c.get("set")) + for c in _sub(split, "flow")["boundary_conditions"]] + assert flow == [("slip", "v", "walls", "free_slip"), + ("dirichlet", "p", "top", "atmosphere"), + ("open", "v", "top", "atmosphere")] + slip = _sub(split, "flow")["boundary_conditions"][0] + assert slip["components"] == ["v_0", "v_1"] and slip["value"] == "0.0" + assert _sub(split, "ncls")["boundary_conditions"] == [] + + +def test_a_condition_on_a_coefficient_goes_to_its_carrier(split): + (inflow,) = _sub(split, "vof")["boundary_conditions"] + assert inflow == {"type": "inflow", "variable": "χ", "region": "top", + "where": "(abs(y - 1.8) <= tol)", "carried_by": "χh", + "set": "atmosphere", "components": ["χh"], "value": "1.0"} + + +def test_a_condition_reaches_every_sub_problem_that_solves_for_its_variable(dambreak): + doc = copy.deepcopy(dambreak) + doc["model"]["strong_form"]["boundary_conditions"] = [ + {"region": "top", "variable": "φ", "type": "dirichlet", "value": 1.0}, + {"region": "top", "variable": "v", "type": "inflow", "formula": "(0, -t)"}, + {"region": "walls", "variable": "p", "type": "characteristic", + "series": [[0.0, 0.0], [1.0, 2.0]]}] + split = split_problem(doc, "level_set_five") + on_phi = [s["name"] for s in split["subproblems"] + if any(c["variable"] == "φ" for c in s["boundary_conditions"])] + assert on_phi == ["ncls", "rdls"] # not mcorr, which only reads φ + flow = {(c["variable"], c["region"]): c for c in _sub(split, "flow")["boundary_conditions"]} + assert flow[("v", "top")]["value"] == {"v_0": "0", "v_1": "-t"} # replaces open + assert flow[("p", "walls")]["series"] == [[0.0, 0.0], [1.0, 2.0]] diff --git a/tests/test_state_coefficients.py b/tests/test_state_coefficients.py new file mode 100644 index 0000000..8fb02cb --- /dev/null +++ b/tests/test_state_coefficients.py @@ -0,0 +1,196 @@ +"""Coefficients that are functions of the state: a formula may name the unknowns. + +A density ρw (1 - χ) + ρa χ, a conductivity 1 + T², a friction law in h and +q: each is a coefficient whose formula names an unknown. It is substituted +into the equations like any other coefficient, so the ADR form carries its +dependence on the unknowns, and it is not a number, so nothing that wants a +number (a scale, a dimensionless number, boundary data) may lean on it. +""" + +from pathlib import Path + +import pytest + +sympy = pytest.importorskip("sympy") +pytest.importorskip("strictyaml") + +from ymf.schema import load_ymf # noqa: E402 +from ymf.symbolic import SymbolicError, adr_problem, state_coefficients # noqa: E402 +from ymf.symbolic.problem import coefficient_numbers # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + + +def _heat(**coefficients): + """The heat example without its exact solution (it is not exact once κ + depends on T), with coefficients replaced or added.""" + doc = load_ymf(EXAMPLES / "heat_equation.ymf") + doc["solution_paths"]["analytical"] = [] + doc["model"]["strong_form"]["coefficients"].update(coefficients) + return doc + + +def _equation(problem, component): + return next(e for e in problem["adr"]["equations"] if e["component"] == component) + + +def test_a_coefficient_of_the_state_is_the_inlined_expression(): + named = adr_problem(_heat(κ={"formula": "1e-3 (1 + T²)"})) + inlined = _heat() + inlined["model"]["strong_form"]["equations"] = ["∂T/∂t = ∇·(1e-3 (1 + T²)∇T) in Ω"] + assert named["adr"] == adr_problem(inlined)["adr"] + diffusion = _equation(named, "T")["diffusion"]["T"] + assert diffusion["depends_on"] == {"T": "nonlinear"} + assert diffusion["da"]["T"] == ["0.002*T", "0.002*T"] + + +def test_a_chain_of_coefficients_reaches_the_state(): + doc = _heat(κ={"formula": "κ0 c"}, c={"formula": "1 + T²"}, κ0={"value": 1.0e-3}) + inlined = _heat() + inlined["model"]["strong_form"]["equations"] = ["∂T/∂t = ∇·(1e-3 (1 + T²)∇T) in Ω"] + assert adr_problem(doc)["adr"] == adr_problem(inlined)["adr"] + + +def test_a_state_dependent_storage_term_is_a_nonlinear_mass(): + doc = _heat(ρc={"formula": "1 + T"}) + doc["model"]["strong_form"]["equations"] = ["∂(ρc T)/∂t = ∇·(κ∇T) in Ω"] + mass = _equation(adr_problem(doc), "T")["mass"] + assert mass["depends_on"] == {"T": "nonlinear"} + + +def test_a_coefficient_may_name_a_vector_unknown(): + doc = load_ymf(EXAMPLES / "navier_stokes" / "plane_poiseuille.ymf") + doc["solution_paths"]["analytical"] = [] + doc["vvuq"] = {} + strong = doc["model"]["strong_form"] + strong["coefficients"]["μ"] = {"formula": "1 + v·v"} + # μΔv is the divergence form only while μ is constant + strong["equations"][0] = "ρ ∂v/∂t + ∇·(ρ v⊗v) − ∇·(μ∇v) + ∇p = f in Ω" + momentum = _equation(adr_problem(doc), "v_0") + assert set(momentum["diffusion"]["v_0"]["depends_on"]) == {"v_0", "v_1"} + + +def test_a_cycle_is_still_an_error(): + with pytest.raises(SymbolicError, match="could not be resolved"): + adr_problem(_heat(κ={"formula": "1e-3 c"}, c={"formula": "κ T"})) + + +def test_a_coefficient_may_not_take_the_name_of_an_unknown(): + with pytest.raises(SymbolicError, match="name of an unknown"): + adr_problem(_heat(T={"value": 1.0})) + + +def test_a_coefficient_may_not_depend_on_a_derivative_yet(): + with pytest.raises(SymbolicError, match="not \\(yet\\) on their derivatives"): + adr_problem(_heat(κ={"formula": "1e-3 (1 + ∂T/∂x)"})) + + +@pytest.mark.parametrize("where", ["boundary", "initial", "exact"]) +def test_point_data_may_not_lean_on_the_state(where): + doc = _heat(k={"formula": "300 + T"}) + strong = doc["model"]["strong_form"] + if where == "boundary": + strong["boundary_conditions"][0] = {"region": "walls", "variable": "T", + "type": "dirichlet", "formula": "k"} + elif where == "initial": + strong["initial_conditions"][0]["formula"] = "k" + else: + doc["solution_paths"]["analytical"] = [ + {"name": "wrong", "provenance": "human_specified", "method": "closed_form", + "solution": {"formula": "T = k"}}] + with pytest.raises(SymbolicError, match=r"depends on the unknowns \(T\)"): + adr_problem(doc) + + +def test_the_exact_solution_of_the_heat_example_is_refused_once_kappa_depends_on_t(): + """Its formula names κ; with κ = κ(T) it is no longer a function of x and t.""" + doc = load_ymf(EXAMPLES / "heat_equation.ymf") + doc["model"]["strong_form"]["coefficients"]["κ"] = {"formula": "1e-3 (1 + T²)"} + with pytest.raises(SymbolicError, match="analytical solution"): + adr_problem(doc) + + +def test_a_state_coefficient_is_not_a_number(): + spec = {"κ": {"formula": "1e-3 (1 + T²)"}, "κ0": {"value": 2.0}} + assert coefficient_numbers(spec) == {"κ0": 2.0} + + +def test_state_coefficients_say_what_they_depend_on(): + strong = _heat(κ={"formula": "κ0 c"}, c={"formula": "1 + T²"}, κ0={"value": 1e-3}) + assert state_coefficients(strong["model"]["strong_form"], 2) == {"κ": ["T"], "c": ["T"]} + + +def test_a_scale_built_on_a_state_coefficient_is_reported(): + pytest.importorskip("pint") + from ymf import check_units + doc = _heat(κ={"formula": "1e-3 (1 + T²)", "units": "m2/s"}) + # the example derives its time scale from κ + issues = [i for i in check_units(doc).issues + if i.location == "characteristic_scales[time]"] + assert len(issues) == 1 and issues[0].severity == "warning" + assert "coefficients[κ] depends on the unknowns (T)" in issues[0].message + assert "constant coefficient" in issues[0].message + + +def test_a_scale_built_on_a_constant_coefficient_is_not_reported(): + pytest.importorskip("pint") + from ymf import check_units + doc = _heat(c={"formula": "1 + T²"}) + assert not [i for i in check_units(doc).issues if "depends on the unknowns" in i.message] + + +def test_a_model_declares_the_coefficient_and_a_case_makes_it_state_dependent(tmp_path): + (tmp_path / "model.ymf").write_text( + (EXAMPLES / "heat_equation.ymf").read_text(encoding="utf-8") + .replace("κ: {value: 1.0e-3, units: m2/s}", "κ: {units: m2/s}") + .replace('time: {derived: true, formula: "1/(2π²κ)", units: s}', + "time: {value: 1.0, units: s}"), + encoding="utf-8") + (tmp_path / "case.ymf").write_text( + "extends: model.ymf\n" + "remove:\n" + ' - "solution_paths.analytical[separation_of_variables]"\n' + "model:\n" + " strong_form:\n" + " coefficients:\n" + ' κ: {formula: "1e-3 (1 + T²)"}\n', + encoding="utf-8") + assert "κ: {units: m2/s}" in (tmp_path / "model.ymf").read_text(encoding="utf-8") + doc = load_ymf(tmp_path / "case.ymf") + assert doc["model"]["strong_form"]["coefficients"]["κ"] == {"units": "m2/s", + "formula": "1e-3 (1 + T²)"} + diffusion = _equation(adr_problem(doc), "T")["diffusion"]["T"] + assert diffusion["depends_on"] == {"T": "nonlinear"} + + +# --- the names, kept beside the substituted ADR form -------------------------- + + +def test_a_state_coefficient_keeps_its_name_in_the_plain_data(): + problem = adr_problem(load_ymf(EXAMPLES / "nonlinear_diffusion.ymf")) + assert problem["coefficients"] == { + "k": {"code": "u**2 + 1", "depends_on": {"u": "nonlinear"}, "d": {"u": "2*u"}}} + + +def test_a_problem_without_one_reads_as_before(): + assert "coefficients" not in adr_problem(load_ymf(EXAMPLES / "poisson.ymf")) + + +def test_every_link_of_a_chain_is_named(): + named = adr_problem(_heat(κ={"formula": "κ0 c"}, c={"formula": "1 + T²"}, + κ0={"value": 1e-3}))["coefficients"] + assert set(named) == {"κ", "c"} + assert named["c"]["d"] == {"T": "2*T"} + + +def test_a_vector_coefficient_of_the_state(): + named = adr_problem(_heat(w={"formula": "(T², T)"}))["coefficients"] + assert named["w"]["code"] == ["T**2", "T"] + assert named["w"]["depends_on"] == {"T": "nonlinear"} + assert named["w"]["d"] == {"T": ["2*T", "1"]} + + +def test_a_sharp_coefficient_is_named_without_a_derivative(): + named = adr_problem(_heat(χ={"formula": "heaviside(T − 300)"}))["coefficients"] + assert named["χ"]["sharp"] is True + assert "d" not in named["χ"] diff --git a/tests/test_subproblems.py b/tests/test_subproblems.py new file mode 100644 index 0000000..2bf511b --- /dev/null +++ b/tests/test_subproblems.py @@ -0,0 +1,154 @@ +"""A discretization split into sub-problems (C42). + +tests/fixtures/two_phase holds the two-phase level-set model and the dam +break split the way proteus solves it: flow, VOF, level set, redistancing +and mass correction. ymf checks the split's structure and names on load; a +solver's dispatcher is what maps each sub-problem onto a module. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.discretizations import subproblem_solver # noqa: E402 +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 +from ymf.state import reached # noqa: E402 + +TWO_PHASE = Path(__file__).parent / "fixtures" / "two_phase" + + +@pytest.fixture +def dambreak(): + return load_ymf(TWO_PHASE / "dambreak.ymf") + + +def _split(doc): + return doc["solution_paths"]["discretizations"][0] + + +def _sub(doc, name): + return next(s for s in _split(doc)["subproblems"] if s["name"] == name) + + +def test_the_model_alone_is_abstract(): + model = load_ymf(TWO_PHASE / "two_phase_level_set.ymf") + assert model["ymf_state"] == "abstract" + + +def test_the_split_dam_break_is_discretized(dambreak): + state, _ = reached(dambreak) + assert state == "discretized" + assert [s["name"] for s in _split(dambreak)["subproblems"]] == [ + "flow", "vof", "ncls", "rdls", "mcorr"] + + +def test_a_sub_problem_refines_the_entry_solver(dambreak): + entry = _split(dambreak) + assert subproblem_solver(entry, _sub(dambreak, "flow"))["max_iterations"] == 100 + assert subproblem_solver(entry, _sub(dambreak, "ncls"))["max_iterations"] == 50 + assert subproblem_solver(entry, _sub(dambreak, "ncls"))["nonlinear_solver"] == "newton" + + +def test_petsc_options_merge_key_by_key(): + entry = {"solver": {"linear_solver": "petsc", + "petsc_options": {"ksp_type": "gmres", "pc_type": "asm"}}} + sub = {"solver": {"petsc_options": {"pc_type": "lu"}}} + assert subproblem_solver(entry, sub)["petsc_options"] == {"ksp_type": "gmres", + "pc_type": "lu"} + + +def test_the_split_is_part_of_the_output_key(dambreak): + from ymf.closure import planned_outputs + before = [k for k, _, _ in planned_outputs(dambreak)] + _sub(dambreak, "flow")["solver"]["max_iterations"] = 99 + after = [k for k, _, _ in planned_outputs(dambreak)] + assert before != after + + +def test_the_whole_model_is_not_one_adr_system(dambreak): + """v·∇φ is a Hamiltonian in another unknown: classified only per sub-problem.""" + pytest.importorskip("sympy") + from ymf.symbolic import SymbolicError, adr_problem + with pytest.raises(SymbolicError): + adr_problem(dambreak) + + +def _broken(doc, change): + doc = copy.deepcopy(doc) + change(doc) + return doc + + +def _set(path, value): + def change(doc): + node = doc + for key in path[:-1]: + node = node[key] if not isinstance(key, tuple) else next( + s for s in node if s["name"] == key[0]) + node[path[-1]] = value + return change + + +SPLIT = ("solution_paths", "discretizations", ("level_set_five",)) + + +@pytest.mark.parametrize("change, message", [ + (_set(SPLIT + ("subproblems", ("ncls",), "solves"), ["interfaces"]), + "names 'interfaces', which is not a model equation"), + (_set(SPLIT + ("subproblems", ("ncls",), "solves"), ["momentum"]), + "momentum is solved by both flow and ncls"), + (_set(SPLIT + ("subproblems", ("ncls",), "solves"), []), + "neither solves a model equation"), + (_set(SPLIT + ("subproblems", ("rdls",), "name"), "flow"), + "flow named more than once"), + (_set(SPLIT + ("subproblems", ("vof",), "introduces"), {"name": "ρ"}), + "ρ is already a name \\(in the model\\)"), + (_set(SPLIT + ("subproblems", ("mcorr",), "introduces"), {"name": "φc", + "correction_of": "χ"}), + "names 'χ', which is not a model unknown"), + (_set(SPLIT + ("subproblems", ("ncls",), "lags"), {"w": "current_step"}), + "names 'w', which is not a model unknown, coefficient or introduced unknown"), + (_set(SPLIT + ("subproblems", ("ncls",), "lags"), {"v": {"via": "rans", "at": "current_step"}}), + "names 'rans', which is not a sub-problem"), + (_set(SPLIT + ("subproblems", ("rdls",), "writes_back"), {"ρ": "replace"}), + "names 'ρ', which is not a model or introduced unknown"), + (_set(SPLIT + ("binds",), {"g": {"value": 1.0}}), + "g is already a name in the model"), + (_set(SPLIT + ("regularize",), {"sqrt": {"by": "smooth"}}), + "only the sharp functions heaviside, sign are regularized"), +]) +def test_a_split_that_does_not_add_up_is_refused(dambreak, change, message): + with pytest.raises(YmfReferenceError, match=message): + resolve(_broken(dambreak, change)) + + +def test_every_model_equation_is_solved_by_some_sub_problem(dambreak): + def drop_ncls(doc): + _split(doc)["subproblems"].remove(_sub(doc, "ncls")) + _sub(doc, "rdls")["lags"]["φs"]["via"] = "vof" # it read φ as ncls left it + doc = _broken(dambreak, drop_ncls) + with pytest.raises(YmfReferenceError, match="no sub-problem solves interface"): + resolve(doc) + + +def test_splitting_needs_named_equations(dambreak): + def unname(doc): + strong = doc["model"]["strong_form"] + strong["equations"] = [e["equation"] for e in strong["equations"]] + doc = _broken(dambreak, unname) + with pytest.raises(YmfReferenceError, match="name them"): + resolve(doc) + + +def test_the_file_level_errors_reach_the_loader(tmp_path): + text = (TWO_PHASE / "dambreak.ymf").read_text(encoding="utf-8") + (tmp_path / "two_phase_level_set.ymf").write_text( + (TWO_PHASE / "two_phase_level_set.ymf").read_text(encoding="utf-8"), encoding="utf-8") + (tmp_path / "dambreak.ymf").write_text( + text.replace("solves: [interface]", "solves: [interfaces]"), encoding="utf-8") + with pytest.raises(YmfReferenceError, match="not a model equation"): + load_ymf(tmp_path / "dambreak.ymf") diff --git a/tests/test_symbolic.py b/tests/test_symbolic.py index bb0386a..dd87086 100644 --- a/tests/test_symbolic.py +++ b/tests/test_symbolic.py @@ -313,10 +313,14 @@ def test_a_piecewise_linear_region_is_carried_with_its_bounding_box(): yaml.safe_dump(problem) # still plain data -def test_periodic_faces_need_a_box(): - doc = _poisson_with({"region": "all_boundaries", "variable": "u", "type": "periodic"}) +def test_periodic_faces_of_a_region_must_be_translations(): + """A region may be periodic (C27 Phase 5, D3) when its two faces match; + the L-shape's x faces do not (lengths 2 and 1).""" + doc = _poisson_with({"region": "ends", "variable": "u", "type": "periodic"}) + doc["model"]["strong_form"]["boundary_regions"].append( + {"name": "ends", "where": "x = 0 or x = 2"}) doc["model"]["strong_form"]["domain"] = L_SHAPE - with pytest.raises(SymbolicError, match="periodic faces are a box's"): + with pytest.raises(SymbolicError, match="faces of axis x do not match: 2 of 2"): adr_problem(doc) diff --git a/tests/test_tag_regions.py b/tests/test_tag_regions.py new file mode 100644 index 0000000..b246d74 --- /dev/null +++ b/tests/test_tag_regions.py @@ -0,0 +1,121 @@ +"""Boundary regions named by the mesh's boundary tags (C27 Phase 5, D1). + +A region is ``where`` (a predicate on the point) or ``tags`` (ints, or names +the domain's ``boundary_tags`` resolve), exactly one. A tag region reaches a +solver as ``flags``, tested against the boundary flag of the point, and has +no ``where`` key: a solver that does not know ``flags`` fails on it rather +than reading a missing predicate as the whole boundary. +""" + +import copy +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf.references import YmfReferenceError, resolve # noqa: E402 +from ymf.schema import load_ymf # noqa: E402 + +EXAMPLES = Path(__file__).parent.parent / "examples" + +L_SHAPE = { + "kind": "piecewise_linear", "dimension": 2, + "vertices": [[0.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [1.0, 2.0], [0.0, 2.0]], + "segments": [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]], + "segment_flags": [1, 2, 3, 3, 4, 5], + "boundary_tags": {"bottom": 1, "right": 2, "notch": 3, "top": 4, "left": 5}, +} + + +@pytest.fixture +def lshape(): + doc = load_ymf(EXAMPLES / "poisson.ymf") + strong = doc["model"]["strong_form"] + strong["domain"] = copy.deepcopy(L_SHAPE) + strong["boundary_regions"] = [{"name": "dry", "tags": ["bottom", "left"]}, + {"name": "wet", "tags": [2, "notch", "top"]}] + strong["boundary_conditions"] = [ + {"region": "dry", "variable": "u", "type": "dirichlet", "value": 0.0}, + {"region": "wet", "variable": "u", "type": "neumann", "formula": "x"}] + doc["solution_paths"]["analytical"] = [] + doc.pop("vvuq", None) + for d in doc["solution_paths"]["discretizations"]: + d["mesh"] = {"h": 0.25} + return doc + + +def test_a_tag_region_reaches_the_solver_as_flags(lshape): + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + resolve(lshape) + problem = adr_problem(lshape) + (dry,) = problem["dirichlet"] + assert dry["flags"] == [1, 5] and "where" not in dry + (wet,) = problem["neumann"] + assert wet["flags"] == [2, 3, 4] and wet["value"] == "x" + + +def test_where_and_tags_mix(lshape): + pytest.importorskip("sympy") + from ymf.symbolic import adr_problem + lshape["model"]["strong_form"]["boundary_regions"][1] = {"name": "wet", "where": "x = 2"} + (wet,) = adr_problem(lshape)["neumann"] + assert wet["where"] == "(abs(x - 2) <= tol)" and "flags" not in wet + + +def _with_region(doc, region): + doc = copy.deepcopy(doc) + doc["model"]["strong_form"]["boundary_regions"][0] = region + return doc + + +@pytest.mark.parametrize("region, message", [ + ({"name": "dry"}, "exactly one; this gives neither"), + ({"name": "dry", "where": "x = 0", "tags": [5]}, "exactly one; this gives where and tags"), + ({"name": "dry", "tags": ["floor"]}, "tag 'floor' is not one of the domain's boundary_tags"), +]) +def test_a_region_says_how_it_is_found(lshape, region, message): + with pytest.raises(YmfReferenceError, match=message): + resolve(_with_region(lshape, region)) + + +def test_a_box_has_no_tags(): + doc = load_ymf(EXAMPLES / "poisson.ymf") + doc["model"]["strong_form"]["boundary_regions"][0] = {"name": "all_boundaries", + "tags": [1]} + with pytest.raises(YmfReferenceError, match="a box carries no boundary tags"): + resolve(doc) + + +def test_a_periodic_region_is_found_by_position(lshape): + lshape["model"]["strong_form"]["boundary_conditions"][0]["type"] = "periodic" + lshape["model"]["strong_form"]["boundary_conditions"][0].pop("value") + with pytest.raises(YmfReferenceError, match="needs where:, not tags:"): + resolve(lshape) + + +def test_a_tag_region_loads_from_a_file(tmp_path): + text = (EXAMPLES / "poisson.ymf").read_text(encoding="utf-8") + text = text.replace(''' domain: "Ω = [0, 1] × [0, 1]"''', ''' domain: + kind: piecewise_linear + dimension: 2 + vertices: [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + segments: [[0, 1], [1, 2], [2, 3], [3, 0]] + segment_flags: [1, 1, 2, 1] + boundary_tags: {sides: 1, lid: 2}''').replace(''' where: "∂Ω"''', ''' tags: [sides, lid]''') + text = text.replace("mesh: {cells: [4, 8, 16]}", "mesh: {h: 0.25}") + (tmp_path / "p.ymf").write_text(text, encoding="utf-8") + region = load_ymf(tmp_path / "p.ymf")["model"]["strong_form"]["boundary_regions"][0] + assert region == {"name": "all_boundaries", "tags": ["sides", "lid"]} + + +def test_split_problem_uses_flags_too(lshape): + pytest.importorskip("sympy") + from ymf.symbolic import split_problem + lshape["model"]["strong_form"]["equations"] = [ + {"name": "poisson", "for": "u", "equation": "-Δu = f in Ω"}] + entry = lshape["solution_paths"]["discretizations"][0] + entry["subproblems"] = [{"name": "one", "solves": ["poisson"]}] + (sub,) = split_problem(lshape, entry["name"])["subproblems"] + assert [c.get("flags") for c in sub["boundary_conditions"]] == [[1, 5], [2, 3, 4]] diff --git a/tests/test_units.py b/tests/test_units.py index 1a0ce50..9c58b66 100644 --- a/tests/test_units.py +++ b/tests/test_units.py @@ -47,6 +47,8 @@ def test_normalize_unit_string(raw, expected_normalized): def test_to_pint_unit_none_and_dimensionless(): assert to_pint_unit(None).dimensionless assert to_pint_unit("dimensionless").dimensionless + # "1", the usual way to write "no units"; pint before 0.25 parses it to an int + assert to_pint_unit("1").dimensionless assert to_pint_unit("").dimensionless