MeanField_SERiF/docs/static/operators.md

Applying equations as discrete operators

An equation describes local physics and its chosen weak formulation. An operator evaluates that description on a discretization to produce a residual R(u) and a Jacobian action J(u) du. New physics using supported operands needs an equation, not another bespoke operator.

The MFEM volume backend supports global scalar VALUE, scalar H1/L2 VALUE, scalar H1 GRADIENT and vector H1 VALUE/DIVERGENCE/GRADIENT operands and tests on three-dimensional conforming serial meshes. Scalar equations, mass and angular-momentum constraints, and primal gravity use the same volume graph and numerical application functions. Surface integrals, H(div), H(curl), and parallel assembly remain outside the current scope.

Where things belong

Responsibility	Implementation
Physical laws, parameters, admissibility, local partial derivatives	equations
Existing field and term identities, spaces, supports, operands	discretization
Bind a spatial term to an MFEM space or a global scalar to its mesh	TermBinding
Supported coefficient indices, orientation, extraction, assembly	DofMap
Traverse an equation's additive contributions in declaration order	operators/contributions.hpp
Interpret operand order, row support, and quadrature forms	EquationPlan
Bind physical terms to coefficient vectors in an ordinary dependency context	CoefficientView, make_term_coefficients
Verify and reuse canonical tensor factors	TensorBasisCache
Contract a dense basis or tensor factors with pointwise maps	BasisEvaluation
Evaluate a supported operation and apply its basis or transpose	OperandEvaluator
Name volume stages and describe their fixed configuration and results	volume_data.hpp
Sample bases and operands, map geometry, and evaluate one element's equation	element_preparation.hpp
Apply these stages over one contribution's supported elements	volume_stages.hpp
Validate and describe current volume sources	volume_inputs.hpp
Declare volume graph nodes and return an ordinary dependency cache	volume_computations.hpp
Consume accepted volume data and apply its Jacobian or transpose	volume_application.hpp
Integrate and differentiate accepted element data	element_evaluation.hpp
Compare exact MFEM coefficient values and persistent binding identities	CoefficientTracking
Track changes, retain accepted results, and publish successful preparations	utils::dependencies
Expose the frozen linearization to MFEM solvers	JacobianOperator

The volume path separates input construction, graph declarations, accepted data, and numerical application. The dependency cache is their only preparation owner. Operators consume the equation interface; concrete physics headers do not depend on operators or MFEM. Geometry and quadrature continue to use the existing mapping and discretization modules.

A binding identifies a term, not a new field catalog and not an operand operation. Gravity's potential and gradient are different terms. Multiple operations on a term share one coefficient vector and one packed Jacobian column. Plan::Operands lists unique sampled operations; Plan::Columns lists unique coefficient terms.

Reading the calculation

make_contribution_computations in volume_computations.hpp reads in numerical order: reference bases, geometry, operand samples, then equation values. Its short lambdas call the ordinary numerical functions in volume_stages.hpp and element_preparation.hpp. The declarations themselves do no numerical work; try_prepare runs the required nodes later.

The contribution index distinguishes repeated equations, including repeated instances of the same type. Each node returns data for the contribution's supported elements; there is no separate C++ graph node for every mesh point.

Result key (Basis, Geometry, and sample keys are in backend::volume_keys)	Declared numerical dependencies
Basis<I>	MeshAndCoefficientBindings: fixed-mesh coordinate sequence and binding identities; equation and integration options are fixed configuration.
Geometry<I>	Basis<I>, displacement coefficients, compactification coefficients.
ReferenceValues<I, Operand>	Basis<I> and the coefficient source identified by Operand::term_type.
PhysicalValues<I, Operand>	That operand's reference values and Geometry<I>; present only when the operand's physical representation depends on geometry.
QuadratureLinearization<I>	Basis, geometry, and each declared operand's reference or physical values. Residuals, derivatives, weights, and contribution scaling are computed together.
Linearization	Mesh and binding metadata and every contribution's QuadratureLinearization; retain ordered batches, assemble the residual, and validate complete integrals.

There are no manual change flags inside these numerical stages. The dependency engine follows the declared inputs to decide which stages need execution. For scalar VALUE,

