Solvers

Stokes Flow

SNES_Stokes

class underworld3.systems.solvers.SNES_Stokes[source]

Bases: SNES_Stokes_SaddlePt

Stokes equation solver for incompressible viscous flow.

This class provides functionality for a discrete representation of the Stokes flow equations assuming an incompressibility (or near-incompressibility) constraint.

The momentum equation is:

\[-\nabla \cdot \left[ \boldsymbol{\tau} - p \mathbf{I} \right] = \mathbf{f}\]

with the incompressibility constraint:

\[\nabla \cdot \mathbf{u} = 0\]

The flux term is a deviatoric stress (\(\boldsymbol{\tau}\)) related to velocity gradients (\(\nabla \mathbf{u}\)) through a viscosity tensor \(\eta\), and a volumetric (pressure) part \(p\):

\[\boldsymbol{\tau} = \frac{\eta}{2}\left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right)\]

The constraint equation gives incompressible flow by default but can be set to any function of the unknown \(\mathbf{u}\) and \(\nabla \cdot \mathbf{u}\).

Parameters:

  • mesh (uw.discretisation.Mesh) – The computational mesh.

  • velocityField (uw.discretisation.MeshVariable, optional) – Pre-existing velocity field. If None, one is created automatically.

  • pressureField (uw.discretisation.MeshVariable, optional) – Pre-existing pressure field. If None, one is created automatically.

  • degree (int, optional) – Polynomial degree for velocity interpolation. Default is 2.

  • p_continuous (bool, optional) – If True (default), pressure is continuous. Set False for discontinuous pressure.

  • verbose (bool, optional) – Enable verbose output during solving. Default is False.

  • DuDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Material derivative operator for velocity (used in derived classes).

  • DFDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Material derivative operator for flux (used in viscoelastic models).

Notes

Viscosity model: The viscosity tensor \(\boldsymbol{\eta}\) is provided by setting the constitutive_model property to one of the uw.constitutive_models classes. It may be constant, spatially varying, non-linear, or anisotropic.

Augmented Lagrangian: Setting penalty to a non-zero value adds \(\lambda \nabla \cdot \mathbf{u}\) to the weak form, improving convergence for incompressible flow (in addition to the constraint equation).

Mixed finite elements: The pressure field interpolation order determines the integration order of the mixed method and is typically lower than the velocity field order.

Viscoelastic models: For viscoelastic behaviour, the flux term contains stress history tracked on a particle swarm. See SNES_VE_Stokes.

See also

SNES_VE_Stokes

Viscoelastic Stokes solver with flux history.

SNES_NavierStokes

Navier-Stokes solver with inertial terms.

uw.constitutive_models

Available viscosity models.

Examples

>>> import underworld3 as uw
>>> mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0,0), maxCoords=(1,1), cellSize=0.1)
>>> stokes = uw.systems.Stokes(mesh, degree=2)
>>> stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel()
>>> stokes.constitutive_model.Parameters.viscosity = 1.0
>>> stokes.bodyforce = [0, -1]  # gravity
>>> stokes.solve()
instances = 0
__init__(mesh, velocityField=None, pressureField=None, degree=2, p_continuous=True, verbose=False, DuDt=None, DFDt=None)[source]
Parameters:
set_jacobian_F1_source(F1_source, linesearch='cp')[source]

Override the F1 expression used to build the Jacobian blocks.

By default, the Stokes Jacobian’s uu / up G2, G3 blocks are autodiff’d from the residual F1. Some problems benefit from differentiating a different but related expression — e.g. a smooth (softmin) viscosity formula for the Jacobian while the residual F1 keeps a sharp Min, so Newton sees a continuous derivative even when the iterate sits exactly on the yield kink.

Setting F1_source triggers a JIT recompile (the Jacobian symbols change). Pass None to revert to autodiff of F1.

Parameters:

  • F1_source (sympy.Matrix or None) – Alternative expression of the same shape as F1.sym.

  • linesearch (str or None, default "cp") – SNES linesearch type to install when F1_source is set. Defaults to "cp" (critical-point) because inexact-Newton steps don’t reliably reduce the residual norm and PETSc’s default bt (backtracking) consequently rejects useful steps with DIVERGED_LINE_SEARCH. cp accepts the predicted step at the optimum of the local linearisation and converges cleanly on the same problems where bt flails. Set to None to leave the linesearch type untouched (e.g. if you’ve already configured one via petsc_options). Has no effect when F1_source is None.

solve(zero_init_guess=True, timestep=None, _force_setup=False, verbose=False, debug=False, debug_name=None, evalf=False, order=None, picard=0, divergence_retries=0)[source]

Solve the Stokes system, with optional viscoelastic stress history.

When a constitutive model with stress history is active (DFDt is not None), the solve includes pre/post hooks for advecting stress history, updating BDF coefficients, and projecting the actual stress after the solve.

