Limits, build time, and worked examples

Known limitations

Almost every boundary below is an explicit, fail-loud NotImplementedError, raised when you assemble a weak form or solve a transient problem through the time route. The residual-PINN path is unaffected.

Two limits are silent — every other one raises

The assembler cannot distinguish these from a legitimate modelling choice, so nothing stops you.

Affine geometry on a curved boundary — by default the mesh is straight-sided, so the domain is approximated to O(h²), which caps every element order above it. The solve is simply suboptimal. shape.curved() fixes it (order 2, simplices) — see Curved geometry.
The 2πr measure on an axisymmetric vector form — exact for scalars, wrong for vectors.
Host geometry inside a jno.derived rule on a moving mesh — ray tables, view factors and neighbour lists are built once, outside the rule, and closed over. They are correct for the mesh they were built on and are never rebuilt; on a mesh that moves or is remeshed they silently describe the old one. Rebuild them and build a new jno.fem.

All three are detailed below, and stated again at the point where you make the choice.

Area	The limit	How it fails
Transient mass, non-nodal	a trainable network coefficient on `u_t * phi` — the mass block is assembled once. A scalar parameter, a coordinate function, or a net on a nodal form are all fine	raises
Temporal order inside a group	a group nested in a product CHAIN (`rho * c(T) * (T_t + u.grad T) * phi`) or handed to an OPERATOR (`inner(u_t + K*u, v)`) stays atomic, and its spatial part is stripped into the mass matrix	raises
A march that never solves	an ABSOLUTE `atol` above the form's own residual scale passes Newton's first test at the incoming iterate, so every step returns its input	raises (the residual is bit-identical before and after)
`d.tag(name, pred)` over a geometry face	the predicate filters point sampling only; facets still map through the geometry, so naming a DIFFERENT face splits the region	warns
Second order in time	nodal Lagrange only, 1-D/2-D/3-D, scalar or vector; the temporal side must stay linear	raises
Reduced-order `basis=`	steady + first-order transient only	raises
Runtime Dirichlet parameters	steady linear, steady nonlinear, linear transient	raises
Time-varying Dirichlet `g(x, t)`	first-order transients, linear or nonlinear, single-field or coupled, a trainable parameter inside the value included (differentiable); not a trainable net inside the value, nor a net/parameter-valued Dirichlet beside it, not a parameter inside it on a `u_tt` form or the τ load path, not a nonlinear `u_tt` form. `sdirk(3)` and `ros2` lose some order to it, and `ros2` holds the wall value only to its order (measured below)	raises
Slip `n·u = 0` on a surface moved by `.trainable()` coordinates	steady nonlinear only (its `P` is rebuilt per solve there)	raises
Affine parameter lowering	one trainable scalar per additive term, not nested	raises
Enclosure radiation	2-D / axisymmetric, needs a direct solve; you write the radiosity yourself	manual composition
`jno.derived` nonlocal fields	nodal Lagrange `on=` only; no chained derived fields; host geometry closed over by the rule is frozen at build, so it goes stale on a moving mesh	raises, except the moving mesh — silent
Plasticity	small-strain, isotropic, linear-hardening	raises
`symmetric=True` matrix fields	native 2-D/3-D Lagrange only — not 1-D, VPINN, periodic ties, complex, or mixed with a non-nodal family; a time-varying whole-tensor Dirichlet value is not mapped to the stored values (write it per entry); an asymmetric whole-tensor value is refused	raises
Constant unknowns (`constant=True`)	native 2-D/3-D Lagrange only; a tie `u(region) - U` on steady linear/nonlinear forms only, and not together with a periodic tie, a moving slip surface, a node that both slips and is tied in the same component, or a hanging-node mesh; `U - g` takes a value that does not depend on position; `jno.fdm` has no constants	raises
`d.unknown(...)` in `jno.fdm`	a real P1 nodal unknown with full storage; `order=2`, a non-nodal `space=`, `complex=True`, `symmetric=True` or `constant=True` make it FEM-only	raises
Adapting a moving mesh	`relocate` fails on two separate bodies (use `alpha`); `enrich()` (p) is not wired; `checkpoint=` on a moving mesh needs `adapt=`, and a checkpointed march cannot be differentiated; runtime connectivity is affine-simplex only	raises
Interpolation covers (`space="cover"`)	first order, simplices only; the layout is padded so memory scales by `1+dim` even where enrichment is off; `jno.solve.enrich` is steady-only; a `u.gap` contact search may not read a cover field	raises
FEM functionals (`fem.eval(F, u)` with a test-free `F`, `expr.integrate(fem)`)	steady, native Lagrange; the whole volume or a tagged boundary region — not a sub-region; not transient, complex, 1-D, non-nodal or VPINN (use `sum(fem.eval(F * phi, u))` there)	raises
VPINN (network trial)	steady only, single field (scalar or vector), no periodic ties; a boundary coefficient must carry a coordinate	raises
`dom.cell_size` / `dom.cell_metric`	native 2-D/3-D volume terms only — a 1-D form or a non-nodal family packs no element Jacobian	raises
Element order on RT / N1E / P0 / Hermite / Argyris / Morley	each family has one intrinsic order	raises
`eigs` on a non-symmetric pencil	eigenvalues differentiate, eigenvectors do not	NaN, not a silent zero
Curved-boundary geometry	straight-sided by default; `shape.curved()` is the fix	silent
Axisymmetric vector forms	the `2πr` measure is wrong for vectors	silent