u(q) = sum_j E(q,j) c(j)
z(q) = W(q) f(u(q), x(q))
R_element = E_test^T z

element_residual performs the last contraction. element_jacobian_action differentiates the same expression. Coefficients, explicit physical position, integration weight, and (for gradients and H1 divergence) both operand and test bases vary. Scalar and vector H1 value bases stay fixed. The evaluator's evaluate method and its transpose contraction work for all supported operations by flattening a vector sample's components into three rows per point.

OperandValue<Operand> describes the sampled shape using both term rank and operation: Value(u) of a vector is std::array<double, 3>, while Divergence(u) is double. PointLinearization<Operands, Result> supplies scalar/vector/matrix partials of the appropriate shape; the default result remains scalar for existing equations.

There are two maps. MFEM first maps from its canonical element to the fixed computational mesh: RT CalcVShape(transformation, ...) and CalcPhysDivShape(...) already include this map. OperandEvaluator then applies SERiF's additional map A = dx/dX, j = det(A): g_x = A g_X / j, div_x(g_x) = div_X(g_X) / j. The same transformation is applied to the test basis. Geometry variations include delta(A g_X / j) and delta(div_X(g_X) / j) at fixed coefficients, on both sides of the contraction. Cached fixed-mesh samples make these variations independent of subsequent changes to the caller's coefficient vectors.

H1 divergence uses a different transformation: div_x(u) = trace(G A^-1), where G[c,k] = partial_X_k(u_c). Its evaluator retains all nine entries of G internally, then produces one scalar per physical quadrature point. A scalar reference divergence alone would lose the information needed under shear or anisotropic dilation. The geometry derivative is trace(G delta(A^-1)), using the existing mapping variation.

Consequently an evaluator's reference sample and physical value can have different shapes. It declares both sizes, reads its own reference samples, and advertises geometry_dependent. Element assembly uses those declarations without branching on a particular finite-element family. Coefficient evaluation and test integration still share the same matrix-vector and transpose contractions.

prepare_volume validates the current coefficients, constructs their dependency context, and asks the ordinary dependency cache to prepare the graph. The final volume_keys::Linearization result owns the assembled residual and each contribution's immutable element data. Ordinary application functions consume this result; there is no second prepared-operator lifecycle.

Residual and Jacobian calls do not evaluate the EOS or silently prepare new values. Packed columns remain unique terms in first-appearance order (Plan::Columns), followed by displacement.

Model ownership and startup

The model owns caches directly. Fixed discretization, current coefficient inputs, and integration choices have separate names. With bindings, coefficient vectors, and the mapping already constructed, one volume row can be set up as follows:

namespace dependencies = serif::utils::dependencies;
namespace backend = serif::operators::backend::mfem;
namespace equations = serif::equations;

const auto closure_equation = equations::BarotropicDensityBalance{equation_of_state};
using Closure = std::remove_cvref_t<decltype(closure_equation)>;

const backend::VolumeDiscretization closure_discretization{
    .residual_binding = density_binding,
    .displacement_binding = displacement_binding,
    .compactification = compactification,
    .mapper = mapper
};
const auto term_coefficients = backend::make_term_coefficients(
    backend::bind(density_binding, density),
    backend::bind(potential_binding, potential),
    backend::bind(surface_level_binding, surface_level),
    backend::bind(angular_velocity_binding, angular_velocity)
);
const auto closure_inputs = backend::make_volume_inputs<Closure>(
    closure_discretization, term_coefficients, displacement
);
auto closure_cache = backend::make_volume_equation_cache(
    closure_equation, closure_discretization, closure_inputs
);
auto& closure = closure_cache;

// Repeat preparation with the current coefficients during nonlinear iterations.
backend::prepare_volume(closure, closure_discretization, term_coefficients, displacement);
const auto& linearization = closure.get<backend::volume_keys::Linearization>();
backend::volume_residual(linearization, density_residual);
backend::apply_jacobian(linearization, direction, displacement_direction, workspace, action);
backend::JacobianOperator jacobian(closure);

The snippets use <type_traits> for std::remove_cvref_t and <utility> for std::move. The equation sandbox supplies the complete manufactured equilibrium. VolumeOptions optionally supplies named quadrature_policy and contraction_mode members as the factory's fourth argument.

