Skip to content

Differential operators¤

Note

For a more detailed mathematical guide (notation, shapes, and backend behavior), see Guides → Differential operators.

Backends (AD / jet / FD / basis)¤

Many differential operators accept a backend keyword:

  • backend="ad": autodiff (works for point batches and coord-separable batches).
  • backend="jet": Taylor-mode AD ("jets") for higher-order derivatives with respect to a single variable. See the guide for mathematical details (including Faà di Bruno / Bell polynomial structure).
  • backend="fd": finite differences on coord-separable grids (falls back to AD on point batches). Use periodic=True for periodic stencils.
  • backend="basis": basis-aware methods on coord-separable grids (falls back to AD on point batches). The basis keyword selects a 1D method per axis:
  • basis="fourier": FFT-based spectral derivatives on periodic grids,
  • basis="sine" / basis="cosine": FFT-based spectral derivatives via odd/even extension (sine \(\leftrightarrow\) odd extension, cosine \(\leftrightarrow\) even extension),
  • basis="poly": polynomial (barycentric) derivatives for generic 1D grids.

Note

FFT-based bases (fourier/sine/cosine) assume a uniformly-spaced coordinate axis. For non-uniform axes, prefer basis="poly" or backend="ad".

Common keyword arguments¤

Many operators share these keywords:

  • var: label to differentiate with respect to (e.g. "x" or "t"). If omitted, the variable is inferred when possible.
  • mode: autodiff mode ("reverse" uses jax.jacrev, "forward" uses jax.jacfwd) for AD-based paths.
  • backend: "ad", "jet", "fd", or "basis" (see above).
  • basis: basis method used when backend="basis".
  • periodic: periodic treatment for FD stencils (used when backend="fd").

See the guide for operator shape conventions and for the math behind surface and fractional operators.

Core derivatives¤