Parameters:

  • zero_init_guess (bool) – If True, use zero initial guess. Otherwise use current field values.

  • timestep (float, optional) – Advection timestep. Required when stress history is active.

  • _force_setup (bool) – Force rebuild of pointwise functions.

  • verbose (bool) – Enable verbose output.

  • evalf (bool) – Force numerical evaluation during history updates.

  • order (int, optional) – Override the VE time integration order.

  • picard (int, default=0) – Number of Picard iterations before switching to Newton. Picard uses a simplified Jacobian and can help convergence for strongly nonlinear problems like VEP at yield onset.

  • divergence_retries (int, default=0) – If SNES returns a DIVERGED reason after the main solve, re-call the underlying Newton up to this many times with a warm start (zero_init_guess=False) to try to rescue. Each retry uses the just-computed iterate plus the freshly-advected stress history, which is often enough for VEP at yield onset (Min/softmin kinks) to step off a bad Newton iterate. 0 preserves legacy behaviour (divergence is terminal). Typical useful value is 1. Only applies in the VE/VEP branch (DFDt is not None).

  • debug (bool)

  • debug_name (str)

property tau

Deviatoric stress from the most recent solve.

When stress history is active (VEP), returns psi_star[0] which contains the actual projected stress. Otherwise falls through to the base class lazy projection.

F0

Velocity equation body force term (pointwise).

The \(\mathbf{f}_0\) term represents external body forces acting on the fluid, typically gravity. This is the negative of the bodyforce parameter.

F1

Velocity equation flux/stress term (pointwise).

The \(\mathbf{F}_1\) tensor represents the stress response of the fluid, combining deviatoric stress \(\boldsymbol{\tau}\), pressure \(p\), and penalty term for weak incompressibility.

PF0

Pressure equation constraint term (continuity).

The \(h_0\) term enforces the incompressibility constraint \(\nabla \cdot \mathbf{u} = 0\). Additional constraints can be added via add_condition().

stokes_problem_description()[source]

Build residual terms for Stokes FEM assembly (deprecated).

property CM_is_setup

Whether the constitutive model is configured for this solver.

property strainrate

Symmetric strain rate tensor from velocity gradients.

The strain rate tensor \(\dot{\varepsilon}\) is computed as:

\[\dot{\varepsilon}_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\]
Returns:

Symmetric tensor of shape (dim, dim) where dim is the mesh dimensionality.

Return type:

sympy.Matrix

See also

strainrate_1d

Voigt notation (vector) form.

stress_deviator

Deviatoric stress computed from strain rate.

property strainrate_1d

Strain rate in Voigt notation (vector form).

Converts the symmetric strain rate tensor to Voigt notation for use in constitutive model calculations. In 2D, returns a 3-component vector; in 3D, a 6-component vector.

Returns:

Strain rate in Voigt notation.

Return type:

sympy.Matrix

See also

strainrate

Full tensor form.

property strainrate_star_1d

Historical strain rate in Voigt notation (for viscoelastic models).

Used in viscoelastic formulations where the stress depends on both current and historical strain rates.

Returns:

Historical strain rate in Voigt notation.

Return type:

sympy.Matrix

See also

strainrate_1d

Current strain rate.

property stress_deviator

Deviatoric stress tensor from the constitutive model.

The deviatoric stress \(\boldsymbol{\tau}\) is the traceless part of the stress tensor, computed by the constitutive model from the strain rate:

\[\boldsymbol{\tau} = \boldsymbol{\eta} : \dot{\boldsymbol{\varepsilon}}\]

For a Newtonian fluid: \(\boldsymbol{\tau} = 2\eta\dot{\boldsymbol{\varepsilon}}\)

Returns:

Deviatoric stress tensor of shape (dim, dim).

Return type:

sympy.Matrix

See also

stress

Total stress (deviatoric + pressure).

constitutive_model

Provides the viscosity relationship.

property stress_deviator_1d

Deviatoric stress in Voigt notation.

Returns:

Deviatoric stress in Voigt notation.

Return type:

sympy.Matrix

See also

stress_deviator

Full tensor form.

property stress

Total Cauchy stress tensor.

The total stress combines the deviatoric stress and pressure:

\[\boldsymbol{\sigma} = \boldsymbol{\tau} - p\mathbf{I}\]

where \(\boldsymbol{\tau}\) is the deviatoric stress and \(p\) is the pressure (positive in compression).

Returns:

Total stress tensor of shape (dim, dim).

Return type:

sympy.Matrix

See also

stress_deviator

Deviatoric (traceless) part.

property stress_1d

Total stress in Voigt notation.

Returns:

Total stress in Voigt notation.

Return type:

sympy.Matrix

See also

stress

Full tensor form.

property div_u

Velocity divergence.

For incompressible flow, this should be zero:

\[\nabla \cdot \mathbf{u} = 0\]
Returns:

Scalar divergence expression.

Return type:

sympy.Expr

Notes