The detail

Second order in time (u_tt) — what composes and what refuses

A nonlinear spatial operator (sine-Gordon, cubic Klein–Gordon, large-deformation elastodynamics) is supported — Newton on the augmented [u; v] block. The temporal side must stay linear: a state-dependent mass or damping c(u)·u_tt is refused, since M2/C are extracted by differentiating at u=0 and would otherwise be frozen there.

A coupled 2-D/3-D system flows through the same assembler as a single field, so damping, the nonlinear path and a driven boundary g(x,t) all compose with coupling. A coupled form still refuses:

a field with no u_tt term — its velocity rows would be singular;
runtime parameters — the parametric coupled steady assembly underneath is not wired;
periodic ties.

Time-varying Dirichlet is refused on nonlinear u_tt forms (a first-order nonlinear transient takes it). A nonlocal jno.Coupling is refused on any u_tt form, and on a single-field transient with a periodic tie; both routes used to drop it silently. A complex coefficient on any u_tt form is refused by name — it used to be silently cast to real. Write the problem first-order in time instead; the complex transient is supported.

A coupled 1-D system carries u_tt on narrower terms (linear, undamped): the augmented state is [u_all; v_all], so fem.offsets lists the displacement blocks then the velocity blocks.

A trainable net on a non-nodal transient mass — measured boundary

The mass block of a non-nodal transient form is assembled once, so a coefficient that trains cannot enter it. The limit is narrower than it sounds — this is the whole map, measured:

element	`c` in `c * ui.t * vi`
P1 nodal	constant, `jno.np.parameter`, `1 + xi`, trainable net	all solve
Hermite / non-nodal	constant, `jno.np.parameter`, `1 + xi`	all solve
Hermite / non-nodal	trainable net	raises

u, phi = d.fem_symbols(space="Hermite")
rho = 1.0 + jno.nn(net)(xi, yi)                      # trainable, spatially varying
jno.fem([rho * ui.t * vi + (ui.x * vi.x + ui.y * vi.y), ...])
# NotImplementedError: jno.fem (non-nodal): a trainable neural coefficient on the mass (u_t)
# term is not supported -- the mass block is assembled once. Use it on spatial terms, or
# .freeze() the network.

jno.fem([ui.t * vi + rho * (ui.x * vi.x + ui.y * vi.y), ...])   # instead: on the stiffness

Reduced-order solves — what basis= covers

Steady and first-order transient, linear and nonlinear. Second-order-in-time (u_tt), complex, and periodic-tied problems each refuse with their own reason. Nonlinear reduces, but without hyper-reduction that is a memory win, not a speed one.