Cache construction declares work without evaluating physics. The graph owns the equation, integration options, and reusable preparation scratch. Meshes, spaces, term bindings, compactification, and mapper are borrowed and must outlive its uses. Coefficient contexts own views, not the coefficient vectors. Keep those vectors alive while preparing; accepted numerical results do not borrow their contents. The workspace used for Jacobian application remains caller-owned and is not a tracked dependency.

make_volume_inputs captures current mesh/binding metadata and borrows current coefficients. Rebuild this context when the discretization changes. Calling prepare_volume or try_prepare_volume does that explicitly on every call and validates before requesting the complete linearization. Calling cache.try_prepare directly is also available with a current, validated context; it is the generic engine's interface and permits partial result requests.

Moving a cache after an adapter borrows it invalidates that reference. The cache must outlive its adapter. The adapter reads the cache's latest accepted linearization on every call, so rejection or a request that omits that linearization makes application unavailable. References returned by cache.get must be reacquired after preparation or invalidation.

Two adapters can share one cache, while separate caches retain independent linearizations. Separate caches are not connected or published as one transaction; a controller prepares them explicitly. Shared intermediates can be declared in a common graph when that becomes useful.

Equation contributions

An equation may describe itself as an additive tuple returned by contributions(); the primal Poisson balance does so for its flux and source. Contributions share the same residual term, while leaves retain independent test operations, supports, and quadrature rules. ContributionsOf collects leaf types for concept checks and operand planning. operators::for_each_contribution visits the corresponding leaf objects in declaration order, borrowing the owned objects, including equations that cannot be copied.

Global scalars: mass, angular momentum, C, and Ω

barotropic_constant::mass_normalization stores the unknown Bernoulli constant C and receives the mass constraint. rotation::angular_velocity stores the signed angular speed Ω about a prescribed axis and center. They are explicitly declared GlobalScalarTerms, with one coefficient each and no finite-element space. Other nonspatial descriptors, such as surface-deformation coefficients, are not implicitly admitted as scalars.

using C = serif::discretization::blocks::barotropic_constant::mass_normalization;
using Omega = serif::discretization::blocks::rotation::angular_velocity;
using serif::discretization::backend::mfem::TermBinding;
namespace backend = serif::operators::backend::mfem;

TermBinding<C, Schema> c_binding(mesh);
TermBinding<Omega, Schema> omega_binding(mesh);
mfem::Vector c(1), omega(1);
c = surface_level;
omega = angular_velocity;
const auto term_coefficients = backend::make_term_coefficients(
    backend::bind(density_binding, density),
    backend::bind(potential_binding, potential),
    backend::bind(c_binding, c),
    backend::bind(omega_binding, omega)
);
const backend::VolumeDiscretization mass_discretization{
    .residual_binding = c_binding,
    .displacement_binding = displacement_binding,
    .compactification = compactification,
    .mapper = mapper
};
const auto mass_inputs = backend::make_volume_inputs<serif::equations::MassNormalization>(
    mass_discretization, term_coefficients, displacement
);
auto mass = backend::make_volume_equation_cache(
    serif::equations::MassNormalization{target_mass}, mass_discretization, mass_inputs
);
backend::prepare_volume(mass, mass_discretization, term_coefficients, displacement);

const backend::VolumeDiscretization angular_momentum_discretization{
    .residual_binding = omega_binding,
    .displacement_binding = displacement_binding,
    .compactification = compactification,
    .mapper = mapper
};
const auto angular_momentum_inputs = backend::make_volume_inputs<serif::equations::AngularMomentum>(
    angular_momentum_discretization, term_coefficients, displacement
);
auto angular_momentum = backend::make_volume_equation_cache(
    serif::equations::AngularMomentum{target_momentum, axis, center},
    angular_momentum_discretization, angular_momentum_inputs
);
backend::prepare_volume(angular_momentum, angular_momentum_discretization, term_coefficients, displacement);