Non-zero divergence indicates compressibility or mass sources/sinks.

property constraints

Constraint equation for the saddle-point system.

By default, this is the incompressibility constraint \(\nabla \cdot \mathbf{u} = 0\). Can be modified for compressible or other constrained formulations.

Returns:

Constraint expression.

Return type:

sympy.Expr

property bodyforce

Body force vector (source term).

The volumetric body force \(\mathbf{f}\) appears on the right-hand side of the momentum equation:

\[-\nabla \cdot \boldsymbol{\sigma} = \mathbf{f}\]

Common examples include gravity (\(\rho\mathbf{g}\)) or buoyancy forces.

Returns:

Body force vector expression.

Return type:

UWexpression

property saddle_preconditioner

Preconditioner for the Schur complement in the saddle-point system.

For the Stokes system, the default preconditioner is \(1/\eta\) (inverse viscosity), which approximates the Schur complement \(\mathbf{S} \approx \mathbf{B}\mathbf{A}^{-1}\mathbf{B}^T\).

Returns:

Preconditioner expression (typically inverse viscosity).

Return type:

sympy.Expr

Notes

A good preconditioner significantly improves convergence of the iterative solver. For variable viscosity, use the local viscosity.

property penalty

Augmented Lagrangian penalty parameter.

The penalty \(\lambda\) adds a term to the weak form that penalizes non-zero divergence:

\[\lambda \int (\nabla \cdot \mathbf{u})(\nabla \cdot \mathbf{v}) \, dV\]

This improves convergence for incompressible flow without changing the solution (since \(\nabla \cdot \mathbf{u} = 0\) at convergence).

Returns:

Augmented Lagrangian penalty parameter (typically \(O(1)\)).

Return type:

UWexpression

Notes

Set to zero for standard Stokes without augmentation. Unlike classical penalty methods that require very large values, the Augmented Lagrangian approach uses modest penalties of \(O(1)\) to improve solver convergence.

References

Glowinski, R., & Le Tallec, P. (1989). Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics. SIAM. https://doi.org/10.1137/1.9781611970838

estimate_dt()[source]

Calculates an appropriate advective timestep for the Stokes solver.

The Stokes equations are quasi-static (no time derivative ∂v/∂t), so there is no diffusive CFL constraint. The only relevant timescale is the advective one: how long it takes material to cross an element.

This method computes a per-element timestep:

\[\delta t_i = h_i / \|\mathbf{v}_i\|\]

where \(h_i\) is the element radius and \(\mathbf{v}_i\) is the velocity at the element centroid, then returns the global minimum. This is more accurate than using global max velocity with global min element size, especially for non-uniform meshes with spatially varying velocity.

Returns:

Pint Quantity or float: The advective timestep with physical time units if a model with reference scales is available, otherwise nondimensional.

SNES_VE_Stokes

Viscoelastic extension of the Stokes solver.

class underworld3.systems.solvers.SNES_VE_Stokes[source]

Bases: SNES_Stokes

Viscoelastic Stokes solver (backward-compatibility wrapper).

Deprecated since version Use: uw.systems.Stokes directly with a ViscoElasticPlasticFlowModel constitutive model. The Stokes solver now creates stress history infrastructure automatically when the constitutive model is assigned (the lazy-creation pathway also reads stress_history_ddt_kwargs from the model — required for integrator='etd' to allocate forcing_star). VE_Stokes pre-creates the DDt at solver __init__ time, before the model exists, so it can’t see those kwargs and is incompatible with integrator='etd'.

Migration: replace:

stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=2)
stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel

with:

stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(
    stokes.Unknowns, order=2,
)
Parameters:

  • mesh (Mesh) – The computational mesh.

  • order (int, default=2) – BDF time integration order for stress history.

  • **kwargs – All other arguments are passed to SNES_Stokes.

instances = 0
__init__(mesh, velocityField=None, pressureField=None, degree=2, order=2, p_continuous=True, verbose=False, DuDt=None, DFDt=None)[source]
Parameters:
property delta_t

Elastic timestep from the constitutive model.

Scalar Equations

SNES_Poisson

class underworld3.systems.solvers.SNES_Poisson[source]

Bases: SNES_Scalar

Poisson equation solver.

Provides a discrete representation of the Poisson equation:

\[\nabla \cdot \left[ \boldsymbol{\kappa} \nabla u \right] = f\]

where \(\mathbf{F} = \boldsymbol{\kappa} \nabla u\) relates the flux to gradients in the unknown \(u\).

Parameters:

  • mesh (Mesh) – The computational mesh.

  • u_Field (MeshVariable, optional) – Pre-existing mesh variable for the solution. If None, one is created.

  • verbose (bool, optional) – Enable verbose output during solve.

  • degree (int, optional) – Polynomial degree for the solution field (default: 2).

  • DuDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Time derivative operator for time-dependent problems.

  • DFDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Time derivative operator for the flux.

