From 9c0ccbc73416ee7957af4768e1125fe2e571f8d2 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 09:26:38 +0100 Subject: [PATCH 1/5] Closure writes out the split and names the coefficients of the state MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The alias expansion (#9) made a closed spec say what was solved for its boundary conditions. Two more things were still left to a reader's defaults: a split discretization's sub-problems (a lag written `φ: previous_step`, a coefficient read through the sub-problem that carries it, a sub-problem's solver refining the entry's, its regularization replacing the entry's function by function, its `for` defaulting to what it introduces), and which coefficients are functions of the state. ymf.explicit writes them out, and closure.write and closure.input_slice use it in place of the bare expansion: - each coefficient whose formula names an unknown, directly or through another coefficient, in the map form with `depends_on` (ρ, μ, χ: [φ] in the dam break; τb, Fd, Fr, β: [h, q] in Torino); - each sub-problem with its `for`, merged `solver` and `regularize`, and every lag as {of, via, at, carried_by} (flow's χ: {via: vof, at: previous_step, carried_by: χh}). The names are read without sympy: the notation's tokenizer moves to ymf.names (ymf.symbolic.notation re-exports it) with the parser's rule for a run of letters (a declared name, else the product of its letters; a function call is not a name), so the core closes a spec as the front end does. It agrees with ymf.symbolic.problem.state_coefficients on every coefficient that resolves (14 in examples and fixtures) and also names the abstract model's ρ and μ, which the symbolic reader leaves out for lack of values. A test checks the name reader's function list is the parser's. `depends_on` and `carried_by` are checked on load (a false one is refused) and `depends_on` does not count towards a key. Measured, closing every example and fixture with closure.write on cekees/ymf main b77262a and on this commit: 9 specs without coefficients of the state or sub-problems close to byte-identical documents with the same keys; examples/nonlinear_diffusion gains `depends_on: [u]` and keeps its 6 keys; the dam break, lake at rest and Torino re-key (the split is written out; no solver has run a split spec, Proteus refuses one), and each closed spec reloads with its spec's keys. split_problem gives equal results for the dam break spec and its closed form. ymf suite 845 passed, 2 skipped (was 827; 16 new in tests/test_explicit.py, 2 core-import cases), macOS arm64, petsc-dev env. Co-Authored-By: Claude Opus 5.5 --- docs/changes.md | 15 +++ docs/problem-specification.md | 25 +++- src/ymf/closure.py | 28 ++-- src/ymf/explicit.py | 160 ++++++++++++++++++++++ src/ymf/names.py | 174 ++++++++++++++++++++++++ src/ymf/references.py | 29 ++++ src/ymf/schema.py | 8 +- src/ymf/symbolic/notation.py | 75 +---------- tests/test_core_without_frontend.py | 2 +- tests/test_explicit.py | 200 ++++++++++++++++++++++++++++ 10 files changed, 627 insertions(+), 89 deletions(-) create mode 100644 src/ymf/explicit.py create mode 100644 src/ymf/names.py create mode 100644 tests/test_explicit.py diff --git a/docs/changes.md b/docs/changes.md index f168403..cb03c0a 100644 --- a/docs/changes.md +++ b/docs/changes.md @@ -7,6 +7,21 @@ rediscover why a thing is the way it is. Open gaps are in ## 2026-11 +- **2026-11-05: closure writes out the split and names the coefficients of + the state (C42 round 2).** The same reasoning as the alias expansion: a + closed dam break said `φ: previous_step` and left the sub-problems' + solvers, regularizations and carriers to `split_problem`'s defaults, and + said nothing of which coefficients follow the state. `ymf.explicit` + (closure.write, closure.input_slice) now writes each sub-problem's + `for`, merged `solver` and `regularize`, and its lags as `{of, via, at, + carried_by}`, and gives every coefficient whose formula names an unknown + `depends_on`. The names are read without sympy (`ymf.names`, the + notation's tokenizer moved out of `ymf.symbolic.notation`, with the + parser's rule for a run of letters), so the core closes a spec as the + front end does. `depends_on` is checked on load and left out of the key; + split specs re-key (the dam break and shallow-water fixtures; no solver + ran them), every other spec keeps its keys and, without coefficients of + the state, closes to the same document. - **2026-11-04: aliases expand into the abstract forms, and closure records the expansion.** cekees's review note: `slip`, `open`, `inflow`, `characteristic` and condition sets had their meaning only in code diff --git a/docs/problem-specification.md b/docs/problem-specification.md index aeed823..3e63fc1 100644 --- a/docs/problem-specification.md +++ b/docs/problem-specification.md @@ -413,6 +413,14 @@ 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. +A closed spec (an archive) names them too: closure writes each coefficient +of the state in the map form with `depends_on`, the unknowns it follows, +directly or through another coefficient (`ρ: {formula: "ρw (1 − χ) + ρa χ", +depends_on: [φ]}`). It is read from the formulas' text (`ymf.names`, no +sympy, so every environment closes a spec alike), checked on load (a +`depends_on` that is not what the formula says is refused), and does not +count towards an output's key: the formula it comes from already does. + 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. @@ -783,8 +791,21 @@ 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`). +replaces the entry's for the functions it names (redistancing smooths its +sign over `0.75 h`, the rest their Heaviside over `1.5 h`). + +**A closed spec writes the split out.** What the rules above leave implicit +is written into each sub-problem when a spec is closed (`ymf.explicit`, +for the reason the boundary-condition aliases are expanded: two readers +must not close one spec differently): its `for` (the unknown it +introduces, when it gives an equation and no `for`); its `solver`, the +entry's refined by its own; its `regularize`, the entry's with its own +replacing function by function; and every lag as `{of, via, at, +carried_by}`, where `carried_by` names the introduced unknown a lagged +coefficient is read as (`χ: {via: vof, at: previous_step, carried_by: +χh}`). A `carried_by` that is not the carrier is refused on load. Written +out or not, the split classifies the same (`split_problem` gives equal +results for both) and keys its outputs alike. **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 diff --git a/src/ymf/closure.py b/src/ymf/closure.py index 972af06..958fba8 100644 --- a/src/ymf/closure.py +++ b/src/ymf/closure.py @@ -163,13 +163,15 @@ def input_slice(doc: Dict[str, Any], discretization: str, realization: Dict[str, Any]) -> Dict[str, Any]: """The part of the input one output depends on (see the module docstring). - Taken from the document with its boundary-condition aliases and condition - sets expanded (ymf.conditions): what is solved is the expansion, so a - spec and the archive it closes into, which holds the expansion, key - their outputs alike; a spec with no shorthand is unchanged by it. + Taken from the document as a closed spec records it (ymf.explicit: the + boundary-condition aliases and condition sets expanded, the sub-problems + written out), so a spec and the archive it closes into key their outputs + alike; a spec with no shorthand is unchanged by it. What closure derives + from the input alone (a coefficient's ``depends_on``) is left out: the + formula it is read from is in the slice already. """ - from ymf.conditions import expand - doc = expand(doc) + from ymf.explicit import explicit, without_derived + doc = without_derived(explicit(doc)) paths = doc.get("solution_paths") or {} entry = _discretization(doc, discretization) return plain({ @@ -320,13 +322,15 @@ def archive_stem(path: str | Path) -> str: def write(path: str | Path, spec: Dict[str, Any], outputs: Dict[str, Any]) -> None: """Write an archive: ``ymf:``, its ``ymf_state``, the realized spec, then ``outputs``. The state is the one the archive reaches (see ymf.state), - whatever its input declared. The realized spec has its boundary-condition - aliases and condition sets expanded (ymf.conditions).""" - from ymf.conditions import expand + whatever its input declared. The realized spec is written out + (ymf.explicit): boundary-condition aliases and condition sets expanded, + the coefficients of the state named with what they follow, every + sub-problem with its solver, regularization and lags in full.""" + from ymf.explicit import explicit from ymf.state import reached - # the closure records what was solved: the conditions expanded, each - # with the alias and the set it came from (ymf.conditions) - document = dict(plain(expand(spec))) + # the closure records what was solved, with nothing left to a reader's + # defaults (ymf.explicit) + document = dict(plain(explicit(spec))) document.pop("outputs", None) document.pop("ymf_state", None) document.pop("generated", None) # it describes the input file, not this one diff --git a/src/ymf/explicit.py b/src/ymf/explicit.py new file mode 100644 index 0000000..44ff9d9 --- /dev/null +++ b/src/ymf/explicit.py @@ -0,0 +1,160 @@ +"""What a closed spec records: the input, with its shorthand written out. + +A spec may leave things to rules: a boundary condition may be an alias, a +sub-problem's solver refines the entry's, its regularization replaces the +entry's function by function, a lag may be a bare ``previous_step``, a +coefficient read through the sub-problem that carries it is that carrier. +Each rule has one meaning, but if it lives only in a reader's code, two +readers can close the same spec differently and the archive cannot show +which. So the closed spec writes each one out (cekees, 2026-11-04, for the +boundary conditions; C42 round 2 for the rest): + +* boundary-condition aliases and condition sets, expanded + (:func:`ymf.conditions.expand`); +* every **coefficient of the state** -- one whose formula names an unknown, + directly or through another coefficient (``ρ = ρw (1 − χ) + ρa χ`` with + ``χ = heaviside(φ)``) -- in the map form, with ``depends_on``: the + unknowns it follows. A solver of one system substitutes such a + coefficient; a kernel generator emits it by name, as one inline function + of the state, and this is where the closed spec says which names those + are; +* every **sub-problem** of a split discretization with what it solves for + (``for``, by default the unknown it introduces), its solver (the entry's, + refined by its own), its regularization (the entry's, function by + function replaced by its own), and each lag as ``{at, via, of, + carried_by}``, where ``carried_by`` names the introduced unknown that + carries a lagged coefficient (``χ`` via VOF is read as ``χh``). + +Applying it twice changes nothing, and a spec with none of the shorthand +(no aliases, no coefficients of the state, no sub-problems) is returned +unchanged. Readers may read either form; what they solve is the same. +""" + +from __future__ import annotations + +import copy +from typing import Any, Dict, List + +__all__ = ["explicit", "without_derived", "state_coefficients", "split_out", "carriers", + "lag_of", "DERIVED_KEYS"] + +#: keys explicit() derives from what the input already says, which therefore +#: do not count towards an output's key (ymf.closure.input_slice) +DERIVED_KEYS = ("depends_on",) + + +def _strong(doc: Dict[str, Any]) -> Dict[str, Any]: + return ((doc.get("model") or {}).get("strong_form") or {}) + + +def _unknowns(strong: Dict[str, Any]) -> List[str]: + return [u if isinstance(u, str) else u.get("name") + for u in strong.get("unknowns") or []] + + +def state_coefficients(doc: Dict[str, Any]) -> Dict[str, List[str]]: + """``{coefficient: [unknown, ...]}``: the coefficients of the state, read + from the formulas' text (ymf.names; no sympy).""" + from ymf.names import state_dependence + strong = _strong(doc) + return state_dependence(strong.get("coefficients") or {}, _unknowns(strong)) + + +def _name_state(doc: Dict[str, Any]) -> Dict[str, Any]: + state = state_coefficients(doc) + if not state: + return doc + doc = dict(doc) + model = dict(doc["model"]) + strong = dict(model["strong_form"]) + coefficients = {} + for name, entry in strong["coefficients"].items(): + if name in state: + entry = dict(entry) if isinstance(entry, dict) else {"formula": str(entry)} + entry["depends_on"] = list(state[name]) + coefficients[name] = entry + strong["coefficients"] = coefficients + model["strong_form"] = strong + doc["model"] = model + return doc + + +def carriers(entry: Dict[str, Any]) -> Dict[str, str]: + """``{model quantity: introduced unknown}`` for one split entry: what each + introduced carrier approximates (``χ`` -> ``χh``).""" + out = {} + for sub in entry.get("subproblems") or []: + new = sub.get("introduces") or {} + if new.get("approximates"): + out[new["approximates"]] = new["name"] + return out + + +def lag_of(name: str, lag: Any, carried: Dict[str, str]) -> Dict[str, Any]: + """One lag written out: ``{of?, via?, at, carried_by?}``.""" + lag = dict(lag) if isinstance(lag, dict) else {"at": lag} + out = {key: lag[key] for key in ("of", "via") if lag.get(key) is not None} + out["at"] = lag["at"] + target = lag.get("of", name) + if target in carried: + out["carried_by"] = carried[target] + elif lag.get("carried_by") is not None: + out["carried_by"] = lag["carried_by"] # refused on load if it is wrong + return out + + +def split_out(entry: Dict[str, Any]) -> Dict[str, Any]: + """One discretization entry with every sub-problem written out (see the + module docstring); an entry without sub-problems is returned as it is.""" + if not entry.get("subproblems"): + return entry + from ymf.discretizations import subproblem_solver + entry = dict(entry) + carried = carriers(entry) + out = [] + for sub in entry["subproblems"]: + sub = dict(sub) + if sub.get("equation") and not sub.get("for") and sub.get("introduces"): + sub["for"] = sub["introduces"]["name"] + if sub.get("lags"): + sub["lags"] = {name: lag_of(name, lag, carried) + for name, lag in sub["lags"].items()} + regularize = dict(entry.get("regularize") or {}) + regularize.update(sub.get("regularize") or {}) + if regularize: + sub["regularize"] = regularize + solver = subproblem_solver(entry, sub) + if solver: + sub["solver"] = solver + out.append(sub) + entry["subproblems"] = out + return entry + + +def explicit(doc: Dict[str, Any]) -> Dict[str, Any]: + """``doc`` as a closed spec records it (see the module docstring).""" + from ymf.conditions import expand + doc = _name_state(expand(doc)) + entries = (doc.get("solution_paths") or {}).get("discretizations") or [] + if any(isinstance(e, dict) and e.get("subproblems") for e in entries): + doc = dict(doc) + paths = dict(doc["solution_paths"]) + paths["discretizations"] = [split_out(e) if isinstance(e, dict) else e + for e in entries] + doc["solution_paths"] = paths + return doc + + +def without_derived(doc: Dict[str, Any]) -> Dict[str, Any]: + """``doc`` with the keys :func:`explicit` derives from the rest removed.""" + strong = _strong(doc) + coefficients = strong.get("coefficients") or {} + if not any(isinstance(e, dict) and any(k in e for k in DERIVED_KEYS) + for e in coefficients.values()): + return doc + doc = copy.deepcopy(doc) + for entry in doc["model"]["strong_form"]["coefficients"].values(): + if isinstance(entry, dict): + for key in DERIVED_KEYS: + entry.pop(key, None) + return doc diff --git a/src/ymf/names.py b/src/ymf/names.py new file mode 100644 index 0000000..24aa60a --- /dev/null +++ b/src/ymf/names.py @@ -0,0 +1,174 @@ +"""The names a formula uses, read without sympy. + +The notation (``ymf.symbolic.notation``) parses a formula into sympy; this +module holds the part of it that needs no sympy -- the tokenizer, and the +rule by which a run of letters is a name -- so that the core (closure, +references) can say which names a formula uses in an environment without +the ``symbolic`` extra, and get the same answer the parser would. + +The rule, as ``Parser._name`` applies it: a run of letters followed by +``(`` that is a function of the notation is that function; a declared name +(a coefficient, an unknown, a vector unknown's component ``v_0``, a +coordinate, ``t``, ``π``) is that name; any other run of letters is the +product of its letters, each of which must be a name. +""" + +from __future__ import annotations + +import re +import unicodedata +from typing import Dict, Iterable, List, Set, Tuple + +__all__ = ["NotationError", "tokens", "FUNCTION_NAMES", "names_in", "state_dependence"] + + +class NotationError(ValueError): + """A string in a spec does not follow the notation.""" + + +_SUPERSCRIPTS = {"²": 2, "³": 3, "⁴": 4} +_ALIASES = {"−": "-", "⋅": "·", "∙": "·", "–": "-", "∆": "Δ"} +_OPERATORS = ["**", "^", "+", "-", "*", "/", "(", ")", "[", "]", ",", "=", + "·", "⊗", "∇", "Δ", "∂", "×", "∈", ";"] + list(_SUPERSCRIPTS) +_NUMBER = re.compile(r"\d+(\.\d*)?([eE][+-]?\d+)?|\.\d+([eE][+-]?\d+)?") + + +def _is_letter(ch: str) -> bool: + # Δ is a Greek capital letter to Unicode, but here it is the Laplacian + return ch != "Δ" and (ch == "_" or unicodedata.category(ch).startswith("L")) + + +def tokens(text: str) -> List[Tuple[str, str, int]]: + """(kind, value, position) for each token; kinds NUM, NAME, OP, NL.""" + out = [] + i = 0 + while i < len(text): + ch = _ALIASES.get(text[i], text[i]) + if ch == "\n": + out.append(("NL", "\n", i)); i += 1; continue + if ch.isspace(): + i += 1; continue + m = _NUMBER.match(text, i) + if m and (ch in "0123456789" or ch == "."): + out.append(("NUM", m.group(0), i)); i = m.end(); continue + if _is_letter(ch): + j = i + 1 + # ASCII digits only: '²'.isdigit() is True, and ² is an operator + while j < len(text) and (_is_letter(text[j]) or text[j] in "0123456789"): + j += 1 + out.append(("NAME", text[i:j], i)); i = j; continue + for op in _OPERATORS: + if text.startswith(op, i) or ch == op: + out.append(("OP", op if text.startswith(op, i) else ch, i)) + i += len(op) if text.startswith(op, i) else 1 + break + else: + # The `else` of the `for` above: no operator matched. + if _SUBSCRIPTS_START <= ch <= _SUBSCRIPTS_END: + # The superscripts ² ³ ⁴ are powers, so a subscript looks as + # if it ought to work too. It does not, and the fix is to + # spell the index into the name and declare it. + run = ch + while i + len(run) < len(text) and \ + _SUBSCRIPTS_START <= text[i + len(run)] <= _SUBSCRIPTS_END: + run += text[i + len(run)] + stem = out[-1][1] if out and out[-1][0] == "NAME" else "" + plain = "".join(_SUBSCRIPT_PLAIN.get(c, c) for c in run) + raise NotationError( + "%r at %d in %r is a subscript, and the notation has none; " + "write %r and declare it (a coefficient or an unknown)" + % (run, i, text, stem + plain)) + raise NotationError("unexpected character %r at %d in %r" % (text[i], i, text)) + return out + + +#: Unicode subscripts, which the notation does not have. Reported on their +#: own rather than as "unexpected character", because the superscripts +#: ² ³ ⁴ *are* notation (powers), so a subscript is a reasonable guess. +_SUBSCRIPTS_START, _SUBSCRIPTS_END = "\u2080", "\u209c" +_SUBSCRIPT_PLAIN = {chr(0x2080 + d): str(d) for d in range(10)} +_SUBSCRIPT_PLAIN.update({ + "\u2090": "a", "\u2091": "e", "\u2092": "o", "\u2093": "x", + "\u2095": "h", "\u2096": "k", "\u2097": "l", "\u2098": "m", + "\u2099": "n", "\u209a": "p", "\u209b": "s", "\u209c": "t", +}) + + +#: the notation's function names (``ymf.symbolic.notation``'s ``_FUNCTIONS`` +#: and the parser's own; a test checks the two agree) +FUNCTION_NAMES = frozenset([ + "sin", "cos", "tan", "exp", "log", "sqrt", "sinh", "cosh", "tanh", "abs", + "heaviside", "sign", "min", "max", + "div", "grad", "lap", "dt", "dot", "outer", "vec", "norm", "sym", "dx", "dy", "dz", +]) + +#: names every formula may use +BUILTIN_NAMES = frozenset(["x", "y", "z", "t", "π", "pi"]) + + +def names_in(text: str, declared: Iterable[str]) -> Set[str]: + """The declared names ``text`` uses, as the parser reads it. + + A function call is not a name; a run of letters that is not declared + counts for each of its letters (``ρU`` is ``ρ·U``). Names that are not + declared are left out: a formula that uses one is an error the parser + reports, not this. + """ + declared = set(declared) + out: Set[str] = set() + toks = [t for t in tokens(str(text)) if t[0] != "NL"] + for i, (kind, value, _) in enumerate(toks): + if kind != "NAME": + continue + following = toks[i + 1][1] if i + 1 < len(toks) else "" + if following == "(" and value in FUNCTION_NAMES: + continue + if value in declared: + out.add(value) + elif value not in FUNCTION_NAMES: + out.update(ch for ch in value if ch in declared) + return out + + +def _component_of(name: str, unknowns: Iterable[str]) -> str: + stem, _, index = name.rpartition("_") + return stem if stem in unknowns and index.isdigit() else name + + +def state_dependence(coefficients: Dict[str, object], unknowns: Iterable[str] + ) -> Dict[str, List[str]]: + """``{coefficient: [unknown, ...]}`` for each coefficient whose formula + names an unknown (or a component of one), directly or through another + coefficient; sorted, in the order the unknowns are declared. + + Read from the formulas' text alone (no sympy). A coefficient with a + ``value`` is a number, whatever else it has. + """ + unknowns = list(unknowns) + formulas = {} + for name, entry in (coefficients or {}).items(): + if isinstance(entry, dict): + if entry.get("value") is None and entry.get("formula") is not None: + formulas[name] = str(entry["formula"]) + elif isinstance(entry, str): + formulas[name] = entry + components = {"%s_%d" % (u, i) for u in unknowns for i in range(3)} + declared = set(coefficients or {}) | set(unknowns) | components | BUILTIN_NAMES + direct = {name: {_component_of(n, unknowns) for n in names_in(text, declared)} + for name, text in formulas.items()} + out: Dict[str, List[str]] = {} + for name in formulas: + seen, todo, state = set(), [name], set() + while todo: + current = todo.pop() + if current in seen: + continue + seen.add(current) + for used in direct.get(current, ()): + if used in unknowns: + state.add(used) + elif used in direct: + todo.append(used) + if state: + out[name] = [u for u in unknowns if u in state] + return out diff --git a/src/ymf/references.py b/src/ymf/references.py index 8d6adfc..95f85e7 100644 --- a/src/ymf/references.py +++ b/src/ymf/references.py @@ -243,6 +243,15 @@ def _check_split(doc: Dict[str, Any], entry: Dict[str, Any]) -> None: "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 isinstance(lag, dict) and lag.get("carried_by") is not None: + from ymf.explicit import carriers + carrier = carriers(entry).get(lag.get("of", field)) + if lag["carried_by"] != carrier: + raise YmfReferenceError( + "%s.lags.%s.carried_by is %s, but %s; closure writes carried_by, " + "so remove it" % (here, field, lag["carried_by"], + "%s carries it" % carrier if carrier else + "no sub-problem introduces a carrier of it")) if "for" in s: _check(s["for"], sorted(set(unknowns) | set(introduced)), here + ".for", "model or introduced unknown") @@ -414,9 +423,29 @@ def _check_regions(doc: Dict[str, Any]) -> None: "needs where:, not tags:" % (where,)) +def _check_state(doc: Dict[str, Any]) -> None: + """A coefficient's ``depends_on`` (which closure writes, ymf.explicit) must + be what its formula says: a closed spec extended with a changed formula + would otherwise carry a false record.""" + from ymf.explicit import state_coefficients + coefficients = (_problem(doc).get("strong_form") or {}).get("coefficients") or {} + stated = {name: list(entry["depends_on"]) for name, entry in coefficients.items() + if isinstance(entry, dict) and entry.get("depends_on") is not None} + if not stated: + return + derived = state_coefficients(doc) + for name, given in stated.items(): + if sorted(given) != sorted(derived.get(name, [])): + raise YmfReferenceError( + "model.strong_form.coefficients.%s.depends_on says %s, but its formula " + "follows %s; closure writes depends_on, so remove it (or fix the formula)" + % (name, given, derived.get(name) or "no unknown")) + + def resolve(doc: Dict[str, Any]) -> None: """Check every reference, raising :class:`YmfReferenceError` on the first bad one.""" _check_equations(doc) + _check_state(doc) _check_regions(doc) _check_conditions(doc) _check_gauges(doc) diff --git a/src/ymf/schema.py b/src/ymf/schema.py index 77e435a..17ce381 100644 --- a/src/ymf/schema.py +++ b/src/ymf/schema.py @@ -170,8 +170,11 @@ # 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 φ) + # carried_by: the introduced unknown a lagged coefficient is read as + # (χ via VOF is χh); written by closure (ymf.explicit), checked on load Optional("lags"): MapPattern(Str(), LagAtEnum | Map({ - Optional("via"): Str(), "at": LagAtEnum, Optional("of"): Str()})), + Optional("via"): Str(), "at": LagAtEnum, Optional("of"): Str(), + Optional("carried_by"): Str()})), # what it hands back: name -> replace | add | mass_quadrature | ... Optional("writes_back"): MapPattern(Str(), Str()), Optional("pseudo_time"): Enum(["steady"]), @@ -435,6 +438,9 @@ Optional("value"): Float(), Optional("units"): Str(), Optional("formula"): Str(), + # the unknowns a formula follows, directly or through another + # coefficient; written by closure (ymf.explicit), checked on load + Optional("depends_on"): Seq(Str()), } ) diff --git a/src/ymf/symbolic/notation.py b/src/ymf/symbolic/notation.py index 41e6483..26bd426 100644 --- a/src/ymf/symbolic/notation.py +++ b/src/ymf/symbolic/notation.py @@ -47,84 +47,13 @@ from __future__ import annotations import re -import unicodedata from typing import Any, Dict, List, Optional, Tuple import sympy -__all__ = ["NotationError", "Parser", "tokens"] - +from ymf.names import _SUPERSCRIPTS, NotationError, tokens -class NotationError(ValueError): - """A string in a spec does not follow the notation.""" - - -_SUPERSCRIPTS = {"²": 2, "³": 3, "⁴": 4} -_ALIASES = {"−": "-", "⋅": "·", "∙": "·", "–": "-", "∆": "Δ"} -_OPERATORS = ["**", "^", "+", "-", "*", "/", "(", ")", "[", "]", ",", "=", - "·", "⊗", "∇", "Δ", "∂", "×", "∈", ";"] + list(_SUPERSCRIPTS) -_NUMBER = re.compile(r"\d+(\.\d*)?([eE][+-]?\d+)?|\.\d+([eE][+-]?\d+)?") - - -def _is_letter(ch: str) -> bool: - # Δ is a Greek capital letter to Unicode, but here it is the Laplacian - return ch != "Δ" and (ch == "_" or unicodedata.category(ch).startswith("L")) - - -def tokens(text: str) -> List[Tuple[str, str, int]]: - """(kind, value, position) for each token; kinds NUM, NAME, OP, NL.""" - out = [] - i = 0 - while i < len(text): - ch = _ALIASES.get(text[i], text[i]) - if ch == "\n": - out.append(("NL", "\n", i)); i += 1; continue - if ch.isspace(): - i += 1; continue - m = _NUMBER.match(text, i) - if m and (ch in "0123456789" or ch == "."): - out.append(("NUM", m.group(0), i)); i = m.end(); continue - if _is_letter(ch): - j = i + 1 - # ASCII digits only: '²'.isdigit() is True, and ² is an operator - while j < len(text) and (_is_letter(text[j]) or text[j] in "0123456789"): - j += 1 - out.append(("NAME", text[i:j], i)); i = j; continue - for op in _OPERATORS: - if text.startswith(op, i) or ch == op: - out.append(("OP", op if text.startswith(op, i) else ch, i)) - i += len(op) if text.startswith(op, i) else 1 - break - else: - # The `else` of the `for` above: no operator matched. - if _SUBSCRIPTS_START <= ch <= _SUBSCRIPTS_END: - # The superscripts ² ³ ⁴ are powers, so a subscript looks as - # if it ought to work too. It does not, and the fix is to - # spell the index into the name and declare it. - run = ch - while i + len(run) < len(text) and \ - _SUBSCRIPTS_START <= text[i + len(run)] <= _SUBSCRIPTS_END: - run += text[i + len(run)] - stem = out[-1][1] if out and out[-1][0] == "NAME" else "" - plain = "".join(_SUBSCRIPT_PLAIN.get(c, c) for c in run) - raise NotationError( - "%r at %d in %r is a subscript, and the notation has none; " - "write %r and declare it (a coefficient or an unknown)" - % (run, i, text, stem + plain)) - raise NotationError("unexpected character %r at %d in %r" % (text[i], i, text)) - return out - - -#: Unicode subscripts, which the notation does not have. Reported on their -#: own rather than as "unexpected character", because the superscripts -#: ² ³ ⁴ *are* notation (powers), so a subscript is a reasonable guess. -_SUBSCRIPTS_START, _SUBSCRIPTS_END = "\u2080", "\u209c" -_SUBSCRIPT_PLAIN = {chr(0x2080 + d): str(d) for d in range(10)} -_SUBSCRIPT_PLAIN.update({ - "\u2090": "a", "\u2091": "e", "\u2092": "o", "\u2093": "x", - "\u2095": "h", "\u2096": "k", "\u2097": "l", "\u2098": "m", - "\u2099": "n", "\u209a": "p", "\u209b": "s", "\u209c": "t", -}) +__all__ = ["NotationError", "Parser", "tokens"] #: words that join predicates; never names, never implicit factors _KEYWORDS = {"and", "or", "in"} diff --git a/tests/test_core_without_frontend.py b/tests/test_core_without_frontend.py index 3d71a30..301942f 100644 --- a/tests/test_core_without_frontend.py +++ b/tests/test_core_without_frontend.py @@ -24,7 +24,7 @@ #: Importable with pyyaml and nothing else. Proteus imports from every one. CORE = ["ymf", "ymf.archive", "ymf.xdmf", "ymf.cli", "ymf.closure", - "ymf.domain", "ymf.references"] + "ymf.domain", "ymf.references", "ymf.names", "ymf.explicit"] #: Front end, and honestly so. Asserted as well as CORE: if the blocker below #: silently stopped blocking, every CORE assertion would pass for the wrong diff --git a/tests/test_explicit.py b/tests/test_explicit.py new file mode 100644 index 0000000..6ece518 --- /dev/null +++ b/tests/test_explicit.py @@ -0,0 +1,200 @@ +"""A closed spec writes out what the input leaves to rules (C42 round 2). + +ymf.explicit: the coefficients of the state named with the unknowns they +follow (depends_on), and every sub-problem of a split discretization with +its unknown (for), solver, regularization and lags in full (carried_by +for a coefficient read through its carrier). Like the boundary-condition +expansion (#9), it is idempotent, leaves a spec without shorthand alone, +and keys a spec and the archive it closes into alike. +""" + +import copy +import glob +from pathlib import Path + +import pytest + +pytest.importorskip("strictyaml") + +from ymf import closure # noqa: E402 +from ymf.compose import load_composed # noqa: E402 +from ymf.explicit import explicit, state_coefficients, without_derived # noqa: E402 +from ymf.names import FUNCTION_NAMES, names_in # noqa: E402 +from ymf.references import YmfReferenceError, resolve # noqa: E402 + +ROOT = Path(__file__).parent.parent +FIXTURES = Path(__file__).parent / "fixtures" +DAMBREAK = FIXTURES / "two_phase" / "dambreak.ymf" +TORINO = FIXTURES / "shallow_water" / "torino.ymf" + + +def _documents(): + for path in sorted(glob.glob(str(ROOT / "examples" / "**" / "*.ymf"), recursive=True) + + glob.glob(str(FIXTURES / "**" / "*.ymf"), recursive=True)): + try: + yield Path(path), load_composed(path, check_state=False) + except Exception: # an input made to fail, or a fragment + continue + + +def _subproblem(doc, name, discretization="level_set_five"): + entry = next(d for d in doc["solution_paths"]["discretizations"] + if d["name"] == discretization) + return next(s for s in entry["subproblems"] if s["name"] == name) + + +# -- the name reader ----------------------------------------------------------- + +def test_names_in_reads_a_formula_as_the_parser_does(): + declared = {"ρw", "ρa", "χ", "φ", "h", "q", "n", "g", "x", "y", "t"} + assert names_in("ρw (1 − χ) + ρa χ", declared) == {"ρw", "ρa", "χ"} + assert names_in("heaviside(φ)", declared) == {"φ"} + # a function call is not a product of letters (h, n in norm and heaviside) + assert names_in("g n² norm(q) q / h^(7/3)", declared) == {"g", "n", "q", "h"} + # an undeclared run of letters is the product of its letters + assert names_in("gh", declared) == {"g", "h"} + + +def test_the_name_readers_functions_are_the_notations(): + sympy = pytest.importorskip("sympy") # noqa: F841 + from ymf.symbolic.language import Space + from ymf.symbolic.notation import Parser + assert set(Parser(Space(3, {})).functions) == set(FUNCTION_NAMES) + + +def test_state_coefficients_agree_with_the_symbolic_reader_where_it_resolves(): + pytest.importorskip("sympy") + from ymf.symbolic.problem import state_coefficients as symbolic + checked = 0 + for path, doc in _documents(): + strong = doc["model"]["strong_form"] + unknowns = [u if isinstance(u, str) else u["name"] for u in strong["unknowns"]] + try: + resolved = symbolic(strong) + except Exception: + continue + text = state_coefficients(doc) + # the symbolic reader leaves out what does not resolve (an abstract + # model's ρ without ρw); where it resolves, the two agree + for name, components in resolved.items(): + assert text[name] == sorted({c.rsplit("_", 1)[0] if c not in unknowns else c + for c in components}, key=unknowns.index), path + checked += 1 + assert checked >= 9 + + +def test_an_abstract_model_names_its_state_coefficients_without_values(): + doc = load_composed(FIXTURES / "two_phase" / "two_phase_level_set.ymf") + assert state_coefficients(doc) == {"χ": ["φ"], "ρ": ["φ"], "μ": ["φ"]} + + +# -- explicit() ------------------------------------------------------------------ + +def test_explicit_is_idempotent_and_leaves_plain_specs_alone(): + plain, written = 0, 0 + for path, doc in _documents(): + once = explicit(doc) + assert explicit(once) == once, path + if not state_coefficients(doc) and not any( + d.get("subproblems") for d in + (doc.get("solution_paths") or {}).get("discretizations") or []): + from ymf.conditions import expand + assert once == expand(doc), path + plain += 1 + else: + written += 1 + assert plain >= 8 and written >= 4 + + +def test_the_closed_dam_break_names_its_coefficients_of_the_state(): + closed = explicit(load_composed(DAMBREAK)) + coefficients = closed["model"]["strong_form"]["coefficients"] + assert coefficients["χ"] == {"formula": "heaviside(φ)", "units": "1", "depends_on": ["φ"]} + assert coefficients["ρ"]["depends_on"] == ["φ"] + assert coefficients["μ"]["depends_on"] == ["φ"] + assert "depends_on" not in coefficients["ρw"] + + +def test_the_closed_dam_break_writes_its_split_out(): + spec = load_composed(DAMBREAK) + closed = explicit(spec) + flow = _subproblem(closed, "flow") + assert flow["lags"] == {"φ": {"at": "previous_step"}, + "χ": {"via": "vof", "at": "previous_step", "carried_by": "χh"}} + # the entry's solver, refined by the sub-problem's own + assert flow["solver"] == {"type": "nonlinear", "nonlinear_solver": "newton", + "linear_solver": "petsc", "max_iterations": 100} + assert flow["regularize"] == {"heaviside": {"by": "smoothed_heaviside", "width": "ε"}} + vof = _subproblem(closed, "vof") + assert vof["for"] == "χh" + assert vof["solver"]["max_iterations"] == 50 + # redistancing's own sign regularization joins the entry's heaviside + rdls = _subproblem(closed, "rdls") + assert set(rdls["regularize"]) == {"heaviside", "sign"} + assert rdls["regularize"]["sign"]["width"] == "εr" + assert rdls["lags"]["φs"] == {"of": "φ", "via": "ncls", "at": "current_step"} + assert _subproblem(closed, "mcorr")["for"] == "φc" + # the input is not changed + assert _subproblem(spec, "flow")["lags"]["φ"] == "previous_step" + + +def test_the_closed_shallow_water_case_names_friction_and_drag(): + closed = explicit(load_composed(TORINO)) + coefficients = closed["model"]["strong_form"]["coefficients"] + for name in ("τb", "Fd", "Fr", "β"): + assert coefficients[name]["depends_on"] == ["h", "q"], name + assert "depends_on" not in coefficients["θ"] + swe = _subproblem(closed, "swe", "convex_limited") + assert swe["solver"] == {"type": "time_marching", "linear_solver": "petsc"} + + +@pytest.mark.parametrize("path", [DAMBREAK, TORINO, FIXTURES / "shallow_water" / "lake_at_rest.ymf", + ROOT / "examples" / "nonlinear_diffusion.ymf"]) +def test_a_closed_spec_reloads_with_its_specs_keys(tmp_path, path): + spec, outputs, _ = closure.load(path, check_state=False) + closure.write(tmp_path / "closed.ymf", spec, outputs) + again, _, _ = closure.load(tmp_path / "closed.ymf", check_state=False) + assert [k for k, _, _ in closure.planned_outputs(again)] == \ + [k for k, _, _ in closure.planned_outputs(spec)] + assert explicit(again)["solution_paths"] == explicit(spec)["solution_paths"] + assert explicit(again)["model"]["strong_form"]["coefficients"] == \ + explicit(spec)["model"]["strong_form"]["coefficients"] + + +def test_depends_on_does_not_count_towards_a_key(): + # nonlinear_diffusion's keys before C42 round 2 (cekees/ymf main b77262a) + spec, _, _ = closure.load(ROOT / "examples" / "nonlinear_diffusion.ymf") + keys = [k for k, _, _ in closure.planned_outputs(spec)] + assert keys[:3] == ["P1/cells=8/55bb9712a3fc", "P1/cells=16/2b0926f8753a", + "P1/cells=32/f7ada584d9ca"] + assert without_derived(explicit(spec))["model"] == spec["model"] + + +def test_split_problem_reads_the_closed_spec_as_the_spec(): + pytest.importorskip("sympy") + from ymf.symbolic import split_problem + spec = load_composed(DAMBREAK) + assert split_problem(explicit(spec), "level_set_five") == \ + split_problem(spec, "level_set_five") + + +# -- what load refuses ------------------------------------------------------------ + +def test_a_false_depends_on_is_refused(): + doc = copy.deepcopy(load_composed(DAMBREAK)) + doc["model"]["strong_form"]["coefficients"]["ρ"]["depends_on"] = ["v"] + with pytest.raises(YmfReferenceError, match="ρ.depends_on says"): + resolve(doc) + doc["model"]["strong_form"]["coefficients"]["ρ"]["depends_on"] = ["φ"] + resolve(doc) + + +def test_a_false_carried_by_is_refused(): + doc = copy.deepcopy(load_composed(DAMBREAK)) + _subproblem(doc, "flow")["lags"]["χ"]["carried_by"] = "φc" + with pytest.raises(YmfReferenceError, match="χh carries it"): + resolve(doc) + _subproblem(doc, "flow")["lags"]["φ"] = {"at": "previous_step", "carried_by": "χh"} + _subproblem(doc, "flow")["lags"]["χ"]["carried_by"] = "χh" + with pytest.raises(YmfReferenceError, match="no sub-problem introduces"): + resolve(doc) From d0e38ed0ea5e21d0c4666ade576f75b1b701ff46 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 09:31:30 +0100 Subject: [PATCH 2/5] The plain data keeps a coefficient of the state by name MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit § 7 of the C42 model-family note left this for the next PR: adr_problem and split_problem substitute ρ = ρw (1 − χ) + ρa χ into every slot, which is all a runner needs, but a kernel generator (C19) and the dispatcher want ρ once, as an inline function of the state, and the slots over the name. Beside the substituted form, each coefficient of the state now has - `inline`: its formula over the state, constants substituted and other coefficients of the state kept by name (`uses`); and the problem (or sub-problem) has `named`: the `order` to evaluate those coefficients in, each naming only earlier ones, and the ADR slot values (m, f, a, r, H) over the names. The equations are read a second time with each coefficient of the state a known field rather than a constant, so a divergence of ρ v stays a flux and a gradient of ρ is kept (grad_ρ_0). The derivatives and depends_on stay the substituted form's, which are exact; in `named` a name is held fixed. Both are absent without coefficients of the state. The dam break's flow sub-problem now reads: order χ, ρ, μ; χ := χh (VOF's carrier), ρ := 998.2 - 996.995*χ, μ := 0.0010021928 - 0.0009841178*χ; momentum m = v_0*ρ, f = [p + v_0**2*ρ, v_0*v_1*ρ], a = [2*μ, μ], the vertical's r = 9.8*ρ. Torino's drag law keeps Fr, β and Fd (a vector) by name; nonlinear_diffusion's a = [k, k] with k := u**2 + 1. Measured: evaluating each named slot after its inline coefficients equals the substituted slot at 7 random points per name to 1e-13 (dam-break flow 15 slot values, Torino swe, nonlinear_diffusion; tests/test_named_coefficients.py). ymf suite 850 passed, 2 skipped (was 845; 5 new, one test updated for the new key), macOS arm64, petsc-dev env. Co-Authored-By: Claude Opus 5.5 --- docs/changes.md | 9 +++ docs/problem-specification.md | 17 +++- src/ymf/symbolic/problem.py | 108 +++++++++++++++++++++++++ src/ymf/symbolic/split.py | 11 ++- tests/test_named_coefficients.py | 131 +++++++++++++++++++++++++++++++ tests/test_state_coefficients.py | 6 +- 6 files changed, 277 insertions(+), 5 deletions(-) create mode 100644 tests/test_named_coefficients.py diff --git a/docs/changes.md b/docs/changes.md index cb03c0a..8effae1 100644 --- a/docs/changes.md +++ b/docs/changes.md @@ -7,6 +7,15 @@ rediscover why a thing is the way it is. Open gaps are in ## 2026-11 +- **2026-11-05: the plain data keeps a coefficient of the state by name + (C42 round 2).** § 7 of the model-family note: substitution is enough for + a runner, but the generator (C19) and the dispatcher want `ρ` as one + inline function of the state. `adr_problem` and `split_problem` now give + each coefficient of the state its `inline` code and `uses`, and `named` + (`order` and the slot values over the names) beside the substituted + `adr`. Additive, and absent without coefficients of the state; an archive + Proteus wrote for such a spec before this reports "the ADR form differs" + under `--check`, because `ymf_run` records the whole `adr_problem`. - **2026-11-05: closure writes out the split and names the coefficients of the state (C42 round 2).** The same reasoning as the alias expansion: a closed dam break said `φ: previous_step` and left the sub-problems' diff --git a/docs/problem-specification.md b/docs/problem-specification.md index 3e63fc1..45f7d5a 100644 --- a/docs/problem-specification.md +++ b/docs/problem-specification.md @@ -410,8 +410,17 @@ with an unknown. 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. +the state rather than its expansion in every term. Each also has its +`inline` code -- its formula over the state, constants substituted, other +coefficients of the state kept by name (`uses`) -- and the problem gains +`named`: the `order` to evaluate those coefficients in (each names only +earlier ones) and the ADR form's slot values over the names (`m: v_0*ρ`, +`a: [2*μ, μ]`). The equations are read a second time for it, with each +coefficient of the state a known field, so a gradient of one is kept (as +`grad_ρ_0`) rather than lost as a constant's. The derivatives and +`depends_on` are only in the substituted form, which is exact; in `named` a +name is held fixed. A problem with no such coefficient has neither key and +reads exactly as before. A closed spec (an archive) names them too: closure writes each coefficient of the state in the map form with `depends_on`, the unknowns it follows, @@ -819,7 +828,9 @@ says so), which the plain data names (`v_0`) and whose gradient it names 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 +carries what it reads and when, its coefficients of the state by name +(with `inline` code and its ADR form over them in `named`, as for +`adr_problem`), 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 diff --git a/src/ymf/symbolic/problem.py b/src/ymf/symbolic/problem.py index 1f32fae..144c011 100644 --- a/src/ymf/symbolic/problem.py +++ b/src/ymf/symbolic/problem.py @@ -721,6 +721,17 @@ def adr_problem(doc: Dict[str, Any], hamiltonian_gradients=()) -> Dict[str, Any] # slot. Absent when there are none, so a problem without them reads # exactly as before. out["coefficients"] = named + # and the problem over those names: rho once, m = rho*v_0 in the + # mass, for the generator and the dispatcher (C42 round 2) + from ymf.normalize import equations_in_order + state = [c for c in (strong.get("coefficients") or {}) if c in named] + form, inline = _named_form(equations_in_order(strong), + strong.get("coefficients") or {}, coefficients, state, + space.dim, ranks, + hamiltonian_gradients=hamiltonian_gradients) + for name, record in inline.items(): + named[name].update(record) + out["named"] = form return out @@ -767,6 +778,103 @@ def _state_coefficient_data(coefficients, space: Space) -> Dict[str, Any]: return out +#: the keys of an ADR slot that hold its value (its derivatives and flags +#: are the substituted form's, see _named_form) +_VALUE_KEYS = {"mass": ("m",), "advection": ("f",), "reaction": ("r",), + "hamiltonian": ("H",)} + + +def _values_only(adr: Dict[str, Any]) -> Dict[str, Any]: + """An ADR form with each slot's value and nothing else.""" + out = [] + for entry in adr["equations"]: + kept: Dict[str, Any] = {"component": entry["component"]} + for slot, keys in _VALUE_KEYS.items(): + if slot in entry: + kept[slot] = {k: entry[slot][k] for k in keys} + if "diffusion" in entry: + kept["diffusion"] = {phi: {k: d[k] for k in ("rowptr", "colind", "a")} + for phi, d in entry["diffusion"].items()} + out.append(kept) + return {"dim": adr["dim"], "components": adr["components"], "equations": out} + + +def _named_form(texts: Sequence[str], spec: Dict[str, Any], values: Dict[str, Any], + state: Sequence[str], dim: int, unknowns: Dict[str, int], + fields: Optional[Dict[str, int]] = None, + parameters: Sequence[sympy.Symbol] = (), + functions: Optional[Dict[str, Any]] = None, + preset: Optional[Dict[str, Any]] = None, + hamiltonian_gradients=()) -> Tuple[Dict[str, Any], Dict[str, Any]]: + """The problem with each coefficient of the state kept by name. + + The ADR form substitutes ``ρ = ρw (1 − χ) + ρa χ`` into every slot; a + kernel generator wants ``m = ρ v`` and ``ρ`` once, as an inline function + of the state. So the equations are read a second time with each + coefficient of the state a known field (so ``∇ρ`` is not lost, as it + would be for a constant), and each coefficient's own formula with the + others likewise. Returns ``(named, inline)``: + + * ``named``: ``{"order": [...], "equations": [...]}``, the coefficients + of the state in an order in which each names only earlier ones, and + the ADR form's slot values (``m``, ``f``, ``a``, ``r``, ``H``) over + those names. The derivatives and ``depends_on`` are the substituted + form's (``adr``), which are exact: here a name is held fixed. + * ``inline``: ``{coefficient: {"inline": code, "uses": [...]}}``, its + formula over the unknowns, the fields it reads and the other + coefficients of the state it ``uses``; constants are substituted. + """ + from sympy.core.function import AppliedUndef + from ymf.symbolic.language import parse_expression + preset = preset or {} + state = list(state) + ranks = {c: (1 if isinstance(values[c], sympy.MatrixBase) else 0) for c in state} + constants = {k: v for k, v in values.items() if k not in ranks} + named_fields = dict(fields or {}) + named_fields.update(ranks) + + def stem(name): + base, _, index = name.rpartition("_") + return base if base in ranks and index.isdigit() else name + + inline: Dict[str, Any] = {} + uses: Dict[str, List[str]] = {} + for name in state: + if name in preset: + value = preset[name] + else: + entry = spec[name] + text = entry.get("formula") if isinstance(entry, dict) else str(entry) + others = {k: r for k, r in named_fields.items() if k != name} + value = parse_expression(text, Space(dim, unknowns, constants, others, + parameters, functions)) + entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] + entries = [sympy.sympify(e) for e in entries] + used = sorted({stem(f.func.__name__) for e in entries for f in e.atoms(AppliedUndef)} + & set(ranks), key=state.index) + uses[name] = used + inline[name] = {"inline": ([to_code(e) for e in entries] + if isinstance(value, sympy.MatrixBase) else to_code(entries[0]))} + if used: + inline[name]["uses"] = used + order: List[str] = [] + todo = list(state) + while todo: + ready = [c for c in todo if set(uses[c]) <= set(order)] + if not ready: + raise SymbolicError("coefficients %s name each other in a cycle" + % ", ".join(todo)) + order.extend(ready) + todo = [c for c in todo if c not in ready] + + space = Space(dim, unknowns, constants, named_fields, parameters, functions) + residuals: List[sympy.Expr] = [] + for text in texts: + residuals.extend(parse_equation(text, space)) + adr = adr_form(residuals, space, hamiltonian_gradients=hamiltonian_gradients) + return dict({"order": order}, **_values_only(adr)), inline + + def _up_to_a_constant(residuals, space, ranks, dirichlet) -> List[str]: from ymf.symbolic.language import D fixed = {e["component"] for e in dirichlet} diff --git a/src/ymf/symbolic/split.py b/src/ymf/symbolic/split.py index 417e2b5..3aa1cd1 100644 --- a/src/ymf/symbolic/split.py +++ b/src/ymf/symbolic/split.py @@ -83,7 +83,8 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: raise SymbolicError(str(exc)) from exc from ymf.discretizations import subproblem_solver from ymf.normalize import normalize_equations - from ymf.symbolic.problem import _coefficients, _frame, _state_coefficient_data + from ymf.symbolic.problem import (_coefficients, _frame, _named_form, + _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: @@ -212,7 +213,15 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: named = _state_coefficient_data( {k: v for k, v in values.items() if k in spec and k in used}, space) if named: + # by name too: each one inline over the state, and the + # sub-problem's ADR form over those names (problem._named_form) + state = [c for c in spec if c in named] + form, inline = _named_form([t for _, t in texts], spec, values, state, dim, + own, fields, supplied, functions, preset) + for name, extra in inline.items(): + named[name].update(extra) record["coefficients"] = named + record["named"] = form if sub.get("initial_condition"): record["initial"] = _initial(sub, own, ranks, values, functions, dim, where) if sharp: diff --git a/tests/test_named_coefficients.py b/tests/test_named_coefficients.py new file mode 100644 index 0000000..c5c4b1a --- /dev/null +++ b/tests/test_named_coefficients.py @@ -0,0 +1,131 @@ +"""The plain data keeps a coefficient of the state by name (C42 round 2). + +adr_problem and split_problem substitute ρ = ρw (1 − χ) + ρa χ into every +slot, which is all a solver of one system needs. A kernel generator (C19) +wants ρ once, as an inline function of the state, and the slots over the +name. So beside the substituted form they give, for each coefficient of the +state, its `inline` code over the state and the other names it `uses`, and +`named`: the order to evaluate them in and the ADR slot values over them. +Evaluating the named slots after the inline coefficients gives the +substituted slots, which these tests check at random points. +""" + +from pathlib import Path + +import pytest + +numpy = pytest.importorskip("numpy") +pytest.importorskip("sympy") +pytest.importorskip("strictyaml") + +from ymf.compose import load_composed # noqa: E402 +from ymf.symbolic import adr_problem, split_problem # noqa: E402 + +ROOT = Path(__file__).parent.parent +FIXTURES = Path(__file__).parent / "fixtures" +DAMBREAK = FIXTURES / "two_phase" / "dambreak.ymf" +TORINO = FIXTURES / "shallow_water" / "torino.ymf" + + +class _Random(dict): + """Names a code string reads, each a random positive array on first use.""" + + def __init__(self, rng): + super().__init__() + self.rng = rng + + def __missing__(self, name): + if name == "numpy" or name.startswith("__"): # a global, or pytest's own + raise KeyError(name) + self[name] = self.rng.uniform(0.5, 1.5, 7) + return self[name] + + +def _evaluate(code, names): + return eval(code, {"numpy": numpy}, names) # noqa: S307 -- our own generated code + + +def _slots(equation): + """(where, code) for every slot value of one equation of an ADR form.""" + for slot, key in (("mass", "m"), ("reaction", "r"), ("hamiltonian", "H")): + if slot in equation: + yield (equation["component"], slot), equation[slot][key] + if "advection" in equation: + for j, f in enumerate(equation["advection"]["f"]): + yield (equation["component"], "advection", j), f + for phi, d in (equation.get("diffusion") or {}).items(): + for j, a in enumerate(d["a"]): + yield (equation["component"], "diffusion", phi, j), a + + +def _check_named_equals_substituted(problem): + names = _Random(numpy.random.default_rng(42)) + for coefficient in problem["named"]["order"]: + inline = problem["coefficients"][coefficient]["inline"] + if isinstance(inline, list): + for i, code in enumerate(inline): + names["%s_%d" % (coefficient, i)] = _evaluate(code, names) + else: + names[coefficient] = _evaluate(inline, names) + checked = 0 + for substituted, named in zip(problem["adr"]["equations"], problem["named"]["equations"]): + assert dict(_slots(substituted)).keys() == dict(_slots(named)).keys() + for (where, code), (_, named_code) in zip(_slots(substituted), _slots(named)): + a = numpy.broadcast_to(_evaluate(code, names), (7,)) + b = numpy.broadcast_to(_evaluate(named_code, names), (7,)) + assert numpy.allclose(a, b, rtol=1e-13, atol=1e-13), (where, code, named_code) + checked += 1 + return checked + + +def _sub(problem, name): + return next(s for s in problem["subproblems"] if s["name"] == name) + + +def test_the_dam_break_flow_keeps_density_and_viscosity_by_name(): + flow = _sub(split_problem(load_composed(DAMBREAK), "level_set_five"), "flow") + assert flow["named"]["order"] == ["χ", "ρ", "μ"] + coefficients = flow["coefficients"] + assert coefficients["χ"]["inline"] == "χh" # read through VOF's carrier + assert coefficients["ρ"]["inline"] == "998.2 - 996.995*χ" + assert coefficients["ρ"]["uses"] == ["χ"] + assert "uses" not in coefficients["χ"] + momentum = flow["named"]["equations"][0] + assert momentum["component"] == "v_0" + assert momentum["mass"] == {"m": "v_0*ρ"} + assert momentum["advection"] == {"f": ["p + v_0**2*ρ", "v_0*v_1*ρ"]} + assert momentum["diffusion"]["v_0"]["a"] == ["2*μ", "μ"] + assert flow["named"]["equations"][1]["reaction"] == {"r": "9.8*ρ"} + # the derivatives stay the substituted form's, which are exact + assert "dm" not in momentum["mass"] and "depends_on" not in momentum["mass"] + assert _check_named_equals_substituted(flow) == 15 # 6 + 7 + 2 slot values + + +def test_sub_problems_without_coefficients_of_the_state_have_no_named_form(): + problem = split_problem(load_composed(DAMBREAK), "level_set_five") + for name in ("vof", "ncls", "rdls", "mcorr"): + assert "named" not in _sub(problem, name), name + + +def test_torino_keeps_its_drag_law_by_name(): + swe = _sub(split_problem(load_composed(TORINO), "convex_limited"), "swe") + order = swe["named"]["order"] + assert set(order) == set(swe["coefficients"]) + # each names only earlier ones + for i, name in enumerate(order): + assert set(swe["coefficients"][name].get("uses", [])) <= set(order[:i]) + assert isinstance(swe["coefficients"]["Fd"]["inline"], list) # a vector + assert _check_named_equals_substituted(swe) > 0 + + +def test_adr_problem_names_its_coefficients_of_the_state(): + problem = adr_problem(load_composed(ROOT / "examples" / "nonlinear_diffusion.ymf")) + assert problem["named"]["order"] == ["k"] + assert problem["coefficients"]["k"]["inline"] == "u**2 + 1" + assert problem["named"]["equations"][0]["diffusion"]["u"]["a"] == ["k", "k"] + assert _check_named_equals_substituted(problem) > 0 + + +def test_a_problem_without_coefficients_of_the_state_reads_as_before(): + problem = adr_problem(load_composed(ROOT / "examples" / "heat_equation.ymf")) + assert "named" not in problem and "coefficients" not in problem diff --git a/tests/test_state_coefficients.py b/tests/test_state_coefficients.py index 8fb02cb..78f8b83 100644 --- a/tests/test_state_coefficients.py +++ b/tests/test_state_coefficients.py @@ -169,7 +169,8 @@ def test_a_model_declares_the_coefficient_and_a_case_makes_it_state_dependent(tm 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"}}} + "k": {"code": "u**2 + 1", "depends_on": {"u": "nonlinear"}, "d": {"u": "2*u"}, + "inline": "u**2 + 1"}} def test_a_problem_without_one_reads_as_before(): @@ -181,6 +182,9 @@ def test_every_link_of_a_chain_is_named(): κ0={"value": 1e-3}))["coefficients"] assert set(named) == {"κ", "c"} assert named["c"]["d"] == {"T": "2*T"} + # by name: κ over c, c over T, in that order + assert named["κ"]["inline"] == "0.001*c" and named["κ"]["uses"] == ["c"] + assert named["c"]["inline"] == "T**2 + 1" def test_a_vector_coefficient_of_the_state(): From 60a25337c9b07c306ab169003e5d3a802922df54 Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 09:37:14 +0100 Subject: [PATCH 3/5] Regularization: the wet/dry cutoff of a reciprocal (SW2DCV's one_over_hReg) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The model-family note's regularization list (§ 3.1 item 5) had three entries: heaviside and sign, done in round 1, and u = q/h at h = 0, which the shallow-water fixture left to a comment. A split discretization (or one sub-problem) may now give regularize: inverse: {by: wet_dry_cutoff, of: h, width: hε} and every negative power of the scalar unknown h it solves for (q⊗q/h, Manning's h^(-7/3), Fr's h^(-3/2)) is read through Proteus SW2DCV.h's one_over_hReg, 2h/(h² + max(h, hε)²): 1/h wherever h ≥ hε, bounded below, 0 at h = 0. The rewrite happens after parsing and before differentiation, so the derivatives are the regularized expression's, and it reaches the named form's inline coefficients too. `inverse` is the one regularization that is not a function of the notation, so it names its unknown (`of`); checked on load (`of` required for inverse, an unknown, and refused on heaviside or sign). The shallow-water numerics bind hε (meaning: dry_depth, the solver's 1e-5 max(h0)) and use it; split_problem records it under `regularized` with its `of`. Measured: the lake-at-rest flux q_0²/h + 4.905 h² is the unregularized value to 1e-12 at h = 1, 2hε and hε, equals 4.905 h² + 0.09·2h/(h² + hε²) at h = hε/4, and is 0 at h = 0 (hε = 1e-5); Torino's Froude drag still matches the deck to 1e-12 at h = 0.1. ymf suite 854 passed, 2 skipped (was 850; 3 refusal cases, 1 dry-node test; 3 shallow-water tests updated for the regularized code), macOS arm64, petsc-dev env. Co-Authored-By: Claude Opus 5.5 --- docs/changes.md | 7 +++ docs/problem-specification.md | 9 ++- src/ymf/references.py | 33 +++++++---- src/ymf/schema.py | 10 +++- src/ymf/symbolic/problem.py | 12 ++-- src/ymf/symbolic/split.py | 59 ++++++++++++++++--- .../shallow_water/proteus_numerics.ymf | 6 +- .../fixtures/shallow_water/shallow_water.ymf | 3 +- tests/test_shallow_water.py | 32 ++++++++-- tests/test_subproblems.py | 8 ++- 10 files changed, 143 insertions(+), 36 deletions(-) diff --git a/docs/changes.md b/docs/changes.md index 8effae1..2c2601b 100644 --- a/docs/changes.md +++ b/docs/changes.md @@ -7,6 +7,13 @@ rediscover why a thing is the way it is. Open gaps are in ## 2026-11 +- **2026-11-05: the wet/dry cutoff (C42 round 2).** The last + regularization the model-family note lists (§ 3.1 item 5, `u = q/h` at + `h = 0`): `regularize: {inverse: {by: wet_dry_cutoff, of: h, width: hε}}` + reads each negative power of `h` through SW2DCV's `one_over_hReg`. The + shallow-water numerics fixture binds `hε` (`meaning: dry_depth`, the + solver's `1e-5 max(h0)`) and uses it; exact above `hε`, so Torino's drag + and porosity still match the deck to 1e-12. - **2026-11-05: the plain data keeps a coefficient of the state by name (C42 round 2).** § 7 of the model-family note: substitution is enough for a runner, but the generator (C19) and the dispatcher want `ρ` as one diff --git a/docs/problem-specification.md b/docs/problem-specification.md index 45f7d5a..ecad534 100644 --- a/docs/problem-specification.md +++ b/docs/problem-specification.md @@ -787,7 +787,14 @@ previous step left it, or as an earlier sub-problem of this step left it `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`). +functions it replaces with smooth ones (`heaviside`, `sign`), and the +reciprocal of an unknown that may vanish: `inverse: {by: wet_dry_cutoff, +of: h, width: hε}` reads every negative power of the depth `h` (`q⊗q/h`, +Manning's `h^(-7/3)`) through Proteus SW2DCV's `one_over_hReg`, +`2h/(h² + max(h, hε)²)`, which is `1/h` wherever `h ≥ hε` and 0 at `h = 0`. +It is applied before anything is differentiated, so the derivatives are the +regularized expression's; `of` must name a scalar unknown the sub-problem +solves for. 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 diff --git a/src/ymf/references.py b/src/ymf/references.py index 95f85e7..ea2571e 100644 --- a/src/ymf/references.py +++ b/src/ymf/references.py @@ -160,7 +160,24 @@ def _check_equations(doc: Dict[str, Any]) -> None: #: sharp functions of the notation a discretization may regularize -SHARP_FUNCTIONS = ("heaviside", "sign") +SHARP_FUNCTIONS = ("heaviside", "sign", "inverse") + + +def _check_regularize(rules: Dict[str, Any], unknowns, where: str) -> None: + for function, rule in (rules or {}).items(): + if function not in SHARP_FUNCTIONS: + raise YmfReferenceError( + "%s.regularize.%s: only %s are regularized" + % (where, function, ", ".join(SHARP_FUNCTIONS))) + of = rule.get("of") if isinstance(rule, dict) else None + if function == "inverse": + if of is None: + raise YmfReferenceError( + "%s.regularize.inverse: say which unknown's reciprocal (of: h)" % (where,)) + _check(of, sorted(unknowns), "%s.regularize.inverse.of" % (where,), "unknown") + elif of is not None: + raise YmfReferenceError( + "%s.regularize.%s: of is for inverse alone" % (where, function)) def _check_split(doc: Dict[str, Any], entry: Dict[str, Any]) -> None: @@ -179,11 +196,9 @@ def _check_split(doc: Dict[str, Any], entry: Dict[str, Any]) -> None: 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))) + introduced_names = [(s.get("introduces") or {}).get("name") for s in subproblems] + _check_regularize(entry.get("regularize"), + set(unknowns) | {n for n in introduced_names if n}, where) if not subproblems: return names = [s["name"] for s in subproblems] @@ -255,11 +270,7 @@ def _check_split(doc: Dict[str, Any], entry: Dict[str, Any]) -> None: 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))) + _check_regularize(s.get("regularize"), set(unknowns) | set(introduced), here) for field in s.get("writes_back") or {}: _check(field, sorted(set(unknowns) | set(introduced)), "%s.writes_back" % (here,), "model or introduced unknown") diff --git a/src/ymf/schema.py b/src/ymf/schema.py index 17ce381..67e8a0a 100644 --- a/src/ymf/schema.py +++ b/src/ymf/schema.py @@ -148,6 +148,11 @@ Optional("proteus"): Any(), # an interim hint, see SubproblemDef }) +#: how one sharp function is made smooth: heaviside: {by: smoothed_heaviside, +#: width: ε}; inverse: {by: wet_dry_cutoff, of: h, width: hε} reads every +#: negative power of the unknown h through its regularized reciprocal +RegularizeDef = Map({"by": Str(), Optional("width"): Str(), Optional("of"): Str()}) + SubproblemDef = Map({ "name": Str(), # model equations (by name) this sub-problem solves @@ -181,7 +186,7 @@ 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("regularize"): MapPattern(Str(), RegularizeDef), Optional("finite_element"): FiniteElementDef, Optional("time_integration"): Str(), Optional("stabilization"): Any(), @@ -212,8 +217,7 @@ 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("regularize"): MapPattern(Str(), RegularizeDef), Optional("step_control"): Map({ Optional("type"): Str(), Optional("cfl"): Float(), diff --git a/src/ymf/symbolic/problem.py b/src/ymf/symbolic/problem.py index 144c011..06ff280 100644 --- a/src/ymf/symbolic/problem.py +++ b/src/ymf/symbolic/problem.py @@ -805,7 +805,7 @@ def _named_form(texts: Sequence[str], spec: Dict[str, Any], values: Dict[str, An parameters: Sequence[sympy.Symbol] = (), functions: Optional[Dict[str, Any]] = None, preset: Optional[Dict[str, Any]] = None, - hamiltonian_gradients=()) -> Tuple[Dict[str, Any], Dict[str, Any]]: + hamiltonian_gradients=(), post=None) -> Tuple[Dict[str, Any], Dict[str, Any]]: """The problem with each coefficient of the state kept by name. The ADR form substitutes ``ρ = ρw (1 − χ) + ρa χ`` into every slot; a @@ -823,7 +823,11 @@ def _named_form(texts: Sequence[str], spec: Dict[str, Any], values: Dict[str, An * ``inline``: ``{coefficient: {"inline": code, "uses": [...]}}``, its formula over the unknowns, the fields it reads and the other coefficients of the state it ``uses``; constants are substituted. + + ``post`` rewrites each expression after it is read (a split + discretization's regularized reciprocal, ymf.symbolic.split). """ + post = post or (lambda e: e) from sympy.core.function import AppliedUndef from ymf.symbolic.language import parse_expression preset = preset or {} @@ -846,8 +850,8 @@ def stem(name): entry = spec[name] text = entry.get("formula") if isinstance(entry, dict) else str(entry) others = {k: r for k, r in named_fields.items() if k != name} - value = parse_expression(text, Space(dim, unknowns, constants, others, - parameters, functions)) + value = post(parse_expression(text, Space(dim, unknowns, constants, others, + parameters, functions))) entries = list(value) if isinstance(value, sympy.MatrixBase) else [value] entries = [sympy.sympify(e) for e in entries] used = sorted({stem(f.func.__name__) for e in entries for f in e.atoms(AppliedUndef)} @@ -870,7 +874,7 @@ def stem(name): space = Space(dim, unknowns, constants, named_fields, parameters, functions) residuals: List[sympy.Expr] = [] for text in texts: - residuals.extend(parse_equation(text, space)) + residuals.extend(post(r) for r in parse_equation(text, space)) adr = adr_form(residuals, space, hamiltonian_gradients=hamiltonian_gradients) return dict({"order": order}, **_values_only(adr)), inline diff --git a/src/ymf/symbolic/split.py b/src/ymf/symbolic/split.py index 3aa1cd1..e8dcf42 100644 --- a/src/ymf/symbolic/split.py +++ b/src/ymf/symbolic/split.py @@ -51,10 +51,19 @@ def smoothed_sign(phi, eps): return 2 * smoothed_heaviside(phi, eps) - 1 -#: sharp function -> {regularization name: f(argument, width)} +def wet_dry_inverse(h, eps): + """Proteus's SW2DCV ``one_over_hReg``: ``2h / (h² + max(h, ε)²)``, which + is ``1/h`` where ``h ≥ ε`` and goes to 0 with ``h`` instead of blowing up.""" + return 2 * h / (h ** 2 + sympy.Max(h, eps) ** 2) + + +#: sharp function -> {regularization name: f(argument, width)}; ``inverse`` +#: is the reciprocal of the unknown its rule names (``of``), not a function +#: of the notation: every negative power ``h^(-k)`` is read as ``(1/h)_ε^k`` REGULARIZATIONS = { "heaviside": {"smoothed_heaviside": smoothed_heaviside}, "sign": {"smoothed_sign": smoothed_sign}, + "inverse": {"wet_dry_cutoff": wet_dry_inverse}, } @@ -179,7 +188,8 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: 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()} + for fn, (_, width, smoother, _) in smooth.items() if fn != "inverse"} + post = _inverse_cutoff(smooth, own, where) 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, @@ -202,6 +212,12 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: (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) + if post is not None: + regular = [post(r) for r in residuals] + if regular != residuals: + sharp.add("inverse") + residuals = regular + values = {k: (post(v) if k in spec else v) for k, v in values.items()} try: record["adr"] = adr_form(residuals, space) except SymbolicError as exc: @@ -217,7 +233,8 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: # sub-problem's ADR form over those names (problem._named_form) state = [c for c in spec if c in named] form, inline = _named_form([t for _, t in texts], spec, values, state, dim, - own, fields, supplied, functions, preset) + own, fields, supplied, functions, preset, + post=post) for name, extra in inline.items(): named[name].update(extra) record["coefficients"] = named @@ -225,8 +242,9 @@ def split_problem(doc: Dict[str, Any], discretization: str) -> Dict[str, Any]: 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} + record["regularized"] = { + fn: dict({"by": by, "width": to_code(width)}, **({"of": of} if of else {})) + for fn, (by, width, _, of) in smooth.items() if fn in sharp} if sub.get("enforces"): try: parse_equation(sub["enforces"], space) @@ -362,7 +380,8 @@ def _describe_bound(binds, values): def _regularizations(entry, sub, bound, dim, where): - """{sharp function: (by, width, smoother)} for this sub-problem.""" + """{sharp function: (by, width, smoother, of)} for this sub-problem + (``of``: the unknown whose reciprocal ``inverse`` regularizes, else None).""" chosen = dict(entry.get("regularize") or {}) chosen.update(sub.get("regularize") or {}) out = {} @@ -374,10 +393,36 @@ def _regularizations(entry, sub, bound, dim, where): 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"]]) + out[function] = (rule["by"], width, smoothers[rule["by"]], rule.get("of")) return out +def _inverse_cutoff(smooth, own, where): + """The rewrite ``inverse`` asks for, as a function on expressions, or None. + + Every negative power of the scalar unknown it names (``h^(-1)``, + ``h^(-7/3)``) becomes the same power of its regularized reciprocal, before + anything is differentiated, so the derivatives are the regularized + expression's.""" + if "inverse" not in smooth: + return None + _, width, smoother, of = smooth["inverse"] + if of is None or own.get(of) != 0: + raise SymbolicError("%s: regularize.inverse.of must name a scalar unknown this " + "sub-problem solves for (it names %r; it solves for %s)" + % (where, of, ", ".join(own))) + symbol = sympy.Symbol(of, real=True) + regular = smoother(symbol, width) + + def rewrite(expr): + if isinstance(expr, sympy.MatrixBase): + return expr.applyfunc(rewrite) + return sympy.sympify(expr).replace( + lambda p: p.is_Pow and p.base == symbol and p.exp.is_negative is True, + lambda p: regular ** (-p.exp)) + return rewrite + + def _written(texts, smooth) -> set: """The regularized sharp functions these texts use.""" from ymf.symbolic.notation import tokens diff --git a/tests/fixtures/shallow_water/proteus_numerics.ymf b/tests/fixtures/shallow_water/proteus_numerics.ymf index e21c7f7..6283a99 100644 --- a/tests/fixtures/shallow_water/proteus_numerics.ymf +++ b/tests/fixtures/shallow_water/proteus_numerics.ymf @@ -11,8 +11,10 @@ solution_paths: 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} + hε: {meaning: dry_depth} # SW2DCV's hEps = 1e-5 max(h0), set by the solver + regularize: + heaviside: {by: smoothed_heaviside, width: he} # a canopy's edge, over one element + inverse: {by: wet_dry_cutoff, of: h, width: hε} # 1/h -> 2h/(h² + max(h, hε)²) (one_over_hReg) 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 diff --git a/tests/fixtures/shallow_water/shallow_water.ymf b/tests/fixtures/shallow_water/shallow_water.ymf index badbd40..3732992 100644 --- a/tests/fixtures/shallow_water/shallow_water.ymf +++ b/tests/fixtures/shallow_water/shallow_water.ymf @@ -4,7 +4,8 @@ # 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). +# q/h singular, which is the discretization's to regularize (SW2DCV's hEps: +# proteus_numerics.ymf reads 1/h as its one_over_hReg). ymf_state: abstract model: diff --git a/tests/test_shallow_water.py b/tests/test_shallow_water.py index 0f7036f..3a3c02d 100644 --- a/tests/test_shallow_water.py +++ b/tests/test_shallow_water.py @@ -57,8 +57,11 @@ 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"] + # 1/h read through SW2DCV's one_over_hReg (the numerics' inverse cutoff) + assert momentum["advection"]["f"][0] == \ + "4.905*h**2 + 2.0*h*q_0**2/(h**2 + numpy.maximum(h, hε)**2)" + assert "q_0*numpy.sqrt(q_0**2 + q_1**2)/(h**2 + numpy.maximum(h, hε)**2)**(7/3)" in \ + momentum["reaction"]["r"] assert swe["coefficients"]["τb"]["depends_on"] == { "h": "nonlinear", "q_0": "nonlinear", "q_1": "nonlinear"} # the wall set's slip, expanded: q.n = 0 and no tangential traction @@ -87,13 +90,15 @@ 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 set(swe["regularized"]) == {"inverse"} # no heaviside: none is written 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"}} + assert swe["regularized"] == { + "heaviside": {"by": "smoothed_heaviside", "width": "he"}, + "inverse": {"by": "wet_dry_cutoff", "width": "hε", "of": "h"}} 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) @@ -111,9 +116,24 @@ def test_the_torino_drag_is_the_decks_froude_law(torino): 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, + # the cutoff is exact where the depth is above it + assert math.isclose(_evaluate(fd, h=h, q_0=q0, q_1=0.0, x=x, he=0.0725, hε=1e-5), 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 + assert _evaluate(fd, h=h, q_0=q0, q_1=0.0, x=2.0, he=0.0725, hε=1e-5) == 0.0 # outside + + +def test_a_dry_node_has_bounded_fluxes(lake): + """SW2DCV's one_over_hReg: 1/h above hε, 2h/(h² + hε²) below, 0 at h = 0, + so q⊗q/h stays bounded as the depth vanishes.""" + flux = _equation(_swe(lake), "q_0")["advection"]["f"][0] + eps = 1e-5 + for h in (1.0, 2 * eps, eps): + assert math.isclose(_evaluate(flux, h=h, q_0=0.3, hε=eps), 4.905 * h ** 2 + 0.09 / h, + rel_tol=1e-12) + h = eps / 4 + assert math.isclose(_evaluate(flux, h=h, q_0=0.3, hε=eps), + 4.905 * h ** 2 + 0.09 * 2 * h / (h ** 2 + eps ** 2), rel_tol=1e-12) + assert _evaluate(flux, h=0.0, q_0=0.3, hε=eps) == 0.0 def test_both_cases_pass_the_unit_check(lake, torino): diff --git a/tests/test_subproblems.py b/tests/test_subproblems.py index 2bf511b..49a9615 100644 --- a/tests/test_subproblems.py +++ b/tests/test_subproblems.py @@ -119,7 +119,13 @@ def change(doc): (_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"), + "only heaviside, sign, inverse are regularized"), + (_set(SPLIT + ("regularize",), {"inverse": {"by": "wet_dry_cutoff", "width": "h"}}), + "say which unknown's reciprocal"), + (_set(SPLIT + ("regularize",), {"inverse": {"by": "wet_dry_cutoff", "of": "ρ", "width": "h"}}), + "names 'ρ', which is not an unknown"), + (_set(SPLIT + ("regularize",), {"heaviside": {"by": "smoothed_heaviside", "of": "φ"}}), + "of is for inverse alone"), ]) def test_a_split_that_does_not_add_up_is_refused(dambreak, change, message): with pytest.raises(YmfReferenceError, match=message): From 43f3c3b5ecfc508668f784281b53aa3dc179eeed Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 10:11:33 +0100 Subject: [PATCH 4/5] C7.1: a mesh records its fingerprint; a discrete domain loads from a spec file Decided by cekees 2026-11-05 (C7.1 round 2, decision 2): store the mesh's fingerprint in the archive now, in the mesh's provenance, so a restart or `pro gauge` can verify the mesh before reading. ymf carries the definition so every writer computes the same digest; proteus writes and checks it on cekees/proteus#200. - archive.provenance(..., fingerprint=) records {method: ymf-mesh-1, sha256}; the canonicalizer takes the key (anything but {method, sha256} is refused), validation checks the digest's form. - ymf.mesh.fingerprint(nodes, cells): SHA-256 over a tag, the shapes (int64), the nodes (float64, three columns, 2D padded with z = 0) and the cells (int64), little-endian, in the order given. grid_fingerprint(grid) of a stored grid; check_fingerprint(grid, base_dir, digest=None) returns the recorded digest, None if none, raises on another mesh or another method; `digest` is the input-order digest for a writer whose rows are permuted. Two gaps proteus found posing a `discrete` domain from a spec file (proteus#200's end-to-end test worked around both): - the strict spec schema had no `attributes` (or `provenance`) on a discrete domain, though ymf.domain allows them; added. - strictyaml reads the mesh's keys as Any, so inline values came back as text ("'0' is not a number" from the canonicalizer). ymf.domain.type_discrete_block, run on load in compose._check_document, types each data item's values by its data_type (and dimensions, precision, topology counts as integers), refusing text that is not a number. The archive's own refusal of text values (test_inline_values_are_checked) is kept. Measured: tests 764 passed, 11 skipped (was 751 + 13 new: 7 fingerprint, 1 spec-file round trip with attributes and typed values, and the parametrized refusals). The unit square of test_mesh digests to 4631cb7c6e4c040d...; the digest is pinned so a change of definition must change the method. Co-Authored-By: Claude Opus 5.5 --- docs/archive-format.md | 15 +++++++ docs/changes.md | 11 +++++ docs/problem-specification.md | 5 ++- src/ymf/archive.py | 32 +++++++++++++-- src/ymf/compose.py | 3 +- src/ymf/domain.py | 59 +++++++++++++++++++++++++++ src/ymf/mesh.py | 70 ++++++++++++++++++++++++++++++++ src/ymf/schema.py | 3 ++ tests/test_discrete_domains.py | 31 ++++++++++++++ tests/test_mesh_fingerprint.py | 74 ++++++++++++++++++++++++++++++++++ 10 files changed, 296 insertions(+), 7 deletions(-) create mode 100644 tests/test_mesh_fingerprint.py diff --git a/docs/archive-format.md b/docs/archive-format.md index a8192be..7d13d60 100644 --- a/docs/archive-format.md +++ b/docs/archive-format.md @@ -418,6 +418,7 @@ provenance: tool: a mesh converter # free text: the archive names no tools version: '0.1' files: [{name: square.msh, sha256: 3f5a...}] + fingerprint: {method: ymf-mesh-1, sha256: 4631...} # of the mesh itself extra: {...} # anything else the writer keeps ``` @@ -443,6 +444,20 @@ 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. +**The fingerprint** identifies the mesh itself: `ymf.mesh.fingerprint(nodes, +cells)` is the SHA-256 of its coordinates (float64, three columns) and its +cells' node numbers (int64), in the input order, with the shapes and a tag +naming the method (`ymf-mesh-1`, recorded beside the digest so a later +definition can be told apart). A writer that stores the rows in another +order -- a partitioned run writing each rank's rows contiguously -- records +the digest of the input order, so it does not depend on the rank count. A +reader that is about to read data on a mesh (a restart, a gauge) recomputes +it: `ymf.mesh.check_fingerprint(grid, base_dir, digest=None)` returns the +recorded digest, None when there is none, and raises when the grid holds +another mesh; `digest` is the input-order digest the caller computed when the +rows are permuted. Coordinates are hashed exactly: the digest identifies +these numbers, not a mesh up to round-off. + **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 diff --git a/docs/changes.md b/docs/changes.md index f168403..911226e 100644 --- a/docs/changes.md +++ b/docs/changes.md @@ -7,6 +7,17 @@ rediscover why a thing is the way it is. Open gaps are in ## 2026-11 +- **2026-11-05: a mesh records its fingerprint; a discrete domain loads from + a spec file.** cekees, 2026-11-05 (C7.1 round 2): store the mesh's + fingerprint in the archive so a restart or a gauge can check the mesh + before reading. `provenance(..., fingerprint=)` records `{method: + ymf-mesh-1, sha256}`; `ymf.mesh.fingerprint`, `grid_fingerprint` and + `check_fingerprint` define and check it (bytes in `archive-format.md`). + Two gaps proteus found posing a `discrete` domain from a spec file: the + strict schema now takes `attributes` and `provenance` on it (`ymf.domain` + already did), and its inline arrays, which strictyaml reads as text, are + typed by each item's `data_type` on load (`ymf.domain.type_discrete_block`); + the archive itself still refuses text where it expects a number. - **2026-11-04: aliases expand into the abstract forms, and closure records the expansion.** cekees's review note: `slip`, `open`, `inflow`, `characteristic` and condition sets had their meaning only in code diff --git a/docs/problem-specification.md b/docs/problem-specification.md index aeed823..82bbf3a 100644 --- a/docs/problem-specification.md +++ b/docs/problem-specification.md @@ -1158,8 +1158,9 @@ nobody has to rediscover them. Gaps that have been closed, and how, are in [changes.md](changes.md). - **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. + Proteus reads `discrete` domains of triangles or tetrahedra (in memory, + cekees/proteus#200), 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 diff --git a/src/ymf/archive.py b/src/ymf/archive.py index 5623827..5e7b263 100644 --- a/src/ymf/archive.py +++ b/src/ymf/archive.py @@ -667,25 +667,38 @@ def facets( return {"topology": topo, "tags": {"data_item": tags}, "view": view} +#: The method name of :func:`ymf.mesh.fingerprint`'s digest, recorded beside +#: it so that a later definition can be told apart from this one. +FINGERPRINT_METHOD = "ymf-mesh-1" + + def provenance( tool: str, *, version: Optional[str] = None, files: Iterable[Dict[str, str]] = (), extra: Optional[Dict[str, Any]] = None, + fingerprint: Optional[str] = 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. + the files it was made from with their SHA-256 checksums, and the mesh's + own fingerprint. ``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. + ``fingerprint`` is :func:`ymf.mesh.fingerprint` of the mesh in its input + order (a SHA-256 hex digest), recorded as ``{method, sha256}``; a reader + recomputes it to check it holds the mesh it expects before reading data + on it. """ 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 fingerprint is not None: + block["fingerprint"] = {"method": FINGERPRINT_METHOD, "sha256": str(fingerprint).lower()} if extra: block["extra"] = dict(extra) return block @@ -1075,11 +1088,11 @@ def _canonicalize_facets(f: Dict[str, Any], where: str) -> Dict[str, Any]: 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"}) + unknown = sorted(set(p) - {"tool", "version", "files", "fingerprint", "extra"}) if unknown: raise YmfArchiveError( - f"{where}: provenance takes tool, version, files and extra, not " - f"{', '.join(unknown)}; put anything else under extra") + f"{where}: provenance takes tool, version, files, fingerprint and extra, " + f"not {', '.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"]) @@ -1090,6 +1103,12 @@ def _canonicalize_provenance(p: Dict[str, Any], where: str) -> Dict[str, Any]: 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("fingerprint") is not None: + f = p["fingerprint"] + if not isinstance(f, dict) or set(f) != {"method", "sha256"}: + raise YmfArchiveError( + f"{where}.fingerprint: a fingerprint is {{method, sha256}}, got {f!r}") + canonical["fingerprint"] = {"method": str(f["method"]), "sha256": str(f["sha256"]).lower()} if p.get("extra"): if not isinstance(p["extra"], dict): raise YmfArchiveError(f"{where}.extra: must be a mapping") @@ -1355,6 +1374,11 @@ def _validate_grid_like(g: Dict[str, Any], where: str) -> None: raise YmfArchiveError( f"{where}.provenance.files[{i}].sha256: {entry['sha256']!r} is " "not a SHA-256 digest (64 hexadecimal digits)") + f = g["provenance"].get("fingerprint") + if f is not None and not _SHA256.match(f["sha256"]): + raise YmfArchiveError( + f"{where}.provenance.fingerprint.sha256: {f['sha256']!r} is not a " + "SHA-256 digest (64 hexadecimal digits)") def _validate_facets(g: Dict[str, Any], where: str) -> None: diff --git a/src/ymf/compose.py b/src/ymf/compose.py index 288e0fc..5459ea1 100644 --- a/src/ymf/compose.py +++ b/src/ymf/compose.py @@ -66,7 +66,7 @@ import yaml from ymf.archive import is_archive_text -from ymf.domain import validate_domain_block +from ymf.domain import type_discrete_block, validate_domain_block from ymf.references import apply_defaults, resolve from ymf.schema import _ALLOW_FLOW_STYLE, PARTIAL_YMF_SCHEMA, YMF_SCHEMA from ymf.state import check as check_state_of @@ -476,6 +476,7 @@ def _check_document(doc: Dict[str, Any]) -> None: strong = (doc.get("model") or {}).get("strong_form") or {} if "domain" in strong: validate_domain_block(strong["domain"], "model.strong_form.domain") + type_discrete_block(strong["domain"], "model.strong_form.domain") def load_composed(path: str | Path, check_state: bool = True) -> Dict[str, Any]: diff --git a/src/ymf/domain.py b/src/ymf/domain.py index 50a9d5c..54649c8 100644 --- a/src/ymf/domain.py +++ b/src/ymf/domain.py @@ -174,3 +174,62 @@ def box_bounds(domain: Any) -> Optional[List[List[float]]]: if isinstance(domain, dict) and domain.get("kind") == "box": return [list(domain["lower"]), list(domain["upper"])] return None + + +_INTEGER_TYPES = {"Int", "UInt", "Char", "UChar"} + + +def _number(text: Any, integral: bool, where: str) -> Any: + if not isinstance(text, str): + return text + try: + value = float(text) + except ValueError: + raise YmfDomainError("%s: %r is not a number" % (where, text)) from None + if integral: + if value != int(value): + raise YmfDomainError("%s: %r is not an integer" % (where, text)) + return int(value) + return value + + +def type_discrete_block(domain: Any, where: str = "domain") -> Any: + """A ``discrete`` domain's mesh with its numbers as numbers. + + The spec schema reads the mesh's keys as ``Any``, so strictyaml gives + every scalar of an inline mesh back as text; the archive, which the mesh + is written in the vocabulary of, refuses text where it expects numbers. + Each data item's ``values`` are read as its ``data_type`` says, and its + ``dimensions`` and ``precision``, and a topology's counts, as integers. + Anything else is returned as it is. + """ + if not isinstance(domain, dict) or domain.get("kind") != "discrete": + return domain + def item(d: Any, at: str) -> None: + if not isinstance(d, dict): + return + integral = str(d.get("data_type", "Float")) in _INTEGER_TYPES + if isinstance(d.get("values"), list): + flat = d["values"] + d["values"] = [[_number(v, integral, "%s.values" % at) for v in row] if isinstance(row, list) + else _number(row, integral, "%s.values" % at) for row in flat] + if isinstance(d.get("dimensions"), list): + d["dimensions"] = [_number(v, True, "%s.dimensions" % at) for v in d["dimensions"]] + if "precision" in d: + d["precision"] = _number(d["precision"], True, "%s.precision" % at) + def topology(t: Any, at: str) -> None: + if not isinstance(t, dict): + return + for key in ("number_of_elements", "nodes_per_element"): + if key in t: + t[key] = _number(t[key], True, "%s.%s" % (at, key)) + item(t.get("data_item"), at + ".data_item") + topology(domain.get("topology"), where + ".topology") + item((domain.get("geometry") or {}).get("data_item"), where + ".geometry.data_item") + for i, a in enumerate(domain.get("attributes") or []): + item((a or {}).get("data_item"), "%s.attributes[%d].data_item" % (where, i)) + facets = domain.get("facets") + if isinstance(facets, dict): + topology(facets.get("topology"), where + ".facets.topology") + item((facets.get("tags") or {}).get("data_item"), where + ".facets.tags.data_item") + return domain diff --git a/src/ymf/mesh.py b/src/ymf/mesh.py index 8ea51a2..8b768cf 100644 --- a/src/ymf/mesh.py +++ b/src/ymf/mesh.py @@ -20,6 +20,7 @@ from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple from .archive import ( + FINGERPRINT_METHOD, QUADRATIC_EDGES, YmfArchiveError, check_facet_tags, @@ -30,6 +31,9 @@ __all__ = [ "file_sha256", "file_provenance", + "fingerprint", + "grid_fingerprint", + "check_fingerprint", "edge_node_permutation", "reorder_cells", "edge_midpoint_offsets", @@ -89,6 +93,72 @@ def file_provenance( return provenance(tool, version=version, files=files, extra=extra) +def fingerprint(nodes, cells) -> str: + """The SHA-256 digest that identifies a mesh: its geometry and topology, in + the order given (method :data:`ymf.archive.FINGERPRINT_METHOD`). + + The bytes hashed are the tag ``b"ymf mesh fingerprint 1\\0"``; the shapes + ``(n_nodes, 3, n_cells, nodes_per_cell)`` as little-endian int64; the + nodes as little-endian float64 with three columns (2D nodes get z = 0); + the cells' zero-based node numbers as little-endian int64. Coordinates + are hashed exactly, so the digest identifies these numbers, not a mesh up + to round-off. A writer that stores the rows in another order (a + partitioned run writing each rank's rows contiguously) records the + digest of the input order, so it does not depend on the rank count. + """ + import numpy as np + nodes = np.asarray(nodes, dtype=" str: + """:func:`fingerprint` of a uniform grid's mesh, in the order it is stored.""" + import numpy as np + topo, geom = grid["topology"], grid["geometry"] + nodes = np.asarray(load_array(geom["data_item"], base_dir), dtype=float) + dim = {"XY": 2, "X_Y": 2}.get(str(geom.get("type", "XYZ")).upper(), 3) + cells = np.asarray(load_array(topo["data_item"], base_dir)).reshape( + int(topo["number_of_elements"]), -1) + return fingerprint(nodes.reshape(-1, dim), cells) + + +def check_fingerprint(grid: Dict[str, Any], base_dir: str | Path = ".", + digest: Optional[str] = None) -> Optional[str]: + """Check the fingerprint a grid's provenance records. + + ``digest`` is the fingerprint of the mesh in its input order, computed by + the caller when the grid's rows are in another order; by default the + grid's own (:func:`grid_fingerprint`). Returns the recorded digest, or + None when the grid records none; raises when it records another method + or another mesh. + """ + recorded = (grid.get("provenance") or {}).get("fingerprint") + if recorded is None: + return None + if recorded.get("method") != FINGERPRINT_METHOD: + raise YmfArchiveError("provenance.fingerprint: method %r, but this ymf computes %r" + % (recorded.get("method"), FINGERPRINT_METHOD)) + actual = digest if digest is not None else grid_fingerprint(grid, base_dir) + if actual != recorded["sha256"]: + raise YmfArchiveError("provenance.fingerprint: the grid holds mesh %s..., not the " + "recorded %s..." % (actual[:12], recorded["sha256"][:12])) + return recorded["sha256"] + + 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. diff --git a/src/ymf/schema.py b/src/ymf/schema.py index 77e435a..8883ef6 100644 --- a/src/ymf/schema.py +++ b/src/ymf/schema.py @@ -396,8 +396,11 @@ Optional("segment_flags"): Seq(Int()), Optional("facet_flags"): Seq(Int()), # kind: discrete -- already a mesh, in the archive's vocabulary + # (ymf.domain says which of these it needs; facets is above) Optional("topology"): Any(), Optional("geometry"): Any(), + Optional("attributes"): Any(), + Optional("provenance"): Any(), # kind: file -- someone else's language Optional("format"): Str(), Optional("path"): Str(), diff --git a/tests/test_discrete_domains.py b/tests/test_discrete_domains.py index 4023dd4..5a3ac19 100644 --- a/tests/test_discrete_domains.py +++ b/tests/test_discrete_domains.py @@ -100,3 +100,34 @@ def test_a_discrete_domain_is_not_periodic(): "type": "periodic"}) with pytest.raises(SymbolicError, match="a discrete domain cannot be periodic"): adr_problem(doc) + + +# -------------------------------------------------------------------------- +# a discrete domain written in a spec file (C7.1 round 2, two gaps proteus found) +# -------------------------------------------------------------------------- + + +def test_a_mesh_inline_in_a_spec_file_loads_with_its_numbers_and_attributes(tmp_path): + """strictyaml reads the mesh's keys as Any, so its scalars came back as text, which the archive refuses; + and the schema had no attributes (ymf.domain allows them) or provenance.""" + from ymf.archive import validate_approximation + from ymf.domain import type_discrete_block + doc = _on(dict(SQUARE, attributes=[{ + "name": "elementMaterialTypes", "attribute_type": "Scalar", "center": "Cell", + "data_item": _item([1, 2], [2], "Int")}], provenance={"tool": "hand", "files": []}), {"levels": 1}) + import yaml + path = tmp_path / "square.ymf" + path.write_text(yaml.safe_dump(doc, allow_unicode=True, sort_keys=False, default_flow_style=False), + encoding="utf-8") + loaded = load_ymf(path)["model"]["strong_form"]["domain"] + assert loaded["topology"]["data_item"]["values"] == [0, 1, 2, 0, 2, 3] + assert loaded["geometry"]["data_item"]["values"] == [0.0, 0.0, 2.0, 0.0, 2.0, 1.0, 0.0, 1.0] + assert loaded["geometry"]["data_item"]["dimensions"] == [4, 2] + assert loaded["attributes"][0]["data_item"]["values"] == [1, 2] + grid = dict((k, loaded[k]) for k in ("topology", "geometry", "attributes", "facets")) + validate_approximation({"times": [{"time": 0.0, "spaces": {"mesh": grid}}]}) + with pytest.raises(Exception, match="is not a number"): + type_discrete_block({"kind": "discrete", "geometry": {"data_item": {"values": ["x"]}}}) + with pytest.raises(Exception, match="is not an integer"): + type_discrete_block({"kind": "discrete", "topology": {"data_item": { + "data_type": "Int", "values": ["1.5"]}}}) diff --git a/tests/test_mesh_fingerprint.py b/tests/test_mesh_fingerprint.py new file mode 100644 index 0000000..21c9fd1 --- /dev/null +++ b/tests/test_mesh_fingerprint.py @@ -0,0 +1,74 @@ +"""A mesh's fingerprint: the digest a writer records in its provenance, so a reader can check it holds the +mesh it expects before reading data on it (C7.1 step I, decided by cekees 2026-11-05).""" + +import pytest + +from ymf.archive import ( + FINGERPRINT_METHOD, + YmfArchiveError, + provenance, + read_ymf, + validate_approximation, + write_ymf, +) + +np = pytest.importorskip("numpy") + +from ymf.mesh import check_fingerprint, fingerprint, grid_fingerprint # noqa: E402 + +from test_mesh import CELLS, NODES, mesh_domain, square_space # noqa: E402 + +#: the digest of test_mesh's unit square; a change to the definition changes it, and must change the method +SQUARE_DIGEST = "4631cb7c6e4c040d31e317e86467b298eaea5d0c1cae3ea6e1e2bd8def99a91e" + + +def test_the_digest_is_of_the_numbers_in_their_order(): + digest = fingerprint(NODES, CELLS) + assert len(digest) == 64 and digest == fingerprint(NODES[:, :2], CELLS.astype(np.int64)) + moved = NODES.copy() + moved[2, 0] += 1e-15 + assert fingerprint(moved, CELLS) != digest + assert fingerprint(NODES, CELLS[::-1]) != digest # another cell order + assert fingerprint(NODES, CELLS[:, [1, 2, 0]]) != digest # another node order within a cell + + +def test_the_definition_is_pinned(): + assert fingerprint(NODES, CELLS) == SQUARE_DIGEST + + +def test_a_grid_records_its_fingerprint_and_reads_it_back(tmp_path): + digest = fingerprint(NODES, CELLS) + source = provenance("a mesh converter", fingerprint=digest.upper()) + assert source["fingerprint"] == {"method": FINGERPRINT_METHOD, "sha256": digest} + domain = validate_approximation(mesh_domain(square_space(source=source))) + 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 grid_fingerprint(space) == digest + assert check_fingerprint(space) == digest + + +def test_a_grid_holding_another_mesh_is_refused(): + space = square_space(source=provenance("t", fingerprint=fingerprint(NODES, CELLS[::-1]))) + with pytest.raises(YmfArchiveError, match="not the recorded"): + check_fingerprint(space) + # a writer that stored its rows in another order passes the digest of the input order + assert check_fingerprint(space, digest=fingerprint(NODES, CELLS[::-1])) + + +def test_no_fingerprint_checks_nothing_and_another_method_is_refused(): + assert check_fingerprint(square_space(source=provenance("t"))) is None + assert check_fingerprint(square_space()) is None + space = square_space(source=provenance("t", fingerprint="ab" * 32)) + space["provenance"]["fingerprint"]["method"] = "some-other-1" + with pytest.raises(YmfArchiveError, match="method 'some-other-1'"): + check_fingerprint(space) + + +@pytest.mark.parametrize("block, match", [ + ({"tool": "t", "fingerprint": "ab" * 32}, r"a fingerprint is \{method, sha256\}"), + ({"tool": "t", "fingerprint": {"method": FINGERPRINT_METHOD, "sha256": "12"}}, "not a SHA-256 digest"), +]) +def test_bad_fingerprints_are_refused(block, match): + with pytest.raises(YmfArchiveError, match=match): + validate_approximation(mesh_domain(square_space(source=block))) From 7ae230c9308edba29e0113aed63c14feaf06462d Mon Sep 17 00:00:00 2001 From: cekees <613654+cekees@users.noreply.github.com> Date: Thu, 5 Nov 2026 10:17:04 +0100 Subject: [PATCH 5/5] discrete domains: an inline planar mesh stored as XYZ is posed in its `dimension` Found by proteus once the previous commit let a spec carry an inline mesh: proteus writes 2D meshes with three coordinates (XYZ, z = 0), and adr_problem took the space dimension from the geometry's columns, so `dimension: 2` was ignored: dim 3, upper z = 0, and proteus solved a 3D problem on a 2D mesh (L2 error nan for square4x4_tagged). _inline_nodes now keeps the first `dimension` coordinates when the block gives one, refuses coordinates beyond it that vary (a tilted mesh is not planar), and refuses a dimension above the geometry's width. Measured: ymf tests 765 passed, 11 skipped (+1). proteus (cekees/proteus#200, test_a_discrete_domain_inline_with_its_attributes_solves_as_its_file): square4x4_tagged inline with its node and element tags, before L2 nan at levels 1 and 2, after < 1e-12 and the archive bit-equal to the `kind: file` run. Co-Authored-By: Claude Opus 5.5 --- src/ymf/symbolic/problem.py | 17 +++++++++++++++-- tests/test_discrete_domains.py | 14 ++++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/ymf/symbolic/problem.py b/src/ymf/symbolic/problem.py index 1f32fae..5e43866 100644 --- a/src/ymf/symbolic/problem.py +++ b/src/ymf/symbolic/problem.py @@ -439,8 +439,21 @@ def _inline_nodes(region: Dict[str, Any]) -> Optional[List[List[float]]]: 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)] + width = _geometry_dimension(region) or int((item.get("dimensions") or [0, 3])[-1]) + rows = [flat[i:i + width] for i in range(0, len(flat), width)] + dim = int(region.get("dimension") or width) + if dim > width: + raise SymbolicError("strong_form.domain: dimension %d, but its nodes have %d coordinates" + % (dim, width)) + if dim < width: + # a planar mesh stored with three coordinates (XYZ), as proteus writes + # one: its space is the first `dimension` of them, the rest constant + dropped = set(tuple(r[dim:]) for r in rows) + if len(dropped) > 1: + raise SymbolicError("strong_form.domain: dimension %d, but its nodes vary in " + "coordinate %d" % (dim, dim + 1)) + rows = [r[:dim] for r in rows] + return rows def _frame(doc: Dict[str, Any]): diff --git a/tests/test_discrete_domains.py b/tests/test_discrete_domains.py index 5a3ac19..46c842b 100644 --- a/tests/test_discrete_domains.py +++ b/tests/test_discrete_domains.py @@ -131,3 +131,17 @@ def test_a_mesh_inline_in_a_spec_file_loads_with_its_numbers_and_attributes(tmp_ with pytest.raises(Exception, match="is not an integer"): type_discrete_block({"kind": "discrete", "topology": {"data_item": { "data_type": "Int", "values": ["1.5"]}}}) + + +def test_a_planar_mesh_stored_with_three_coordinates_is_two_dimensional(): + """proteus writes 2D meshes as XYZ with z = 0; `dimension: 2` says the space is the plane.""" + flat = copy.deepcopy(SQUARE) + flat["geometry"] = {"type": "XYZ", "data_item": _item([0, 0, 0, 2, 0, 0, 2, 1, 0, 0, 1, 0], [4, 3])} + flat["dimension"] = 2 + problem = adr_problem(_on(flat, {"levels": 1})) + assert problem["adr"]["dim"] == 2 + assert problem["domain"]["lower"] == [0.0, 0.0] and problem["domain"]["upper"] == [2.0, 1.0] + tilted = copy.deepcopy(flat) + tilted["geometry"]["data_item"]["values"][2] = 0.5 + with pytest.raises(SymbolicError, match="vary in coordinate 3"): + adr_problem(_on(tilted, {"levels": 1}))