The scalar basis is one at every quadrature point. The same basis contractions, geometry derivatives, coefficient validation, frozen snapshots, and packed adapter serve spatial and global terms. A scalar row sums every supported element into its single entry; there is no artificial order-zero FE space or per-element scalar DOF. Global scalar reductions compensate cancellation at quadrature points and across elements, preserving the old column-sum accuracy. Its mesh, support, and lifetime checks still apply. Preparation also rejects overflow in the integrated scalar residual, including its offset, before publishing the candidate. Calling .space() on a global binding is an error; common assembly uses .mesh() and .full_size().

MassNormalization evaluates rho over the stellar domains. Its residual_offset() returns -target_mass, which final volume assembly adds once after integration. This optional hook is restricted to global scalar rows and describes a constant, state-independent offset. Arithmetic scales the offset along with its contribution; the offset has no density or geometry derivative. The mass row depends on density and displacement, not on C directly. C influences the density-balance row.

AngularMomentum receives the Ω row and evaluates rho Omega s^2, where s is distance from the prescribed physical axis through center. Its offset is -target_momentum. The target must be finite and nonnegative; trial angular speeds may be signed or zero. It reuses RigidRotation for s^2 = 2 Psi_unit and its position derivative. The equation supplies dF/drho = Omega s^2, dF/dOmega = rho s^2, and dF/dx = 2 rho Omega grad(Psi_unit); the existing backend differentiates volume.

An equation may also define validation_integrand(point, linearization) and validate_integral(value). During the existing quadrature loop the backend sums that scalar measure with the physical integration weight, then validates its total over the contribution's support before publishing the preparation. Throwing std::domain_error rejects the trial state. Each leaf is validated separately. Angular momentum uses its already computed dF/dOmega to reject a nonfinite or negative total moment of inertia, including at Ω = 0; it does not require every local trial density to be positive.

The density-balance constructor accepts an axis (default z) and center (default zero), not C or Ω. The axis is normalized once. The balance reads Ω and C from the state. It implements d(A(rho) + Phi - Psi - C)/dC = -1 and dPsi/dOmega = 2 Omega Psi_unit. These derivatives are defined at Ω = 0 without dividing by Ω. RigidRotation remains a fixed-vector utility for independent calculations.

Packed density-balance columns are rho, Phi, C, Omega, displacement. Mass columns are density, displacement; angular-momentum columns are density, Omega, displacement. Each global scalar column occupies one entry. Directions must bind the same scalar descriptors, including an explicit zero when a control is held fixed. Reprepare after changing its value, as with spatial coefficients.

The mass port retains the old 1e-12 integral comparison, 5e-7 affine scaling, 3e-8 density finite difference at step 1e-3, 2e-7 geometry finite difference at step 1e-6, 2e-11 independent complete action, and 2e-15 column sum. Serial transpose actions and preparation work reports are available. MPI reductions and coupled root-layout routing remain deferred.

The angular-momentum port retains the old 2e-13 moment comparison, 7e-7 fifth-power dilation law, 3e-13 parallel-axis identity, 2e-11 independent action, 3e-15 column sum, and 2e-15 angular-speed action comparison. Density, geometry, and angular-speed differences retain thresholds 4e-8, 4e-7, and 2e-10 at steps 1e-3, 1e-6, and 1e-6. Tests exercise tilted and displaced axes, both element types and vector orderings, a curved star, frozen snapshots, rejection and recovery, and default quadrature. Serial cache behavior is checked through the shared preparation reports. Distributed/root-layout tests await their corresponding infrastructure. The supported physical part of the old rotating-system test also checks one shared Ω across density-balance and invariant rows (3e-7 and 3e-10 at step 1e-5), its zero columns in mass, the test closure, and gravity, and the different zero-speed limits. Seed generation, surface projection, and root manifests are not part of this volume-row port.

Prescribed surface pressure

Isobaric prescribes a typed PressureValue. It describes the physical condition without selecting enthalpy or any other field. The model's surface compiler selects a formulation from the available thermodynamic fields, solved equation rows, and EOS relations:

1. Use a pressure unknown directly when one is available.
2. Otherwise use a declared inverse EOS relation to resolve a boundary field.
3. Otherwise evaluate a declared forward pressure relation and subtract the target.

An unavailable or ambiguous choice is rejected at compile time. The compiler uses the EOS's declared relations; it does not numerically invert an EOS. A preferred surface field can narrow the choice when the physical model admits several rows. See AutomaticPressureSurface.