__init__(mesh, u_Field=None, verbose=False, degree=2, DuDt=None, DFDt=None)[source]
Parameters:
F0

Source term for the Poisson equation (pointwise).

The \(f_0\) term represents the source/sink in the diffusion equation. Set via the f property (e.g., poisson.f = heat_source).

F1

Diffusive flux term for the Poisson equation (pointwise).

The \(\mathbf{F}_1\) vector represents the flux \(k \nabla u\) from the constitutive model, where \(k\) is the diffusivity.

poisson_problem_description()[source]

Build residual terms for Poisson FEM assembly.

property f

Source term for the Poisson equation.

The source term \(f\) appears on the right-hand side:

\[\nabla \cdot (\kappa \nabla u) = f\]
Returns:

Source term expression (scalar, shape (1, 1)).

Return type:

sympy.Matrix

See also

constitutive_model

Provides the diffusivity \(\kappa\).

property CM_is_setup

Whether the constitutive model is configured for this solver.

SNES_Darcy

class underworld3.systems.solvers.SNES_Darcy[source]

Bases: SNES_Scalar

Darcy flow equation solver for groundwater problems.

Provides a discrete representation of the groundwater flow equations:

\[\underbrace{S_s \frac{\partial h}{\partial t}}_{\dot{u}} - \nabla \cdot \underbrace{\left[ \boldsymbol{\kappa} \nabla h - \boldsymbol{s} \right]}_{\mathbf{F}} = \underbrace{W}_{h}\]

The flux term \(\mathbf{F}\) relates the effective velocity to pressure gradients:

\[\boldsymbol{v} = \boldsymbol{\kappa} \nabla h - \boldsymbol{s}\]
Parameters:

  • mesh (Mesh) – The computational mesh.

  • h_Field (MeshVariable, optional) – Mesh variable for hydraulic head. Created automatically if not provided.

  • v_Field (MeshVariable, optional) – Mesh variable for Darcy velocity. Created automatically if not provided.

  • degree (int, default=2) – Polynomial degree for the finite element discretization.

  • verbose (bool, default=False) – Enable verbose output.

  • DuDt (optional) – Time derivative operator for the unknown.

  • DFDt (optional) – Time derivative operator for the flux.

Notes

  • The unknown is \(h\), the hydraulic head

  • The permeability tensor \(\kappa\) is set via the constitutive_model property using one of the uw.constitutive_models classes

  • \(W\) is a pressure source term

  • \(S_s\) is the specific storage coefficient

  • The time-dependent term \(\dot{f}\) is not implemented in this version

  • The solver returns both the primary field and the Darcy flux (mean-flow velocity)

See also

SNES_Poisson

Related diffusion-only solver.

uw.constitutive_models.DarcyFlowModel

Constitutive model for Darcy flow.

__init__(mesh, h_Field=None, v_Field=None, degree=2, verbose=False, DuDt=None, DFDt=None)[source]
Parameters:

  • mesh (Mesh)

  • h_Field (EnhancedMeshVariable | None)

  • v_Field (EnhancedMeshVariable | None)

  • degree (int)

F0

Source term for the Darcy pressure equation (pointwise).

The \(f_0\) term represents fluid sources/sinks. Set via the f property.

F1

Darcy flux term (pointwise).

The \(\mathbf{F}_1\) vector is the Darcy flux \(K \nabla p\), where \(K\) is the permeability.

darcy_problem_description()[source]

Build residual terms for Darcy flow FEM assembly.

property f

Source term for the Darcy equation.

property darcy_flux

Darcy flux velocity computed from pressure gradient.

property v

Projected Darcy velocity field.

solve(zero_init_guess=True, timestep=None, verbose=False, _force_setup=False, divergence_retries=0)[source]

Solve the Darcy flow system.

Computes the pressure field and Darcy flux velocity.

Parameters:

  • zero_init_guess (bool, optional) – If True (default), start from zero initial guess. If False, use current field values as initial guess.

  • timestep (float, optional) – Timestep value for inertial terms (if applicable).

  • verbose (bool, optional) – If True, print solver progress information.

  • _force_setup (bool, optional) – Force re-setup of solver even if already configured.

  • divergence_retries (int, optional) – If SNES reports DIVERGED, retry with warm start up to this many times. 0 preserves legacy behaviour.

Notes

After solving, the pressure field self.u and velocity field self.v contain the solution.

SNES_AdvectionDiffusion

class underworld3.systems.solvers.SNES_AdvectionDiffusion[source]

Bases: SNES_Scalar

Advection-diffusion equation solver using semi-Lagrangian Crank-Nicolson.

Implements the characteristics-based method described in Spiegelman & Katz (2006):

\[\underbrace{\frac{\partial u}{\partial t} + \left( \mathbf{v} \cdot \nabla \right) u}_{\dot{u}} - \nabla \cdot \underbrace{\left[ \boldsymbol{\kappa} \nabla u \right]}_{\mathbf{F}} = \underbrace{f}_{\mathbf{h}}\]