Runtime Dirichlet parameters — where the gradient flows

A trainable jno.np.parameter may sit in an essential value (u(top) - g, or scaling a coordinate profile u(top) - g*sin(pi*x)), and the boundary value is recovered from data like any other parameter — ∂b/∂g flows through the symmetric elimination (linear), and ∂/∂g through the solve's / each step's custom_root (nonlinear / transient).

A value that is both parametric and t/τ-dependent is held at each step's t/τ WITH the runtime parameters, on a first-order transient and on a τ load-path march (u(top)[1] - g*tau, displacement control with a trainable grip). Refused loudly on a second-order (u.tt) form, whose block evaluates the value and its rate without them: train that amplitude through a Neumann / body term instead. A parametric wall beside a τ-ramped one (u(left) - g with u(right) - delta*tau) is fine on the load path: both are held at every step. A t/τ-dependent essential value needs steps to be held at: on a steady form that reads no step history (.i(k)), and so does not march its domain(tau=...) grid, it is refused, with or without a parameter in the form. A FIELD-sized optimizer-less parameter stays the nodal data-field value (a neighbour's field in a DD solve), gathered per node.

Time-varying Dirichlet g(x, t) — the orders each scheme keeps

The constrained rows are imposed at the time each step or stage lands on, on linear and nonlinear, single-field and coupled first-order transients alike (see Boundary conditions). What that costs in accuracy depends on the scheme. Measured on the Taylor–Green vortex (Navier–Stokes, Taylor–Hood P2/P1, ν = 1, T = 0.5, h = π/6): the walled box [0, π]² has the exact velocity imposed on all four walls; the periodic box [0, 2π]² has no Dirichlet data at all. Velocity error against a same-mesh 256-step ros34pw2 reference (so the spatial error cancels), rates over 8 → 16 → 32 steps:

scheme	walls driven by `g(x, t)`	periodic box	wall rows exact?
`theta(1.0)`	0.98, 0.99	0.96, 0.98	yes
`theta(0.5)`	2.18, 2.01	2.10, 2.01	yes
`bdf2()`	2.23, 1.82 (→ 1.96 by 128 steps)	2.10, 2.04	yes
`sdirk(2)`	2.01, 2.00	2.01, 2.00	yes
`sdirk(3)`	2.63, 2.60	2.95, 2.98	yes
`rosenbrock()` (ros34pw2)	2.87, 2.93	2.95, 2.98	yes
`rosenbrock("ros2")`	1.48, 1.68	1.77, 1.87	no — 5.3e-4 off at 16 steps

The losses in bold are the classical order reduction of one-step methods with low stage order under time-dependent boundary data (Ostermann & Roche, Math. Comp. 59 (1992) 403–420), not a mis-timed boundary value: the same form with the convection dropped, marched as a linear block, gives rates identical to three digits, and the one-step-late control drops BDF2 to first order (1.05) at ~900× the error. "Wall rows exact" means the nodal wall value equals g(x, tⁿ) to round-off at every step: a stiffly accurate scheme's last stage is the step, so it lands on g(tₙ₊₁); ros2 is not stiffly accurate, so its step is a combination of stages and holds the wall value only to its order (1.4e-3 off after the first step, 5.3e-4 at T, falling as dt²). Re-imposing g(tₙ₊₁) on the wall rows after each step was tried and not adopted: the wall becomes exact, but the change it makes to the interior decays only as O(dt) — on the walled Taylor–Green box the ros2 rate over 32 → 64 → 128 → 256 steps fell from 1.81, 1.89, 1.94 to 1.77, 1.83, 1.82, with 19 % more interior error at 256 steps. Use the default ros34pw2 when a time-varying wall value must hold exactly.

Plasticity — what runs today

Deformation theory (monotonic / proportional) and the path-dependent flow-theory tau= load-path march both run. The march assembles on the real, steady native-Lagrange path, single-field or coupled — not transient / complex / 1-D / non-nodal, each rejected with a clear error. Periodic ties, exact slip conditions and hanging-node constraints compose (each step solves in the reduced space u = P ũ); with one of them, tau=jno.solve.arclength(...) refuses by name. A .bounds(...) box composes with a periodic tie whose bound is the same on both faces, and refuses beside a weighted elimination (mortar, hanging nodes, slip).