For a polytrope using density, the selected condition reduces to

R_i = rho_i - rho_EOS(P_surface)
delta R_i = delta rho_i

Other EOS relations may depend on temperature, composition, or other bound fields. Their values and derivatives are evaluated from the current boundary state. Zero pressure also goes through the selected relation; its carrier value need not be zero. A pressure unknown needs no EOS evaluation at all.

The finite-element backend compiles traces of the required fields at the selected carrier nodes on StellarSurfaceBoundary. These fields may have different spaces and polynomial orders. The scalar carrier must have a supported nodal H1 row; Bernstein/positive coefficients are rejected because they are not point values.

For standalone operator work, the shared boundary graph can be used directly. Here Formulation is the selected formulation, eos_source explicitly tracks the EOS parameters, and coefficients binds the fields required by that formulation:

namespace physics = serif::geometry::surface;
namespace surface = physics::backend::mfem;
namespace dependencies = serif::utils::dependencies;

const auto boundary = surface::make_pressure_boundary<Formulation>(coefficients);
const auto inputs = surface::make_pressure_boundary_inputs(
    boundary, eos_source,
    physics::Isobaric{serif::dimensions::PressureValue{0.03125}}, coefficients
);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(
    surface::make_pressure_boundary_graph<Formulation>(), inputs
);
const auto prepared = cache.try_prepare(inputs);
if (prepared) {
    const auto& linearization = cache.get<physics::keys::Linearization>();
    // First assemble the volume residual and Jacobian action into these vectors.
    surface::apply_residual_rows(*boundary, linearization, residual);
    surface::apply_jacobian_rows(*boundary, linearization, direction, action);
}
// Handle prepared.error() before using a result from a rejected trial.

The graph tracks pressure, EOS parameters, and field coefficients. Preparation owns the evaluated residuals and derivatives until the next preparation, so later input changes do not silently change an already prepared action. Invalid numerical states return a SurfaceFailure; invalid or stale bindings remain programming errors. Both coefficients and direction must use the bindings that compiled the traces, and those bindings and borrowed values must remain alive.

Apply boundary replacement after volume assembly. It removes every volume coupling at a selected row and inserts every derivative required by the chosen surface relation. For the transpose, mask selected residual-dual entries before applying the volume transpose, then add the boundary transpose to each affected field. Interior rows remain untouched. There is no surface measure or quadrature weight in this strong condition. The model's assembled operator performs this composition; a standalone volume operator alone does not.

The pressure regressions retain the MeanField residual/action comparisons at 1e-13, the analytic closure/density-balance thresholds of 2e-12, and the mass threshold of 2e-11 * target_mass. The constrained density-balance finite difference retains step 1e-5 and threshold 2e-6. Analytic boundary values and independent MFEM boundary masks check the shared graph on tetrahedra, hexahedra, both vector orderings, and a curved star. General surface tests also cover non-enthalpy EOS relations, state-dependent couplings, transpose actions, cache invalidation, rejection, and recovery.

Selective serial preparation

Preparation defaults to PreparationMode::reuse_unchanged, matching the generic dependency engine. Request a full refresh explicitly with PreparationMode::rebuild:

backend::prepare_volume(closure, closure_discretization, term_coefficients, displacement);
const auto report = backend::volume_work_report(closure);
// report.geometry_points, sampled_operands, physics_points and reused_elements

backend::prepare_volume(
    closure, closure_discretization, term_coefficients, displacement,
    backend::PreparationMode::rebuild
);

The dependency engine owns exact coefficient snapshots and accepted results. make_volume_inputs supplies MeshAndCoefficientBindings, displacement, compactification, and the equation plan's coefficient terms as source roles. CoefficientTracking identifies each source by its persistent binding and compares coefficient contents. Temporary context wrappers and copied equal vectors therefore retain valid work. A changed operand refreshes its samples and the physics of contributions using that operand. Other contributions retain their accepted data. Displacement or compactification changes refresh mapping and mapped bases while retaining unaffected fixed-mesh coefficient samples. A fixed-mesh coordinate change rebuilds those bases too. No coefficient hashes or caller-maintained coefficient revision counters are required.