The flux term \(\mathbf{F}\) relates diffusive fluxes to gradients in \(u\). Advective fluxes along the velocity field \(\mathbf{v}\) are handled in the \(\dot{u}\) term.

The time derivative \(\dot{u}\) involves upstream sampling to find \(u^*\), the value of \(u\) at points which later arrive at mesh nodes. This is achieved using a hidden swarm variable advected backwards from nodal points automatically during solve.

Parameters:

  • mesh (Mesh) – The computational mesh.

  • u_Field (MeshVariable) – Mesh variable for the transported scalar.

  • V_fn (MeshVariable or sympy.Basic) – Velocity field for advection.

  • order (int, default=1) – Time integration order (1 or 2).

  • restore_points_func (callable, optional) – Function to restore particles to valid domain.

  • verbose (bool, default=False) – Enable verbose output.

  • DuDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Time derivative operator for the unknown.

  • DFDt (SemiLagrangian_DDt or Lagrangian_DDt, optional) – Time derivative operator for the flux.

  • monotone_mode (str or None, optional) –

    Monotonicity limiter for the semi-Lagrangian trace-back. Forwarded to the internally-constructed SemiLagrangian_DDt instances for DuDt and DFDt.

    • None (default): pure FE trace-back. Can overshoot at non-nodal upstream points in cells with sharp gradients (e.g. thin boundary layers in deformed cells); legacy behaviour.

    • "clamp": clip the FE trace-back result to [nbr_min, nbr_max] of the k = dim + 1 nearest psi_star DOFs. Bit-identical to pure FE in smooth regions; bounds overshoots where the FE interpolant would otherwise leave the local data range.

    • "pick": keep the FE result where in-bounds; for out-of-bounds DOFs only, re-evaluate via Shepard’s-method RBF. More conservative than clamp at the catastrophe edge.

    When a user supplies a pre-built DuDt, this kwarg is applied to the internally-constructed DFDt only — the user’s DuDt already encodes whatever monotone_mode it was constructed with.

  • theta (float, default=0.5) –

    Adams-Moulton theta for the diffusive flux at order 1. Forwarded to the internally-constructed SemiLagrangian_DDt instances (same forwarding rule as monotone_mode).

    • 0.5 (default): Crank-Nicolson, A-stable but not L-stable. Second-order accurate. Rings on stiff modes with sharp gradients (negative amplification factor).

    • 1.0: Backward Euler, L-stable, monotone for the diffusive flux. First-order accurate. Recommended when SLCN+CN ringing dominates the discretisation error.

    • 0.0: Forward Euler — unstable for stiff diffusion; included for completeness.

Notes

  • The diffusivity \(\kappa\) is set via the constitutive_model property

  • Sources \(f\) can be any sympy expression involving mesh/swarm variables

References

Spiegelman, M., & Katz, R. F. (2006). A semi-Lagrangian Crank-Nicolson algorithm for the numerical solution of advection-diffusion problems. Geochemistry, Geophysics, Geosystems, 7(4). https://doi.org/10.1029/2005GC001073

See also

SNES_Diffusion

Pure diffusion solver without advection.

SNES_Navier_Stokes

Full momentum advection-diffusion.

__init__(mesh, u_Field, V_fn, order=1, restore_points_func=None, verbose=False, DuDt=None, DFDt=None, monotone_mode=None, theta=0.5)[source]
Parameters:
property F0

Pointwise source term including time derivative.

property F1

Pointwise diffusive flux term (Adams-Moulton integration).

adv_diff_slcn_problem_description()[source]

Build residual terms for advection-diffusion FEM assembly.

property f

Source term for the advection-diffusion equation.

The source \(f\) appears on the right-hand side:

\[\frac{\partial u}{\partial t} + \mathbf{V} \cdot \nabla u = \nabla \cdot (\kappa \nabla u) + f\]
Returns:

Source term expression.

Return type:

sympy.Matrix

property V_fn

Velocity field for advection.

The advection velocity \(\mathbf{V}\) transports the scalar field \(u\). Can be a MeshVariable or symbolic expression.

Returns:

Velocity vector expression.

Return type:

sympy.Matrix

property delta_t

Timestep for time integration.

The timestep \(\Delta t\) controls the temporal discretization. For explicit advection, this should satisfy the CFL condition:

\[\Delta t < \frac{h}{|\mathbf{V}|}\]

where \(h\) is the element size and \(|\mathbf{V}|\) is the velocity magnitude.

Returns:

Timestep value.

Return type:

UWexpression

See also

estimate_dt

Computes a stable timestep automatically.

estimate_dt(direction_aware=False)[source]

Estimate an appropriate timestep for the advection-diffusion solver.

This is an implicit solver so the returned \(\delta t\) is the minimum of:

  • \(\delta t_{\textrm{diff}}\): typical time for diffusion across an element

  • \(\delta t_{\textrm{adv}}\): typical element-crossing time for a fluid parcel