An internal state advances on every cell by default; state.evolves(formula, region="strip") restricts it to one region, and outside it the state is frozen at the value it already has — so a zero-initialised plastic strain leaves the rest of the body elastic without a second constitutive branch, and the return map is never evaluated on cells carrying another material's constants. That is what makes a plastic part bearing on an elastic one expressible.

Kinematic / nonlinear (Voce) hardening are separate formulas, not built.

Curved boundaries: P3 buys nothing over P2

With the default straight-sided mesh the domain itself is approximated to O(h²), and that error caps every element order above it. Measured on the unit disk (-Δu = 4, u|∂Ω = 0, exact u* = 1 − r²), L2 rates under refinement:

order	expected	measured	error at `h = 0.05`
P1	2	2.00	1.14e-03 (1 536 dofs)
P2	3	2.02	7.46e-04 (6 015 dofs)
P3	4	2.01	7.40e-04 (13 438 dofs)

P3 buys nothing over P2 — 13 438 dofs for the same answer — and both are barely ahead of P1. The cap comes from imposing the boundary condition at straight-edge nodes that sit on the chord, O(h²) inside the true arc.

On a polygonal domain the advertised rates hold exactly (the suite measures P2/P3 there).

The fix on a curved one is shape.curved(), which places the midside nodes on the true surface and recovers P2's own third order — 570× more accurate at the finest resolution in the same study. It is order 2 and simplices only, and non-nodal families keep affine geometry, so where it does not apply the advice above still stands: prefer h-refinement (or the adaptive loop) over order ≥ 2. See Curved geometry.

Interpolation covers — the padding, and the Dirichlet trace

A cover field's node array is the mesh nodes repeated 1+dim times, so memory scales by that factor even where enrichment is switched off — an unenriched node's slots exist and are pinned to zero (the Dirichlet mechanism, eliminated exactly). That is the price of the padded layout, and it is what buys a rebuild that changes no points and no connectivity.

An inhomogeneous Dirichlet trace stays the P1 interpolant of g: the tangential covers pin to zero, not to dg/ds. On a curved or non-constant boundary value that is a real accuracy limit, not a formality.

jno.solve.enrich(...) is wired on the steady loop only — enriching mid-march would change the DOF layout under the stepper, and the state transfer that carries a solution across an adapt round is written for a change of mesh, not of space. A transient problem is refused by name, as is order= (the cover supplies the extra order) and a non-simplex cell (the cover gradient assumes an affine map). See Interpolation covers.

VPINN (network trial) — where the lowering stops

A network trial is test-projected onto the FE basis, and that projection is spatial — so a form carrying the time coordinate is refused. It used to build and evaluate on a meaningless number: the time grid never reached the residual and the initial condition was silently dropped (measured: u(initial) - 0 and u(initial) - 7 gave bit-identical residuals). Making it real needs a decided treatment of time — space-time test functions, or collocation on the declared grid with the IC as its own loss — which is a design choice rather than an implementation gap.

3-D works. It used to raise -- first on an internal quadrature-pool tag, then by name -- and the restriction turned out to be a conservative guard rather than a limitation: the native assembler builds the context in 3-D exactly as in 2-D. Measured on the cube: Poisson 4.6e-04, a Neumann flux face 8.4e-04, with the same scalar and vector operator surface as 2-D.

One spelling is required on a boundary term: its coefficient must carry a coordinate from that region ((g + 0.0 * xr) * v_r). A bound test keeps its binding on the view rather than in the expression tree, so a bare constant is indistinguishable from a volume term once the form is flattened, and would be integrated over the volume silently. It is refused instead. It is also single-field: a vector unknown is fine (one field with value_shape=(d,), and all four source spellings lower identically), but a coupled multi-field system is not. Periodic ties are refused because a tie is an algebraic reduction of FE trial DOFs and a network trial has none — impose periodicity inside the network instead.