Reuse assumes fixed equation configuration and any external data it references. Request PreparationMode::rebuild, or invalidate the cache, after modifying external EOS data; changes hidden behind a runtime EOS view cannot be discovered automatically. Changing owned equation parameters requires a new graph. preparation_count() counts successful generic cache preparations, including reuse and partial requests. Failed attempts do not increment it; invalidation preserves the count. Numerical work counts are derived from the most recent graph report, rather than separately stored state.

Inspect any cache directly through the dependency API:

namespace stages = backend::volume_keys;
const auto& graph_report = closure.report();
const auto basis_status = graph_report.status<stages::Basis<0>>();
const auto geometry_status = graph_report.status<stages::Geometry<0>>();
const auto physics_status = graph_report.status<stages::QuadratureLinearization<0>>();

Statuses are computed, reused, or not_required. The numerical work report summarizes these statuses in element and quadrature counts. Volume helpers request the whole graph; volume_keys::Linearization retains accepted contributions and the assembled residual together. Direct dependency requests can prepare selected intermediates, but application requires the final linearization to be available.

Rejection makes current cache results unavailable. The generic engine keeps its last accepted storage for possible later reuse; backend preparation helpers additionally invalidate that storage on validation or numerical failure, ensuring a fresh retry. Neither path publishes a partly prepared result. Reacquire accepted data after every successful preparation rather than retaining a reference across cache updates.

Reuse depends only on unchanged tracked inputs. A recomputed node advances its version even when it happens to produce the same value, so required downstream computations also refresh. There is no output-equivalence test. Coefficient lookup uses the ordinary dependency context; operand-list deduplication uses CombinedKeyList.

Tensor contractions and complete transposes

On eligible H1/L2 and RT hexahedra, the backend verifies Cartesian quadrature ordering, native signed DOF maps, and the factorization against MFEM basis values and derivatives. It then performs successive one-dimensional contractions. H1 gradients and RT divergence use the same factors and transpose routine. Owned tables are shared within one volume graph and keyed by finite-element identity, operation and complete quadrature contents. Temporary or modified rules cannot leave borrowed table data behind.

Small pointwise maps apply MFEM's canonical-to-fixed-mesh transformation and SERiF's additional physical transformation. Tetrahedra, non-Cartesian rules and unverified factorizations keep the dense path. Neither path changes the quadrature rule. Supply VolumeOptions{.contraction_mode = ContractionMode::dense} to the cache factory to use the dense reference path explicitly. tensor_element_count(linearization) reports how many accepted contribution-elements use tensor contractions, including geometry.

apply_jacobian_transpose(linearization, residual_dual, workspace, packed_result) returns all coefficient duals in Plan::Columns order followed by displacement. JacobianOperator also implements MFEM's MultTranspose. The transpose includes weight, explicit position, operand mapping, and test mapping derivatives. Geometry is differentiated through twelve local inputs: three displacement values and nine canonical derivatives. The mapping module supplies their exact directional variations; no finite differences or global Jacobian-column probing are used in production.

A strong nodal replacement requires the following transpose order:

backend::mask_residual_dual(boundary, residual_dual, volume_dual);
backend::apply_jacobian_transpose(volume, volume_dual, workspace, packed_dual);
// Obtain the constrained term's block from packed_dual using the volume plan's column order.
backend::add_jacobian_transpose_rows(boundary, residual_dual, constrained_term_dual);

Mask residual rows before the volume transpose, then add the constraint transpose to the appropriate coefficient block. Zeroing coefficient entries after the volume transpose would incorrectly remove contributions from unconstrained residual rows. For centering, the constrained coefficient block is the displacement geometry block. These remain serial coordinate transposes; no mass-matrix identification of duals or MPI true-DOF ownership is implicit.

Tensor tests retain the old scalar 2e-12 and gradient/RT 2e-11 comparisons. Full equation transposes are checked against dense contractions, individual Jacobian columns, the old mass-adapter 2e-12 pairing threshold, both element types/orderings, and curved compactified geometry. Existing tests and numerical tolerances are unchanged.

Restrictions and what depends on them