phydrax.operators.grad(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Gradient/Jacobian of u with respect to a labeled variable.

For a geometry variable \(x\in\mathbb{R}^d\) this constructs \(\nabla_x u\). Concretely:

If \(u:\Omega\to\mathbb{R}\) is scalar-valued, then grad(u) returns a vector field in \(\mathbb{R}^d\):

\[ \nabla_x u = \left(\frac{\partial u}{\partial x_1},\dots,\frac{\partial u}{\partial x_d}\right). \]

If \(u:\Omega\to\mathbb{R}^m\) is vector-valued, then grad(u) returns the Jacobian \(J\in\mathbb{R}^{m\times d}\) with entries \(J_{ij}=\partial u_i/\partial x_j\).

For a scalar domain variable (e.g. time \(t\)), this reduces to the partial derivative \(\partial u/\partial t\).

Arguments:

  • u: Input function.
  • var: Variable label to differentiate with respect to. If None, then the domain must have exactly one differentiable variable.
  • mode: Autodiff mode: "reverse" uses jax.jacrev, "forward" uses jax.jacfwd.
  • backend: Differentiation backend.
  • "ad": autodiff (works for point batches and coord-separable batches).
  • "fd": finite differences on coord-separable grids (falls back to "ad").
  • "basis": spectral/barycentric methods on coord-separable grids (falls back to "ad").
  • basis: Basis method used when backend="basis".
  • periodic: Whether to treat the differentiated axis as periodic (used by backend="fd" tuple paths).

Notes:

  • Coord-separable inputs (tuples of 1D coordinate arrays) are handled by taking JVPs along each coordinate axis and stacking.

Returns:

  • A DomainFunction representing the gradient/Jacobian. For geometry variables, the derivative axis is appended on the right:
  • scalar \(u\): grad(u) has trailing shape (..., d);
  • vector-valued \(u\) with trailing size \(m\): trailing shape (..., m, d).

Example:

import phydrax as phx

geom = phx.domain.Interval1d(0.0, 1.0)

@geom.Function("x")
def u(x):
    return x[0] ** 2

du = phx.operators.grad(u, var="x")

phydrax.operators.div(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Divergence of a vector field.

For a vector field \(u:\Omega\to\mathbb{R}^d\), returns

\[ \nabla\cdot u = \sum_{i=1}^{d}\frac{\partial u_i}{\partial x_i}. \]

The input is expected to be vector-valued with last axis size \(d\) (the geometry dimension for var). If the value has additional leading axes (e.g. multiple vector fields stacked), divergence is applied componentwise over those axes.

Arguments:

  • u: Vector field to differentiate.
  • var: Geometry variable label. If None, the domain must have exactly one geometry variable.
  • mode, backend, basis, periodic: Passed through to grad(u, ...).

Returns:

  • A DomainFunction representing \(\nabla\cdot u\).

phydrax.operators.curl(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Curl of a 3D vector field.

For \(u:\Omega\to\mathbb{R}^3\), returns \(\nabla\times u\):

\[ \nabla\times u = \begin{pmatrix} \partial_y u_z - \partial_z u_y \\ \partial_z u_x - \partial_x u_z \\ \partial_x u_y - \partial_y u_x \end{pmatrix}. \]

Arguments:

  • u: Vector field (must be 3D-valued).
  • var: Geometry variable label. If None, the domain must have exactly one geometry variable.
  • mode, backend, basis, periodic: Passed through to grad(u, ...).

Returns:

  • A DomainFunction representing \(\nabla\times u\) (vector-valued with trailing size 3).

phydrax.operators.hessian(u: DomainFunction, /, *, var: str | None = None, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Hessian of u with respect to a labeled variable.

For a geometry variable \(x\in\mathbb{R}^d\) and scalar-valued \(u\), this returns the matrix of second derivatives

\[ H_{ij}(x) = \frac{\partial^2 u}{\partial x_i\,\partial x_j}. \]

For vector-valued \(u:\Omega\to\mathbb{R}^m\), the Hessian is taken componentwise, producing a trailing shape (..., m, d, d).

For a scalar variable (e.g. time \(t\)) this reduces to \(\partial^2 u/\partial t^2\).

Arguments:

  • u: Input function.
  • var: Variable label to differentiate with respect to. If None, the domain must have exactly one differentiable variable.

Returns:

  • A DomainFunction representing the Hessian/second derivative.

Example:

import phydrax as phx

geom = phx.domain.Square(center=(0.0, 0.0), side=2.0)

@geom.Function("x")
def u(x):
    return x[0] ** 2 + x[1] ** 2

H = phx.operators.hessian(u, var="x")

phydrax.operators.laplacian(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'jet', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Laplacian of a scalar field with respect to a geometry variable.

For \(u:\Omega\to\mathbb{R}\) and \(x\in\mathbb{R}^d\), this returns

\[ \Delta u \;=\; \nabla\cdot\nabla u \;=\; \sum_{i=1}^{d}\frac{\partial^2 u}{\partial x_i^2}. \]

Arguments:

  • u: Input function (typically scalar-valued).
  • var: Geometry label to differentiate with respect to. If None, the domain must have exactly one geometry variable.
  • mode: Autodiff mode for Jacobian/Hessian construction when applicable.
  • backend: Differentiation backend:
  • "ad": sums second partials via partial_n(..., order=2, backend="ad") using AD derivatives (respecting ad_engine).
  • "jet": uses JAX Jet expansions for second directional derivatives (and enables latent Jet fast paths when available).
  • "fd": uses finite differences on coord-separable grids (falls back to "ad").
  • "basis": uses spectral/barycentric methods on coord-separable grids (falls back to "ad").
  • basis: Basis method used when backend="basis".
  • periodic: Whether to treat differentiated axes as periodic (used by backend="fd").
  • ad_engine: AD execution strategy used when backend="ad":
  • "auto": existing mixed strategy (default).
  • "reverse" / "forward": force Jacobian AD mode.
  • "jvp": force directional-JVP derivatives.

Returns:

  • A DomainFunction representing \(\Delta u\).

phydrax.operators.bilaplacian(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'jet', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Bi-Laplacian of a scalar field with respect to a geometry variable.

Returns the fourth-order operator

\[ \Delta^2 u = \Delta(\Delta u). \]

The backend options match laplacian.

Arguments:

  • u: Input function (typically scalar-valued).
  • var: Geometry label to differentiate with respect to. If None, the domain must have exactly one geometry variable.
  • mode, backend, basis, periodic, ad_engine: As in laplacian.

Returns:

  • A DomainFunction representing \(\Delta^2 u\).

phydrax.operators.dt(u: DomainFunction, /, *, var: str = 't', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Alias for \(\partial u/\partial t\) (partial_t).

phydrax.operators.dt_n(u: DomainFunction, /, *, var: str = 't', order: int, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'jet'] = 'ad', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Nth time derivative \(\partial^n u/\partial t^n\).

This is a thin wrapper around partial_n(u, var=time_var, axis=None, order=n, ...).

Arguments:

  • u: Input function.
  • var: Time label.
  • order: Derivative order \(n\ge 0\).
  • mode: Autodiff mode used when backend="ad".
  • backend: "ad" (repeated application) or "jet" (Jet expansions for \(n\ge 2\)).

phydrax.operators.directional_derivative(u: DomainFunction, v: DomainFunction, /, *, var: str = 'x', mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Directional derivative of u along a direction field v.

For a geometry variable \(x\in\mathbb{R}^d\), this computes

\[ D_v u \;=\; v\cdot\nabla_x u, \]

where \(v=v(x)\) is a vector field and \(\nabla_x u\) is computed by grad.

If \(u\) is vector- or tensor-valued, the directional derivative is applied componentwise, i.e. \(D_v u\) has the same value shape as \(u\).

Arguments:

  • u: Input function.
  • v: Direction field (must be vector-valued with trailing size \(d\) for var).
  • var: Geometry label to differentiate with respect to (must be a geometry variable).
  • mode, backend, basis, periodic: Passed through to grad(u, ...).

Returns:

  • A DomainFunction representing \(D_v u\).

phydrax.operators.material_derivative(u: DomainFunction, v: DomainFunction, /, *, spatial_var: str = 'x', time_var: str = 't', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Material derivative (advective derivative) of a field.

Given a velocity field \(v(x,t)\), the material derivative is

\[ \frac{D u}{D t} = \frac{\partial u}{\partial t} + v\cdot\nabla_x u. \]

Arguments:

  • u: Field to differentiate.
  • v: Velocity field (vector-valued over spatial_var).
  • spatial_var: Geometry label (default "x").
  • time_var: Time label (default "t").
  • mode: Autodiff mode used by dt and directional_derivative.
  • ad_engine: AD execution strategy for composed AD derivatives.

phydrax.operators.partial_t(u: DomainFunction, /, *, var: str = 't', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Convenience wrapper for \(\partial u/\partial t\) (scalar variable var).

Equivalent to partial(u, var=var, axis=None, mode=mode).

phydrax.operators.partial_x(u: DomainFunction, /, *, var: str = 'x', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Convenience wrapper for \(\partial u/\partial x\) (axis 0 of geometry variable var).

Equivalent to partial(u, var=var, axis=0, mode=mode).

phydrax.operators.partial_y(u: DomainFunction, /, *, var: str = 'x', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Convenience wrapper for \(\partial u/\partial y\) (axis 1 of geometry variable var).

Equivalent to partial(u, var=var, axis=1, mode=mode).

phydrax.operators.partial_z(u: DomainFunction, /, *, var: str = 'x', mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Convenience wrapper for \(\partial u/\partial z\) (axis 2 of geometry variable var).

Equivalent to partial(u, var=var, axis=2, mode=mode).

phydrax.operators.partial_n(u: DomainFunction, /, *, var: str, axis: int | None = None, order: int, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'jet', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Nth partial derivative with respect to a labeled variable.

Computes \(\partial^n u / \partial x_i^n\) (for geometry variables with axis=i) or \(\partial^n u / \partial t^n\) (for scalar variables).

Arguments:

  • u: Input function.
  • var: Variable label to differentiate with respect to.
  • axis: For geometry variables, which coordinate axis \(i\) to differentiate along. For scalar variables, must be None.
  • order: Derivative order \(n\ge 0\).
  • mode: Autodiff mode used when backend="ad".
  • ad_engine: AD execution strategy used when backend="ad":
  • "auto": default strategy.
  • "reverse" / "forward": force Jacobian AD mode.
  • "jvp": force directional-JVP derivatives.
  • backend:
  • "ad": repeated application of partial.
  • "jet": uses Jet expansions for \(n\ge 2\) (point inputs and coord-separable inputs).
  • "fd": finite differences on coord-separable grids (falls back to "ad" for point inputs).
  • "basis": spectral/barycentric methods on coord-separable grids (falls back to "ad" for point inputs).
  • basis: Basis method used when backend="basis".
  • periodic: Whether to treat the differentiated axis as periodic (used by backend="fd").

Returns:

  • A DomainFunction representing the requested \(n\)th partial derivative.

Surface operators¤

phydrax.operators.surface_grad(u: DomainFunction, component: DomainComponent, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Surface (tangential) gradient on a boundary component.

Let \(n\) be the outward unit normal and \(P = I - n\otimes n\) the tangential projector. For a scalar field \(u\), the surface gradient is

\[ \nabla_{\Gamma} u = P\,\nabla u. \]

Arguments:

  • u: Field to differentiate (typically scalar-valued).
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label for the boundary variable.
  • mode: Autodiff mode passed to grad.

Returns:

  • A DomainFunction representing \(\nabla_\Gamma u\) (tangential vector field).

phydrax.operators.surface_div(v: DomainFunction, component: DomainComponent, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Surface (tangential) divergence on a boundary component.

With tangential projector \(P = I - n\otimes n\), this implements

\[ \nabla_{\Gamma}\cdot v = \text{tr}(P\,\nabla v), \]

where \(\nabla v\) is the Jacobian of \(v\) with respect to the ambient coordinates.

Arguments:

  • v: Tangential or ambient vector field.
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label for the boundary variable.
  • mode: Autodiff mode passed to grad.

Returns:

  • A DomainFunction representing \(\nabla_\Gamma\cdot v\) (scalar field).

phydrax.operators.surface_curl_scalar(u: DomainFunction, component: DomainComponent, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Surface curl of a scalar field on a 3D surface.

For a scalar field \(u\) on a surface in \(\mathbb{R}^3\), returns the tangential vector field

\[ \text{curl}_{\Gamma} u = n \times \nabla_{\Gamma} u. \]

Arguments:

  • u: Scalar field on the surface.
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label (must be 3D).
  • mode: Autodiff mode used by surface_grad.

Returns:

  • A DomainFunction representing the tangential vector field \(\text{curl}_\Gamma u\).

phydrax.operators.surface_curl_vector(v: DomainFunction, component: DomainComponent, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Surface curl of a vector field on a 3D surface.

For a vector field \(v\) on a surface in \(\mathbb{R}^3\), returns the scalar

\[ \text{curl}_{\Gamma} v = n \cdot (\nabla \times v). \]

Arguments:

  • v: Vector field on the surface.
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label (must be 3D).
  • mode: Autodiff mode used by curl.

Returns:

  • A DomainFunction representing the scalar surface curl \(\text{curl}_\Gamma v\).

phydrax.operators.tangential_component(w: DomainFunction, component: DomainComponent, /, *, var: str | None = None) -> DomainFunction ¤

Project a vector field onto the local tangent space.

Given a unit normal field \(n(x)\) on a boundary component, the tangential projection of a vector field \(w\) is

\[ w_{\tau} = w - (w\cdot n)\,n. \]

Arguments:

  • w: Vector field to project (trailing size = ambient dimension).
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label for the boundary variable (defaults to inferred geometry label).

Returns:

  • A DomainFunction representing the tangential projection \(w_\tau\).

phydrax.operators.laplace_beltrami(u: DomainFunction, component: DomainComponent, /, *, var: str | None = None, curvature_aware: bool = False, variant: Literal['contraction', 'divgrad'] | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Laplace–Beltrami operator on a boundary component.

The Laplace–Beltrami operator \(\Delta_{\Gamma}\) is the surface analogue of the Laplacian. Two common realizations are supported:

  • projection–contraction form (default): \(\Delta_{\Gamma} u \approx \text{tr}(P\,\nabla^2 u\,P)\);
  • divergence-of-surface-gradient (via variant="divgrad"): \(\Delta_{\Gamma} u = \nabla_{\Gamma}\cdot\nabla_{\Gamma}u\).

The curvature_aware=True option selects the divgrad variant by default.

Arguments:

  • u: Field to differentiate.
  • component: Boundary DomainComponent used to supply the unit normal field.
  • var: Geometry label.
  • curvature_aware: If True, defaults to variant="divgrad".
  • variant: "contraction" (projection–contraction) or "divgrad".
  • mode: Autodiff mode used by the underlying grad/surface_grad.

Returns:

  • A DomainFunction representing \(\Delta_\Gamma u\).

phydrax.operators.laplace_beltrami_divgrad(u: DomainFunction, component: DomainComponent, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Laplace–Beltrami operator via surface divergence of the surface gradient.

Implements

\[ \Delta_{\Gamma} u = \nabla_{\Gamma}\cdot(\nabla_{\Gamma}u). \]

Arguments:

  • u: Field to differentiate.
  • component: Boundary DomainComponent.
  • var: Geometry label.
  • mode: Autodiff mode used by surface_grad/surface_div.

Returns:

  • A DomainFunction representing \(\Delta_\Gamma u\).

Fractional derivatives¤

phydrax.operators.fractional_laplacian(u: DomainFunction, /, *, alpha: float, var: str | None = None, radius: float | None = None, num_points: int = 4096, eps: float = 1e-06, desingularize: bool = False, angular_points: int | None = None) -> DomainFunction ¤

Fractional Laplacian via a truncated ball quadrature estimator.

The (spectral/integral) fractional Laplacian in \(\mathbb{R}^d\) admits an integral representation (up to a normalization constant \(C_{d,\alpha}\)):

\[ (-\Delta)^{\alpha/2}u(x) \propto \int_{\mathbb{R}^d}\frac{u(x)-u(y)}{\|x-y\|^{d+\alpha}}\,dy, \qquad 0<\alpha<2. \]

This implementation uses a truncated ball integral with radius radius and a Monte Carlo-style quadrature rule on offsets \(y=x+\xi\):

\[ \int_{B_R(0)} \frac{u(x)-u(x+\xi)}{\|\xi\|^{d+\alpha}}\,d\xi. \]

This routine returns a quantity proportional to \((-\Delta)^{\alpha/2}u\); it does not apply the normalization constant \(C_{d,\alpha}\) and truncates the integral to a finite ball.

Arguments:

  • u: Input function \(u\).
  • alpha: Fractional order \(\alpha\in(0,2)\).
  • var: Geometry label to apply the operator to. If None, Phydrax infers a differentiable geometry label when possible.
  • radius: Truncation radius \(R\) (defaults to an AABB-derived scale).
  • num_points: Number of ball quadrature offsets.
  • eps: Excludes offsets with \(\|\xi\|\le \varepsilon\) to avoid the singularity.
  • desingularize: If True and \(\alpha>1\), subtracts a first-order correction term involving \(\nabla u\) to reduce variance.
  • angular_points: Reserved for future deterministic angular quadrature; currently unused.

Returns:

A DomainFunction representing the (truncated, unnormalized) fractional Laplacian.

phydrax.operators.fractional_derivative_gl_mc(u: DomainFunction, /, *, alpha: float, side: Literal['left', 'right'] = 'right', axis: int = 0, N: int = 100, K: int = 1, x_lb: float | None = None, x_ub: float | None = None, sampler: str = 'sobol_scrambled', caputo: bool = False, var: str | None = None, time_var: str | None = 't') -> DomainFunction ¤

Grünwald–Letnikov fractional derivative (Monte Carlo / GMC estimator).

For \(\alpha\in(1,2)\), the Grünwald–Letnikov derivative can be viewed as a limit of weighted finite differences. This routine implements a Monte Carlo / generalized Monte Carlo (GMC) estimator for left- or right-sided fractional derivatives along a chosen spatial axis.

The derivative is taken with respect to a geometry variable var along coordinate axis axis, using bounds x_lb/x_ub to define the one-sided interval.

Arguments:

  • u: Input DomainFunction \(u\).
  • alpha: Order \(\alpha\in(1,2)\).
  • side: "left" or "right".
  • axis: Spatial axis within var used for the one-sided derivative.
  • N: Number of sub-intervals used to define the step size \(h\).
  • K: Sample multiplier; total Monte Carlo samples is N * K.
  • x_lb: Left bound (required when side="left").
  • x_ub: Right bound (required when side="right").
  • sampler: Sampler name for drawing the GMC random variates.
  • caputo: If True, uses a Caputo-style correction.
  • var: Geometry label to differentiate with respect to (defaults to an inferred geometry label).
  • time_var: Optional time label; when present in u.deps, the derivative is conditioned on the current time coordinate.

phydrax.operators.riesz_fractional_derivative_gl_mc(u: DomainFunction, /, *, beta: float, bounds: Sequence[tuple[float, float]], N: int = 100, K: int = 1, sampler: str = 'sobol_scrambled', var: str | None = None, time_var: str | None = 't') -> DomainFunction ¤

Riesz fractional derivative assembled from one-sided GL derivatives.

For \(\beta\in(1,2)\), this constructs a Riesz-type fractional derivative by summing left and right Grünwald–Letnikov derivatives along each spatial axis and applying the standard normalization factor:

\[ \mathcal{D}^{\beta}u \propto \frac{1}{2\cos(\pi\beta/2)}\sum_{i=1}^{d} \left(D_{+,i}^{\beta}u + D_{-,i}^{\beta}u\right). \]

The per-axis bounds are provided via bounds.

Arguments:

  • u: Input DomainFunction \(u\).
  • beta: Order \(\beta\in(1,2)\).
  • bounds: Sequence of per-axis bounds [(x_lb_0, x_ub_0), ..., (x_lb_{d-1}, x_ub_{d-1})].
  • N: Number of sub-intervals used to define step sizes \(h\) per axis.
  • K: Sample multiplier; total Monte Carlo samples per one-sided derivative is N * K.
  • sampler: Sampler name for drawing the GMC random variates.
  • var: Geometry label to differentiate with respect to (defaults to an inferred geometry label).
  • time_var: Optional time label; when present in u.deps, the derivative is conditioned on time.

phydrax.operators.caputo_time_fractional(u: DomainFunction, /, *, alpha: float, time_var: str = 't', mode: str = 'auto', sampler: str = 'sobol_scrambled', order: int = 64, tau_epsilon: float | str = 'auto', cluster_exponent: float | None = None) -> DomainFunction ¤

Caputo fractional derivative in time.

For a time interval starting at \(t_0\), the Caputo derivative of order \(\alpha\in(0,2)\) is defined (for sufficiently smooth \(u\)) by

For \(0<\alpha<1\):

\[ {}^C D_t^{\alpha}u(t)=\frac{1}{\Gamma(1-\alpha)}\int_{t_0}^{t}\frac{u'(s)}{(t-s)^{\alpha}}\,ds; \]

For \(1<\alpha<2\):

\[ {}^C D_t^{\alpha}u(t)=\frac{1}{\Gamma(2-\alpha)}\int_{t_0}^{t}\frac{u''(s)}{(t-s)^{\alpha-1}}\,ds. \]

This implementation uses quadrature/Monte Carlo estimators depending on alpha and mode.

Arguments:

  • alpha: Fractional order \(\alpha\in(0,2)\).
  • time_var: Time label (must correspond to a scalar domain factor).
  • mode: "auto" selects a reasonable default; "gj" uses Gauss–Jacobi nodes; "gl" uses Gauss–Legendre; other values fall back to sampling.
  • order: Number of quadrature/sampling points.
  • tau_epsilon: Regularization for the kernel singularity near \(t=s\).
  • cluster_exponent: Optional clustering of Gauss–Legendre nodes toward \(t_0\) for the \(0<\alpha<1\) branch.

phydrax.operators.caputo_time_fractional_dw(u: DomainFunction, /, *, alpha: float, time_var: str = 't', M: int = 128, sampler: str = 'sobol_scrambled', tau_epsilon: float | str = 'auto', mode: str = 'gj') -> DomainFunction ¤

Convenience wrapper for caputo_time_fractional using order=M.

Named for a common discretization setting in physics-informed fractional models.

Continuum mechanics¤

phydrax.operators.deformation_gradient(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Deformation gradient for a displacement field.

For a displacement field \(u:\Omega\to\mathbb{R}^d\), the deformation gradient is

\[ F = I + \nabla u. \]

Arguments:

  • u: Displacement field (vector-valued).
  • var: Geometry label to differentiate with respect to.
  • mode: Autodiff mode passed to grad.
  • ad_engine: AD execution strategy passed to grad.

Returns:

  • A DomainFunction representing \(F\) with trailing shape (..., d, d).

phydrax.operators.green_lagrange_strain(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Green–Lagrange strain tensor.

With deformation gradient \(F=I+\nabla u\), the Green–Lagrange strain is

\[ E = \tfrac12\left(F^\top F - I\right). \]

Arguments:

  • u: Displacement field.
  • var: Geometry label to differentiate with respect to.
  • mode: Autodiff mode used by deformation_gradient.

Returns:

  • A DomainFunction representing \(E\) with trailing shape (..., d, d).

phydrax.operators.cauchy_strain(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Cauchy (small) strain tensor for a displacement/velocity field.

For a vector field \(u:\Omega\to\mathbb{R}^d\), the small-strain tensor is

\[ \varepsilon(u) = \tfrac12\left(\nabla u + (\nabla u)^\top\right). \]

Arguments:

  • u: Displacement/velocity field (vector-valued with trailing size \(d\)).
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to grad(u, ...).

Returns:

  • A DomainFunction representing \(\varepsilon(u)\) with trailing shape (..., d, d).

phydrax.operators.strain_rate(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Alias for cauchy_strain (common notation \(D(u)\) for velocity fields).

phydrax.operators.strain_rate_magnitude(u: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Magnitude of the strain-rate tensor.

With \(D(u)=\tfrac12(\nabla u+(\nabla u)^\top)\), this returns the scalar field

\[ \|D\| = \sqrt{2\,D{:}D} = \sqrt{2\sum_{i,j} D_{ij}^2}. \]

Arguments:

  • u: Velocity field.
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to strain_rate(u, ...).

Returns:

  • A DomainFunction representing \(\|D(u)\|\).

phydrax.operators.pk1_from_pk2(u: DomainFunction, S: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Convert 2nd PK stress to 1st PK stress.

With deformation gradient \(F\), the first Piola–Kirchhoff stress is

\[ P = F S. \]

Arguments:

  • u: Displacement field (used to compute \(F\)).
  • S: Second PK stress field.
  • var: Geometry label.
  • mode: Autodiff mode used by deformation_gradient.

Returns:

  • A DomainFunction representing \(P\) with trailing shape (..., d, d).

phydrax.operators.cauchy_from_pk2(u: DomainFunction, S: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Convert 2nd PK stress to Cauchy stress.

With deformation gradient \(F\) and \(J=\det F\), the Cauchy stress is

\[ \sigma = \frac{1}{J} F S F^\top. \]

Arguments:

  • u: Displacement field (used to compute \(F\)).
  • S: Second PK stress field.
  • var: Geometry label.
  • mode: Autodiff mode used by deformation_gradient.

Returns:

  • A DomainFunction representing \(\sigma\) with trailing shape (..., d, d).

phydrax.operators.cauchy_stress(u: DomainFunction, /, *, lambda_: DomainFunction | ArrayLike, mu: DomainFunction | ArrayLike, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Linear-elastic Cauchy stress (Hooke's law).

\[ \sigma(u) = 2\mu\,\varepsilon(u) + \lambda\,\operatorname{tr}(\varepsilon(u))\,I. \]

Here \(\varepsilon(u)=\tfrac12(\nabla u + (\nabla u)^\top)\) is the small-strain tensor, and \(\lambda,\mu\) are Lamé parameters (which may be constants or DomainFunctions).

Arguments:

  • u: Displacement field (vector-valued).
  • lambda_: First Lamé parameter \(\lambda\).
  • mu: Shear modulus \(\mu\).
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to cauchy_strain(u, ...).

Returns:

  • A DomainFunction representing \(\sigma(u)\) with trailing shape (..., d, d).

phydrax.operators.viscous_stress(u: DomainFunction, /, *, mu: DomainFunction | ArrayLike, lambda_: DomainFunction | ArrayLike | None = None, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Newtonian viscous stress tensor for a velocity field.

Uses the symmetric strain-rate tensor

\[ D(u) = \tfrac12(\nabla u + \nabla u^\top), \]

and returns

\[ \tau(u) = 2\mu\,D(u) + \lambda\,\text{tr}(D(u))\,I. \]

If lambda_ is not provided, uses Stokes' hypothesis \(\lambda=-\tfrac{2}{3}\mu\).

phydrax.operators.von_mises_stress(sigma: DomainFunction, /, *, var: str | None = None) -> DomainFunction ¤

Von Mises equivalent stress.

For deviatoric stress \(s = \sigma - \tfrac{1}{d}\text{tr}(\sigma)\,I\), returns

\[ \sigma_{\text{vm}} = \sqrt{\frac{3}{2}\,s:s}, \]

with a specialized closed form used in 2D.

Arguments:

  • sigma: Stress tensor field.
  • var: Geometry label used to infer dimension and select the 2D/3D formula.

Returns:

  • A DomainFunction representing the scalar von Mises equivalent stress.

phydrax.operators.hydrostatic_stress(sigma: DomainFunction, /, *, var: str | None = None) -> DomainFunction ¤

Hydrostatic (spherical) part of a stress tensor.

Returns

\[ \sigma_{\text{hyd}} = \frac{\text{tr}(\sigma)}{d}\,I. \]

Arguments:

  • sigma: Stress tensor field.
  • var: Geometry label used to infer \(d\).

Returns:

  • A DomainFunction representing the hydrostatic stress tensor.

phydrax.operators.hydrostatic_pressure(sigma: DomainFunction, /, *, var: str | None = None) -> DomainFunction ¤

Hydrostatic pressure from a stress tensor.

Uses the sign convention

\[ p = -\frac{1}{d}\text{tr}(\sigma). \]

Arguments:

  • sigma: Stress tensor field.
  • var: Geometry label used to infer \(d\).

Returns:

  • A DomainFunction representing the scalar pressure field \(p\).

phydrax.operators.deviatoric_stress(sigma: DomainFunction, /, *, var: str | None = None) -> DomainFunction ¤

Deviatoric part of a stress tensor.

For a stress tensor \(\sigma\), the deviatoric part is

\[ \sigma_{\text{dev}} = \sigma - \frac{\text{tr}(\sigma)}{d}\,I. \]

Arguments:

  • sigma: Stress tensor field (trailing shape (d, d)).
  • var: Geometry label used to infer \(d\).

Returns:

  • A DomainFunction representing \(\sigma_{\text{dev}}\).

phydrax.operators.maxwell_stress(*, E: DomainFunction | None = None, H: DomainFunction | None = None, epsilon: DomainFunction | ArrayLike | None = None, mu: DomainFunction | ArrayLike | None = None) -> DomainFunction ¤

Maxwell stress tensor from electric and/or magnetic fields.

For electric field \(E\) and magnetic field \(H\), the (vacuum-style) Maxwell stress tensor is

\[ T = \epsilon\left(E\otimes E - \tfrac12\|E\|_2^2 I\right) + \mu\left(H\otimes H - \tfrac12\|H\|_2^2 I\right), \]

where \(\epsilon\) and \(\mu\) can be scalars or fields.

At least one of E or H must be provided.

Arguments:

  • E: Electric field (vector-valued) or None.
  • H: Magnetic field (vector-valued) or None.
  • epsilon: Permittivity \(\epsilon\) (scalar/field); defaults to 1.
  • mu: Permeability \(\mu\) (scalar/field); defaults to 1.

Returns:

  • A DomainFunction representing the Maxwell stress tensor \(T\) with trailing shape (..., d, d).

phydrax.operators.neo_hookean_pk1(u: DomainFunction, /, *, mu: DomainFunction | ArrayLike, kappa: DomainFunction | ArrayLike, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Compressible neo-Hookean first Piola–Kirchhoff stress.

With deformation gradient \(F = I + \nabla u\) and \(J=\det F\), this returns

\[ P = \mu\,(F - F^{-T}) + \kappa\,\ln(J)\,F^{-T}, \]

where \(\mu\) is the shear modulus and \(\kappa\) is the bulk modulus.

Arguments:

  • u: Displacement field.
  • mu: Shear modulus \(\mu\) (constant or DomainFunction).
  • kappa: Bulk modulus \(\kappa\) (constant or DomainFunction).
  • var: Geometry label.
  • mode: Autodiff mode used by deformation_gradient.

Returns:

  • A DomainFunction representing the PK1 stress \(P\) with trailing shape (..., d, d).

phydrax.operators.neo_hookean_cauchy(u: DomainFunction, /, *, mu: DomainFunction | ArrayLike, kappa: DomainFunction | ArrayLike, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

Compressible neo-Hookean Cauchy stress.

With \(F = I + \nabla u\), \(J=\det F\), and left Cauchy–Green tensor \(B=FF^\top\), this returns

\[ \sigma = \frac{\mu}{J}(B - I) + \frac{\kappa\ln(J)}{J}\,I. \]

Arguments:

  • u: Displacement field.
  • mu: Shear modulus \(\mu\).
  • kappa: Bulk modulus \(\kappa\).
  • var: Geometry label.
  • mode: Autodiff mode used by deformation_gradient.

Returns:

  • A DomainFunction representing the Cauchy stress \(\sigma\) with trailing shape (..., d, d).

phydrax.operators.svk_pk2_stress(u: DomainFunction, /, *, lambda_: DomainFunction | ArrayLike, mu: DomainFunction | ArrayLike, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

St. Venant–Kirchhoff 2nd Piola–Kirchhoff stress.

With Green–Lagrange strain \(E\), the StVK constitutive law is

\[ S(E) = 2\mu\,E + \lambda\,\text{tr}(E)\,I. \]

Arguments:

  • u: Displacement field.
  • lambda_, mu: Lamé parameters (constants or DomainFunctions).
  • var: Geometry label.
  • mode: Autodiff mode used by green_lagrange_strain.

Returns:

  • A DomainFunction representing the PK2 stress \(S\) with trailing shape (..., d, d).

phydrax.operators.linear_elastic_cauchy_stress_2d(u: DomainFunction, /, *, E: DomainFunction | ArrayLike, nu: DomainFunction | ArrayLike, mode2d: Literal['plane_stress', 'plane_strain'] = 'plane_stress', var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

2D isotropic linear-elastic Cauchy stress (plane stress/strain).

Given Young's modulus \(E\) and Poisson ratio \(\nu\), returns the stress tensor \(\sigma\) computed from the small strain tensor \(\varepsilon(u)\) using either:

  • plane stress, or
  • plane strain,

controlled by mode2d.

Arguments:

  • u: Displacement field.
  • E: Young's modulus \(E\) (constant or field).
  • nu: Poisson ratio \(\nu\) (constant or field).
  • mode2d: "plane_stress" or "plane_strain".
  • var: Geometry label (must be 2D).
  • mode: Autodiff mode used by cauchy_strain.

Returns:

  • A DomainFunction representing the 2D Cauchy stress tensor with trailing shape (2, 2).

phydrax.operators.linear_elastic_orthotropic_stress_2d(u: DomainFunction, /, *, E1: DomainFunction | ArrayLike, E2: DomainFunction | ArrayLike, nu12: DomainFunction | ArrayLike, G12: DomainFunction | ArrayLike, mode2d: Literal['plane_stress', 'plane_strain'] = 'plane_stress', var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse') -> DomainFunction ¤

2D orthotropic linear-elastic Cauchy stress (plane stress/strain).

Uses orthotropic material parameters \((E_1, E_2, \nu_{12}, G_{12})\) and computes the 2D stress response under either plane stress or plane strain assumptions.

Arguments:

  • u: Displacement field.
  • E1, E2: Orthotropic Young's moduli.
  • nu12: Major Poisson ratio.
  • G12: In-plane shear modulus.
  • mode2d: "plane_stress" or "plane_strain".
  • var: Geometry label (must be 2D).
  • mode: Autodiff mode used by cauchy_strain.

Returns:

  • A DomainFunction representing the 2D orthotropic stress tensor with trailing shape (2, 2).

phydrax.operators.div_tensor(T: DomainFunction, /, *, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Divergence of a second-order tensor field.

For a tensor field \(T:\Omega\to\mathbb{R}^{d\times d}\), returns the vector field

\[ (\nabla\cdot T)_j = \sum_{i=1}^{d}\frac{\partial T_{ij}}{\partial x_i}. \]

Arguments:

  • T: Tensor field with trailing shape (d, d) (optionally with extra leading value axes).
  • var: Geometry variable label. If None, the domain must have exactly one geometry variable.
  • mode, backend, basis, periodic: Passed through to grad(T, ...).

Returns:

  • A DomainFunction representing \(\nabla\cdot T\) (vector-valued with trailing size d).

phydrax.operators.div_cauchy_stress(u: DomainFunction, /, *, lambda_: DomainFunction | ArrayLike, mu: DomainFunction | ArrayLike, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Divergence of the linear-elastic Cauchy stress.

Computes the vector field

\[ \nabla\cdot\sigma(u), \]

where \(\sigma(u)\) is given by cauchy_stress.

Arguments:

  • u, lambda_, mu: As in cauchy_stress.
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to cauchy_stress and div_tensor.

Returns:

  • A DomainFunction representing \(\nabla\cdot\sigma(u)\) (vector-valued).

phydrax.operators.navier_stokes_stress(u: DomainFunction, p: DomainFunction | ArrayLike, /, *, mu: DomainFunction | ArrayLike, lambda_: DomainFunction | ArrayLike | None = None, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Cauchy stress for incompressible/compressible Navier–Stokes.

Given a velocity field \(u\) and pressure \(p\), this returns

\[ \sigma(u,p) = \tau(u) - p\,I, \]

where \(\tau(u)\) is the (Newtonian) viscous stress returned by viscous_stress.

Arguments:

  • u: Velocity field.
  • p: Pressure (scalar field or constant).
  • mu, lambda_: Viscosity parameters for viscous_stress (with Stokes' hypothesis used when lambda_ is omitted).
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to viscous_stress.

Returns:

  • A DomainFunction representing \(\sigma(u,p)\) with trailing shape (..., d, d).

phydrax.operators.navier_stokes_divergence(u: DomainFunction, p: DomainFunction | ArrayLike, /, *, mu: DomainFunction | ArrayLike, lambda_: DomainFunction | ArrayLike | None = None, var: str | None = None, mode: Literal['reverse', 'forward'] = 'reverse', backend: Literal['ad', 'fd', 'basis'] = 'ad', basis: Literal['poly', 'fourier', 'sine', 'cosine'] = 'poly', periodic: bool = False, ad_engine: _ADEngine = 'auto') -> DomainFunction ¤

Divergence of the Navier–Stokes stress.

Computes

\[ \nabla\cdot(\tau(u) - p I), \]

where navier_stokes_stress(u, p, ...) returns the stress tensor.

Arguments:

  • u, p, mu, lambda_: As in navier_stokes_stress.
  • var: Geometry label to differentiate with respect to.
  • mode, backend, basis, periodic: Passed through to navier_stokes_stress and div_tensor.

Returns:

  • A DomainFunction representing \(\nabla\cdot(\tau(u) - p I)\) (vector-valued).