Each of the three refuses by name at the point the decision is made. They used to surface further in: a 3-D VPINN died on Tag 'fem_gauss' is not in the mesh pool or context, an internal quadrature-pool name for a scope limit nobody could infer from it.

What does compose is on Formulations, measured against analytic solutions — including a network trial co-trained with a jno.np.parameter.

Element order on a non-nodal family — refused, not applied

RT / N1E / P0 / Hermite / Argyris / Morley each have one intrinsic order. space="N1E", order=2 used to return the same lowest-order space silently; it now raises. Note where it raises: fem_symbols accepts the argument, and jno.fem([...]) is what refuses — the order is not checked until an element is actually built.

u, v = d.fem_symbols(space="N1E", value_shape=(2,), order=2)   # accepted
ui, vi = u.bind(x=xi, y=yi), v.bind(x=xi, y=yi)
jno.fem([inner(ui, vi) - inner(e_x, vi)])                       # <- raises here
# NotImplementedError: jno.fem: order=2 is not selectable on the N1E element — its order is
# intrinsic to the family (lowest order (N1E0)) and jNO builds only that one. Drop order= (the
# default) to get it.

u, v = d.fem_symbols(space="N1E", value_shape=(2,))             # instead: take the family's order

The mesh is the only accuracy knob on an H(curl)/H(div) problem — see Mesh resolution for wave problems for what a given points-per-wavelength buys, measured.

A complex essential value — refused in 1-D, 2-D and 3-D

A complex form is solved as two coupled real legs that share one Dirichlet row set. That set can impose Re u = g with Im u = 0, but not a prescribed Im u, so a complex g is refused rather than silently truncated:

fem = jno.fem([ui.x * vi.x + ui.y * vi.y - k**2 * (1 + 0.1j) * ui * vi - 1.0 * vi,
               u(xb, yb) - (0.5 + 0.25j)])
# NotImplementedError: jno.fem: a COMPLEX essential value on region 'boundary' is not supported
# — a complex form's Re/Im legs share one Dirichlet row set ...

fem = jno.fem([..., u(xb, yb) - 0.5])       # instead: a REAL essential value

The restriction is on the boundary value, not on the problem: the operator and the source stay fully complex, and the solution has an imaginary part everywhere (measured max|Im u| = 2.6e-02 on the 2-D Helmholtz above with g = 0.5). Carry the complex part there.

eigs routing, differentiability, and when linear= pays

jno.solve.eigs / FEM.eigs route on the operator's actual symmetry. A symmetric pencil takes the symmetric reductions (real spectrum, differentiable). A genuinely non-self-adjoint operator goes to ARPACK/Arnoldi (Lehoucq & Sorensen 1996) and returns the complex spectrum it actually has — the case that matters for stability problems (resistive MHD growth rates, drift waves, anything with a mean flow), where the sign of the growth rate is the physics. Neither path ever returns the spectrum of ½(K+Kᵀ) as though it were the answer.

The routing probe is a randomized bilinear test and is concrete-only, so under jit the symmetric path is assumed.

Differentiability of the non-symmetric path. The eigenvalues are differentiable in reverse mode — dλ = wᴴ(dA − λ dB)v / (wᴴBv) for a simple eigenvalue (Wilkinson 1965 ch. 2), with the left eigenvector obtained by inverse iteration on (A − λᵢB)ᴴ, verified against finite differences to 1e-09. The eigenvectors are not, and differentiating through them yields NaN rather than the silent zero a plain callback would give. A defective eigenvalue has no derivative at all (its perturbation series runs in √ε) and is detected via the eigenvalue condition number, giving NaN instead of the enormous finite number the formula would produce. Because that guard is a custom_vjp, forward mode (jax.jvp/jacfwd) is unavailable — use jax.grad / jacrev.