Boundary	Enforcement	Why the numerical code needs it
Authoring interface	equations::Equation	Physics supplies the declared pointwise input/output contract. This does not promise execution support.
Current volume algorithm	operators::VolumeEquation	Every contribution needs a volume-domain support and implemented test/input operations. ScalarVolumeEquation remains available for existing callers.
Requested operand evaluation	OperandEvaluator<Operand> specialization; EvaluableOperand queries its interface	There is no generic fallback that ignores the operation tag. An unsupported operation cannot silently become VALUE sampling.
Actual MFEM element used for sampling	Always-active range/map checks in sample_basis; gradients and H1 divergence additionally require a transformation	Only scalar VALUE bases (one or three copies) are sampled.
Actual space bound to a term	TermBinding: family/rank compile-time checks; dimension, continuity, order, range, mapping runtime checks	Coefficient layout and basis representation must match the declared term. A valid binding still needs a supported operation; vector H1 VALUE and DIVERGENCE are implemented; GRADIENT and CURL remain unsupported.
Supported DOF mapping	DofMap checks serial conforming 3D spaces, sizes, support, orientation, and revisions	Extraction and assembly assume this representation.
Selected boundary rows	select_boundary_dofs and DofSelection::validate	Boundary tags must exist in the schema, have supported DOFs, and remain unchanged during use.
Cross-binding consistency	make_volume_inputs and the basis stage check schema, mesh, binding identities, and input support	Every operand must cover its contribution's integration support, which may be narrower than the row. No operand is silently extended into vacuum.
Inadmissible geometry or physics	Element preparation returns PreparationRejection with element and point; optional integral validation rejects at contribution level	A nonlinear trial state can be rejected without publishing incomplete data. Integral rejections have no single element/point.
Missing or stale preparation	Dependency result availability and validate_volume_data	Cached matrices and geometry must match the borrowed discretization.

Malformed bindings and coefficient sizes are programming errors. Numerical trial-state rejections are values returned from try_prepare_volume. The corresponding prepare_* convenience functions translate a rejection for the caller. Physical admissibility remains the equation/EOS author's responsibility.

Extending the implementation later

For new physics expressible using supported operands, change only the equation and model selection. Keep local physical partial derivatives in the equation; mesh and mapping derivatives belong to numerical execution.

For a new space representation, start with TermBinding and the existing space descriptors. For a new operation, start with OperandEvaluator: its result shape, basis evaluation, coefficient action, test contraction, and geometry dependence must agree. Then examine the element residual and its Jacobian together. In particular, the existing gradient and H1 divergence implementations show why both operand and test variations are needed. A new operation must provide these mathematics before its capability check can admit it.

The volume admission contract should expand only after those calculations, point/cache shapes, and quadrature handling exist. Unsupported operations remain compile-time failures. There is no separate hand-maintained capability registry.

Compile-time checks in operators/plan.cpp cover supported scalar, gradient and H1 divergence operations, sampled result shapes, packed columns, and rejection of unsupported operations. Authoring an equation remains distinct from executing it with this backend. The existing equation tests retain independent quadrature comparisons, analytic dilation checks, centered differences, ownership checks, and their original tolerances.

Gravity without a separate assembly path

primal_gravity.hpp returns its weak Poisson residual as a tuple of contributions: the potential-gradient flux on all domains and the density source on StellarDomains. Contribution<TestOperand, Support> is an optional authoring helper; each contribution supplies Operands, evaluate, and additional_order. All contributions share the row term, but may select different test operations and supports. Existing single-contribution equations need no wrapper. The source's StellarDomains support prevents reading density on vacuum elements. This is volume-operator support, not a new boundary-condition API.

API migration

Previous name	Current name
equations::backend::mfem::EquationField	discretization::backend::mfem::TermBinding
equations::backend::mfem::FieldDofs	discretization::backend::mfem::DofMap
equations::backend::mfem::FieldState	operators::backend::mfem::CoefficientView
EquationState, State	make_term_coefficients(bind(...), ...) returns a dependency context
PreparedEquation, PreparedOperator, VolumeComputations	make_volume_equation_cache, prepare_volume, and free application functions
equations::EquationPlan	operators::EquationPlan
EQWorkspace, EQElementData	operators::backend::mfem::Workspace, ElementData

bind(binding, coefficients) now produces a view with .binding and .values. Pointwise equation inputs and physical derivatives retain their existing interface.