Parameters:

direction_aware (bool, default False) – If True, the advective dt uses the per-cell extent along the local velocity direction — h_eff_c = max_i(s_i) - min_i(s_i) where s_i = (x_i - centroid) · v̂ over the cell vertices. This is the distance material actually traverses through the cell per unit |v|, and is always ≥ the isotropic mesh._radii estimate, by 1.5–3× for equant cells (geometric factor) and up to ~10× for cells that the mover has stretched along the flow direction. On adapted meshes the gain is substantial; on uniform meshes it’s the geometric factor only. Off by default to preserve historical behaviour; safe to enable everywhere once validated.

Returns:

The recommended timestep with physical time units if a model with reference scales is available, otherwise nondimensional.

Return type:

pint.Quantity or float

solve(zero_init_guess=True, timestep=None, _force_setup=False, _evalf=False, verbose=False, divergence_retries=0)[source]

Generates solution to constructed system.

Params

zero_init_guess:

If True, a zero initial guess will be used for the system solution. Otherwise, the current values of self.u will be used.

divergence_retries:

If SNES reports DIVERGED, retry with warm start up to this many times. 0 preserves legacy behaviour.

Parameters:

  • zero_init_guess (bool)

  • timestep (float)

  • _force_setup (bool)

  • divergence_retries (int)

SNES_Diffusion

class underworld3.systems.solvers.SNES_Diffusion[source]

Bases: SNES_Scalar

Diffusion equation solver using mesh-based finite elements.

Solves the scalar diffusion equation:

\[\underbrace{\frac{\partial u}{\partial t}}_{\dot{f}} - \nabla \cdot \underbrace{\left[ \boldsymbol{\kappa} \nabla u \right]}_{\mathbf{F}} = \underbrace{f}_{h}\]

The flux term \(\mathbf{F}\) relates diffusive fluxes to gradients in the unknown \(u\).

Parameters:

  • mesh (Mesh) – The computational mesh.

  • u_Field (MeshVariable) – Mesh variable for the diffusing scalar.

  • order (int, default=1) – Time integration order.

  • theta (float, default=0.0) – Time integration parameter (0=explicit, 0.5=Crank-Nicolson, 1=implicit).

  • evalf (bool, default=False) – Numerically evaluate symbolic expressions during setup.

  • verbose (bool, default=False) – Enable verbose output.

  • DuDt (Eulerian_DDt, SemiLagrangian_DDt, or Lagrangian_DDt, optional) – Time derivative operator for the unknown.

  • DFDt (Eulerian_DDt, SemiLagrangian_DDt, or Lagrangian_DDt, optional) – Time derivative operator for the flux.

Notes

  • The diffusivity \(\kappa\) is set via the constitutive_model property

  • Sources \(f\) can be any sympy expression involving mesh/swarm variables

See also

SNES_AdvectionDiffusion

Adds advection transport.

SNES_Poisson

Steady-state diffusion (no time derivative).

__init__(mesh, u_Field, order=1, theta=0.0, evalf=False, verbose=False, DuDt=None, DFDt=None)[source]
Parameters:
property F0

Pointwise source term including time derivative.

property F1

Pointwise diffusive flux term.

property f

Source term for the diffusion equation.

property delta_t

Timestep for time integration.

estimate_dt()[source]

Estimate an appropriate timestep for the diffusion solver.

This solver only has a diffusive component, so the returned \(\delta t\) is:

  • \(\delta t_{\textrm{diff}}\): typical time for diffusion across an element

Returns:

The diffusive timestep with physical time units if a model with reference scales is available, otherwise nondimensional.

Return type:

pint.Quantity or float

solve(zero_init_guess=True, timestep=None, evalf=False, _force_setup=False, verbose=False, divergence_retries=0)[source]

Generates solution to constructed system.

Params

zero_init_guess:

If True, a zero initial guess will be used for the system solution. Otherwise, the current values of self.u will be used.

divergence_retries:

If SNES reports DIVERGED, retry with warm start up to this many times. 0 preserves legacy behaviour.

Parameters:

Projection Solvers

SNES_Projection

class underworld3.systems.solvers.SNES_Projection[source]

Bases: SNES_Scalar

Scalar projection solver for mapping functions to mesh variables.

Solves \(u = \tilde{f}\) where \(\tilde{f}\) is a function that can be evaluated within an element and \(u\) is a mesh variable with associated shape functions.

Typically used to obtain a continuous representation of a function not well-defined at mesh nodes (e.g., derivatives or flux components). More broadly, it is a projection from one basis to another.

Strong form (the screened-Poisson smoother)

The projection is implemented by solving

\[u - \nabla \cdot \left( \alpha \nabla u \right) = \tilde{f},\]

or equivalently, in the more familiar Helmholtz form

\[u - \alpha \, \nabla^{2} u \;=\; \tilde{f}.\]

