Source code for femorph_solver.validation._benchmark

"""Benchmark-problem + published-value dataclasses.

A :class:`BenchmarkProblem` couples a mesh-building recipe to a
catalogue of :class:`PublishedValue` references.  Subclasses
implement two methods:

* ``build_model(**mesh_params)`` — construct a
  :class:`~femorph_solver.Model` ready to solve.
* ``extract(model, result, name)`` — pull the quantity that
  the named :class:`PublishedValue` predicts.

Running ``problem.validate(**mesh_params)`` builds, solves,
extracts, and returns a :class:`ValidationResult` per published
quantity.
"""

from __future__ import annotations

import resource
import time
from abc import ABC, abstractmethod
from dataclasses import dataclass, field
from typing import TYPE_CHECKING, Any

if TYPE_CHECKING:
    from femorph_solver.model import Model


def _peak_rss_mb() -> float:
    """Return the current process peak RSS in MiB.

    ``resource.ru_maxrss`` is in KiB on Linux (and bytes on macOS,
    which the caller must normalise separately if that matters);
    we convert from the Linux convention since that's where CI
    runs.
    """
    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024.0


[docs] @dataclass(frozen=True) class PublishedValue: """One published reference value associated with a benchmark. The ``source`` string must uniquely identify the publication: ``"Timoshenko 1955 §5.4"``, ``"Rao 2019 Table 8.1"``, ``"Irons & Razzaque 1972"`` — enough that a reviewer can find the citation in the reference list without further lookup. ``tolerance`` is the *acceptance* tolerance for this quantity at the recommended mesh — looser on coarse meshes, tighter on finer ones. Convergence studies may use a stricter ``tolerance_fn(n_dof)`` later; the scalar is the floor. """ name: str value: float unit: str = "" source: str = "" formula: str = "" tolerance: float = 0.05
[docs] @dataclass(frozen=True) class ValidationResult: """One validation attempt at a specific mesh refinement. ``rel_error`` is the signed relative error (``(computed - published) / published``). ``passed`` is set when the absolute relative error is within the :class:`PublishedValue`'s ``tolerance``. All other fields are the raw inputs so the caller can re-render the comparison without re-solving. ``wall_s`` and ``peak_rss_mb`` capture the wall-clock time and peak resident-set size of the solve that produced this row, so validation outputs can feed the :mod:`femorph_solver.estimators` training loader — the same way benchmark rows do. """ problem_name: str published: PublishedValue computed: float mesh_params: dict[str, Any] = field(default_factory=dict) n_dof: int | None = None wall_s: float | None = None peak_rss_mb: float | None = None @property def abs_error(self) -> float: return self.computed - self.published.value @property def rel_error(self) -> float: if self.published.value == 0.0: return float("inf") if self.computed != 0.0 else 0.0 return (self.computed - self.published.value) / self.published.value @property def passed(self) -> bool: # When the published value is exactly zero, an absolute # check against ``tolerance`` is the appropriate semantic # (reference problems with ``value=0`` — like the Irons # patch test's ``max_strain_error`` — want a machine-precision # floor, not a relative comparison that's always ``inf``). if self.published.value == 0.0: return abs(self.computed) <= self.published.tolerance return abs(self.rel_error) <= self.published.tolerance
[docs] class BenchmarkProblem(ABC): """Declarative canonical FEA problem with published reference(s). Subclasses declare: * :attr:`name` — short identifier, ``snake_case``. * :attr:`description` — one-sentence summary. * :attr:`published_values` — tuple of :class:`PublishedValue`. * :meth:`build_model` — construct a ready-to-solve model. * :meth:`extract` — pull a named quantity from the solve output. The base class provides :meth:`validate` which composes ``build_model`` → ``solve/modal_solve`` → ``extract`` and returns a list of :class:`ValidationResult`. """ name: str description: str published_values: tuple[PublishedValue, ...] #: Analysis type — ``"static"`` or ``"modal"``. Subclasses #: override. analysis: str = "static"
[docs] @abstractmethod def build_model(self, **mesh_params: Any) -> Model: """Construct a :class:`~femorph_solver.Model` ready to solve."""
[docs] @abstractmethod def extract(self, model: Model, result: Any, name: str) -> float: """Pull the quantity named ``name`` from the solve output. ``result`` is whatever ``model.solve()`` / ``model.modal_solve()`` returned, unboxed. ``name`` is the :class:`PublishedValue` ``name``. """
[docs] def validate(self, **mesh_params: Any) -> list[ValidationResult]: """Build, solve, compare against every published value. Returns one :class:`ValidationResult` per :attr:`published_values` entry. """ model = self.build_model(**mesh_params) rss_before = _peak_rss_mb() t0 = time.perf_counter() if self.analysis == "static": result = model.solve() elif self.analysis == "modal": # Problem classes can override ``default_n_modes`` to # request more than the single fundamental — useful when # the extract has to filter for a specific mode shape # (e.g. the first *bending* mode on a cantilever, or the # (1, 1) plate mode when spurious in-plane modes may # appear below it on coarse meshes). default_n = int(getattr(self, "default_n_modes", 1)) n_modes = max(1, int(mesh_params.get("n_modes", default_n))) result = model.modal_solve(n_modes=n_modes) else: msg = f"unknown analysis type {self.analysis!r}" raise ValueError(msg) wall_s = time.perf_counter() - t0 # ``ru_maxrss`` is process-wide high-water; report the # delta against the pre-solve snapshot so subsequent # validate() calls in the same process don't inherit the # previous mesh's memory peak as a spurious baseline. peak_rss_mb = max(_peak_rss_mb() - rss_before, 0.0) n_dof = int(model.n_nodes * 3) # translational-only fallback out = [] for pv in self.published_values: computed = float(self.extract(model, result, pv.name)) out.append( ValidationResult( problem_name=self.name, published=pv, computed=computed, mesh_params=dict(mesh_params), n_dof=n_dof, wall_s=wall_s, peak_rss_mb=peak_rss_mb, ) ) return out