linear= is opt-in because it is not always a win. It accepts jno.solve.lu(backend="pardiso"/"cudss"/"host") to drive ARPACK's shift-invert factorization — those kernels are plain numpy functions, so ARPACK can call them from host code — but refuses backend="device" (a JAX primitive), the Krylov solvers, and precond=, rather than ignoring them. ARPACK applies the inverse ~50–70 times, so it trades one fast factorization against per-application overhead. Measured with PARDISO:

pencil size	speed-up
n = 3,000	0.72× (slower)
n = 20,000	10.05× (21.6 s vs 217 s)

It needs k < n-1; smaller pencils take an exact dense scipy.linalg.eig. Order with which="LR" / "SR" (real part — the growth rate) or target an interior region with sigma=.

Axisymmetric vector forms are your responsibility

The 2πr measure is exact for scalars and wrong for vectors — elasticity hoop strain, and for vector Maxwell the cylindrical curl's own 1/r terms plus the meridional/azimuthal decoupling.

For elasticity and flow the extra terms are short, and the recipe is written out and verified in The vector recipe: add ε_θθ = u_r/r to the strains and to the trace (so div u = ∂ᵣu_r + u_r/r + ∂zu_z). It reproduces the Lamé thick-walled cylinder to under 1 % and Poiseuille pipe flow to machine precision; dropping the hoop term is over 10 % wrong and raises nothing (tests/test_fem_axisymmetric_vector.py). Vector Maxwell is the case with no recipe here: use a full 3-D mesh. jNO ships no axisymmetric H(curl)/H(div) element, and multiplying by r is arithmetic the assembler cannot distinguish from a legitimate radial coefficient, so nothing raises.

Use a full 3-D mesh for vector Maxwell, or write the cylindrical operator out yourself. See Axisymmetric (bodies of revolution).

Build time: what to expect, and the one knob

jno.fem([...]) fully assembles the operator — it returns concrete matrix values, which is why fem.solve() is then only the linear solve (measured 7 ms on 3-D Poisson at 27,833 nodes).

Some libraries instead defer assembly into their solve. That makes their "build" look faster and their solve slower, so compare build + solve — and remember jNO assembles once, while a per-solve assembler pays again on every Newton iteration.

Most of a cold build is XLA compilation, and the cost is fixed per problem structure rather than per DOF: a 15× larger mesh still compiles about the same number of programs. Two caches attack it, both on by default.

Across processes — the persistent XLA cache

Enabled at import jno, stored in ~/.cache/jno/xla.

3-D Poisson, 27,833 nodes

	first build	repeat build
no cache	4.75 s	2.48 s
persistent cache (default)	2.22 s	1.51 s

The very first run on a machine is slower — populating the cache costs more than not having one. It pays back from the second process onward, which is jNO's normal life: sweeps, optimisation loops, test suites, re-running a script after an edit.

Opting out

Any one of: JNO_COMPILE_CACHE=0, jno.setup(__file__, compile_cache=False), or per project in .jno.toml:

[jno]
compile_cache = false

Within a process — identical-problem reuse

Rebuilding an identical problem (same mesh content, same terms) reuses the already-compiled assembly kernels outright, keyed on content rather than object identity. A rebuild then costs meshing plus host prep and no XLA work at all — measured 1.94 s → 0.28 s on 3-D Poisson at 29k nodes.

Anything that changes the operator — a different mesh, a different coefficient — recompiles exactly the kernels that bake it. Structure a tokenizer cannot key by value simply never caches (a safe miss); the coverage is measurable, not guessed (jno.utils.solver.fem_utils._ELEM_MAP_STATS).

Worked examples

The FEM tutorials cover every pattern in this guide:


Elliptic	Poisson · mixed Dirichlet/Robin reaction–diffusion · mixed-BC Helmholtz
Nonlinear	Allen–Cahn interface
Solid & fluid	linear-elastic cantilever · Stokes channel flow
Transient	heat · two second-order-in-time wave examples (vibrating membrane, vector-elastodynamics cantilever)
Inverse	a hidden diffusivity field · a transient rate
Non-nodal	H(div) mixed Poisson · H(curl) Maxwell / eddy currents
Neural	a variational PINN — a neural trial in the same weak form