With \(\alpha = 0\) this is a pure pointwise L2 projection of \(\tilde f\) onto the discrete space of \(u\). With \(\alpha > 0\) it is a screened-Poisson smoother: the equation enforces a balance between fidelity to \(\tilde f\) and curvature of \(u\). The natural length scale that emerges from this balance is

\[L \;=\; \sqrt{\alpha},\]

so \(\alpha\) has dimensions of length squared. The free-space Green’s function of the operator decays as \(\exp(-r/L)\) (in 2D it is \(G(r) \propto K_{0}(r/L)/L^{2}\)), so the solution behaves like a Gaussian-like convolution of \(\tilde f\) of width \(L\) — obtained implicitly by one elliptic solve, without ever assembling the kernel. Features in \(\tilde f\) of scale much smaller than \(L\) are attenuated; features much larger than \(L\) pass through essentially unchanged.

Weak form

Multiplying by a test function \(v\) and integrating by parts gives the symmetric weak form actually assembled by PETSc:

\[\int_\Omega (u - \tilde f)\, v \; + \; \int_\Omega \alpha \, \nabla u \cdot \nabla v \;=\; 0,\]

which is exactly minimising \(\tfrac{1}{2}\!\int (u-\tilde f)^2 + \tfrac{\alpha}{2}\!\int |\nabla u|^2\) — a Tikhonov-regularised L2 projection.

Setting the smoothing length

Two equivalent accessors are provided:

  • smoothing — set the coefficient \(\alpha\) directly (units of length²). Historically used with tiny values (e.g. \(10^{-6}\)) as a numerical regulariser, which corresponds to a sub-grid \(L\) and produces no physical smoothing.

  • smoothing_length — set \(L\) directly (length units, unit-aware via underworld3.non_dimensionalise()). This is the recommended path when you actually want the projection to act as a low-pass filter of a chosen physical scale.

See also

SNES_Vector_Projection

Vector field projection.

SNES_Tensor_Projection

Tensor field projection.

param mesh:

The computational mesh.

type mesh:

Mesh

param u_Field:

Target mesh variable for the projection.

type u_Field:

MeshVariable, optional

param scalar_Field:

Alternative name for the target field.

type scalar_Field:

MeshVariable, optional

param degree:

Polynomial degree for the finite element space.

type degree:

int, default=2

param solver_name:

Name for the solver instance.

type solver_name:

str, optional

param verbose:

Enable verbose output.

type verbose:

bool, default=False

__init__(mesh, u_Field=None, degree=2, verbose=False)[source]
Parameters:

  • mesh (Mesh)

  • u_Field (EnhancedMeshVariable)

F0

Projection misfit term (pointwise).

The \(f_0\) term measures the weighted difference between the mesh variable and the target function for L2 projection.

F1

Projection smoothing term (pointwise).

The \(\mathbf{F}_1\) term provides optional regularization via \(\alpha \nabla u\), where \(\alpha\) is the smoothing parameter.

uw_function

Function to project onto mesh

property smoothing

Smoothing coefficient \(\alpha\) of the screened-Poisson form.

The projection solves \(u - \nabla\!\cdot\!(\alpha\,\nabla u) = \tilde f\), so \(\alpha\) has dimensions of length² and the implied smoothing length is \(L = \sqrt{\alpha}\). The free-space Green’s function decays as \(\exp(-r/L)\), giving the projection the action of a Gaussian-like convolution of width \(L\) — without ever forming the kernel.

See the class docstring for the full derivation.

Two usage patterns

  • Pure L2 projection — smoothing = 0 (or omit). No regularisation; u is the best L2 fit to ũ.

  • Genuine length-scale smoother — set smoothing_length to the desired physical length (unit-aware), or equivalently smoothing = L**2. The output is ũ smoothed at scale L.

Note

Historical UW3 code occasionally sets smoothing to a tiny number (e.g. 1e-6) as a numerical regulariser against rank deficiency. That value corresponds to \(L \approx 10^{-3}\) in the problem’s ND length units — almost always well below the cell size, so it does no useful smoothing. Use a true sub-grid value only if you need that regularisation; for filtering, use smoothing_length.

property smoothing_length

Smoothing length scale \(L\) (length units, unit-aware).

L-valued view on smoothing with the convention \(L = \sqrt{\alpha}\). Setting smoothing_length = L is equivalent to setting smoothing = L**2, but is the natural physical knob because \(L\) is what the smoother actually filters by.

Mathematical meaning

The projection then solves the screened-Poisson equation

\[u \;-\; L^{2}\,\nabla^{2} u \;=\; \tilde f,\]

whose Green’s function decays as \(\exp(-r/L)\). In practice this acts like a Gaussian convolution of width \(L\):

  • features in \(\tilde f\) with spatial scale \(\ll L\) are damped (roughly as \(1/(1+k^{2}L^{2})\) for a wavenumber \(k\));

  • features with scale \(\gg L\) are preserved;

  • features at \(\sim L\) are attenuated by a factor near 1/2.

Choosing \(L\) smaller than the local mesh size has essentially no effect (the field is already band-limited by the discretisation). A useful default for light de-noising is \(L \approx 1\!-\!2\,h\), where \(h\) is a representative cell size.

Units

The setter accepts a plain number (assumed already non-dimensional), a pint Quantity with length units, or any unit-aware object understood by underworld3.non_dimensionalise(). Internally the squared non-dimensional value is stored in self._smoothing (so smoothing and smoothing_length stay consistent).

The getter returns a Pint Quantity with length units when a scaling context is configured; otherwise the plain non-dimensional float \(\sqrt{\alpha}\).

property uw_weighting_function

Weighting function applied during projection.

SNES_Vector_Projection

class underworld3.systems.solvers.SNES_Vector_Projection[source]

Bases: SNES_Vector

Vector projection solver for mapping vector functions to mesh variables.

Solves \(\mathbf{u} = \tilde{\mathbf{f}}\) where \(\tilde{\mathbf{f}}\) is a vector function that can be evaluated within an element and \(\mathbf{u}\) is a vector mesh variable with associated shape functions.

Typically used to obtain a continuous representation of a vector function not well-defined at mesh nodes (e.g., gradient or flux vectors), or as a length-scale-aware smoother of an existing vector field.

Strong form (screened-Poisson, vector-valued)

\[\mathbf{u} \;-\; \nabla \cdot \left( \alpha\, \nabla \mathbf{u} \right) \;+\; \lambda \left( \nabla \cdot \mathbf{u} \right) \mathbf{I} \;=\; \tilde{\mathbf{f}} .\]

The \(\alpha\)-term is the same screened-Poisson smoother as in SNES_Projection, applied component-wise: it has the same \(L = \sqrt{\alpha}\) smoothing-length interpretation and the same \(\exp(-r/L)\) Green’s function. The extra \(\lambda\) term is a divergence penalty (see penalty) — set it nonzero when you want an approximately solenoidal projection of \(\tilde{\mathbf{f}}\).

Setting \(\alpha = 0\) (and \(\lambda = 0\)) gives a pure pointwise L2 projection. See SNES_Projection for the full mathematical context; the relationship between smoothing (units length²) and smoothing_length (units length) is the same as for the scalar projection.

param mesh:

The computational mesh.

type mesh:

Mesh

param u_Field:

Target vector mesh variable for the projection.

type u_Field:

MeshVariable, optional

param degree:

Polynomial degree for the finite element space.

type degree:

int, default=2

param verbose:

Enable verbose output.

type verbose:

bool, default=False

See also

SNES_Projection

Scalar field projection (full mathematical detail).

SNES_Tensor_Projection

Tensor field projection.

__init__(mesh, u_Field=None, degree=2, verbose=False)[source]
Parameters:

  • mesh (Mesh)

  • u_Field (EnhancedMeshVariable)

F0

f_0(u)

Type:

Vector projection pointwise misfit term

F1

F_1(u)

Type:

Vector projection pointwise smoothing term

projection_problem_description()[source]

Build residual terms for vector projection FEM assembly.

uw_function

Vector function to project onto mesh

property smoothing

Smoothing coefficient \(\alpha\) (units length²).

Coefficient of the \(\nabla\!\cdot\!(\alpha\,\nabla \mathbf u)\) term in the vector screened-Poisson equation (see the class docstring). Acts component-wise; the smoothing length is \(L = \sqrt{\alpha}\) and the Green’s function decays as \(\exp(-r/L)\).

Use smoothing_length for the L-valued, unit-aware knob. See SNES_Projection.smoothing for the full derivation and usage patterns.

property smoothing_length

Smoothing length \(L\) (length units, unit-aware).

L-valued view on smoothing, with \(L = \sqrt{\alpha}\). The projection then acts as a component-wise screened-Poisson smoother \(\mathbf u - L^{2}\,\nabla^{2}\mathbf u = \tilde{\mathbf f}\), i.e. a Gaussian-like convolution of width \(L\) applied to each Cartesian component of the input vector field.

Choose \(L \gtrsim h\) (a cell size) for noticeable smoothing; \(L < h\) does essentially nothing because the discretisation already band-limits at that scale. See SNES_Projection.smoothing_length for the full mathematical and units discussion.

property penalty

Divergence penalty \(\lambda\) for (approx.) incompressibility.

Coefficient of the \(\lambda (\nabla\!\cdot\!\mathbf u) \mathbf I\) term in the vector projection. Large positive values bias the projection toward \(\nabla\!\cdot\!\mathbf u = 0\), i.e. a solenoidal approximation of \(\tilde{\mathbf f}\). Has no length interpretation — unlike smoothing it does not introduce a filter scale.

property uw_weighting_function

Weighting function applied during projection.