underworld.systems module¶
This module contains routines relating to differential system.
Functions¶
underworld.systems.Solver 
This method simply returns a necessary solver for the provided system. 
Classes¶
underworld.systems.AdvectionDiffusion 
This class provides functionality for a discrete representation of an advectiondiffusion equation. 
underworld.systems.HeatSolver 
Steady State Heat Equation Solver. 
underworld.systems.SteadyStateDarcyFlow 
This class provides functionality for a discrete representation of the steady state darcy flow equation. 
underworld.systems.SteadyStateHeat 
This class provides functionality for a discrete representation of the steady state heat equation. 
underworld.systems.Stokes 
This class provides functionality for a discrete representation of the Stokes flow equations. 
underworld.systems.StokesSolver 
The Block Stokes Schur Complement Solver: This solves the saddlepoint system 
underworld.systems.SwarmAdvector 
Objects of this class advect a swarm through time using the provided velocity field. 
underworld.systems.TimeIntegration 
Abstract class for integrating numerical objects (fields, swarms, etc.) in time. 

underworld.systems.
Solver
(eqs, type='BSSCR', *args, **kwargs)[source]¶ This method simply returns a necessary solver for the provided system.

class
underworld.systems.
AdvectionDiffusion
(phiField=None, velocityField=None, fn_diffusivity=None, fn_sourceTerm=None, method='SUPG', conditions=[], phiDotField=None, allow_non_q1=False, gauss_swarm=None, **kwargs)[source]¶ Bases:
object
This class provides functionality for a discrete representation of an advectiondiffusion equation.
\[\frac{\partial\phi}{\partial t} + {\bf u } \cdot \nabla \phi= \nabla { ( k \nabla \phi ) } + H\]Two methods are available to integrate the scalar \(\phi\) through time:
 SUPG  The Streamline Upwind Petrov Galerkin method. [1]
 SLCN  The SemiLagrangian CrankNicholson method. [2]
SLCN is the preferred method for Q1 elements on a orthogonal cartesian meshes. It is quicker, less diffusive and unconditionally stable. SUPG, the legacy method, is more robust for arbitrarily deformed meshes. Both methods are considered EXPERIMENTAL for non Q1 element meshes.
Parameters:  phiField (underworld.mesh.MeshVariable) – The concentration field, typically the temperature field
 velocityField (underworld.mesh.MeshVariable) – The velocity field.
 fn_diffusivity (underworld.function.Function) – A function that defines the diffusivity within the domain.
 fn_sourceTerm (underworld.function.Function) – A function that defines the heating within the domain. Optional.
 conditions (underworld.conditions.SystemCondition) – Numerical conditions to impose on the system. This should be supplied as the condition itself, or a list object containing the conditions.
 phiDotField (underworld.mesh.MeshVariable) – Only used for SUPG. A MeshVariable that defines the initial time derivative of the phiField. Typically 0 at the beginning of a model, e.g. phiDotField.data[:]=0 When using a phiField loaded from disk one should also load the phiDotField to ensure the solving method has the time derivative information for a smooth restart. No Dirichlet conditions are required for this field as the phiField degrees of freedom map exactly to this field’s Dirichlet conditions, the value of which ought to be 0 for constant values of phi.
 gauss_swarm (underworld.swarm.GaussIntegrationSwarm) – If provided this gauss_swarm will be used for (gaussian) numerical integration rather than a default gauss integration swarm that is automatically generated and dependent on the element order of the mesh.
 allow_non_q1 (Bool (default False)) – Allow the integration to perform over a non Q1 element mesh. (Under Q2 elements instabilities have been observed as the implementation is only for Q1 elements)
Notes
Constructor must be called by collectively all processes.
[1] Brooks, A. N. and Hughes, T. J. R., “Streamline Upwind/PetrovGalerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible NavierStokes Equations”, Comput. Methods Appl. Mech. Eng., Aug 1990, 199259. [2] Spiegelman, M. and Katz, R.F., “A semiLagrangian CrankNicolson algorithm for the numerical solution of advectiondiffusion problems”, Geochemistry, Geophysics, Geosystems, 7(4), 2006.

class
underworld.systems.
HeatSolver
(heatSLE, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
Steady State Heat Equation Solver.

configure
(solve_type='')[source]¶ Configure solver.
solve_type can be one of:
 mumps : MUMPS parallel direct solver.
 superludist : SuperLU parallel direct solver.
 superlu : SuperLU direct solver (serial only).
 lu : LU direct solver (serial only).

solve
(nonLinearIterate=None, nonLinearTolerance=0.01, nonLinearMaxIterations=500, callback_post_solve=None, **kwargs)[source]¶ Solve the HeatEq system
Parameters:  nonLinearIterate (bool) – True will perform non linear iterations iterations, False (or 0) will not
 nonLinearTolerance (float, Default=1.0e2) – Relative tolerance criterion for the change in the velocity field
 nonLinearMaxIterations (int, Default=500) – Maximum number of non linear iteration to perform
 callback_post_solve (func, Default=None) – Optional callback function to be performed at the end of a linear solve iteration. Commonly this will be used to perform operations between non linear iterations, for example, calibrating the solution or removing the system null space.


class
underworld.systems.
SteadyStateDarcyFlow
(pressureField, fn_diffusivity, fn_bodyforce=None, voronoi_swarm=None, conditions=[], velocityField=None, swarmVarVelocity=None, _removeBCs=True, gauss_swarm=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
This class provides functionality for a discrete representation of the steady state darcy flow equation.
The class uses a standard Galerkin finite element method to construct a system of linear equations which may then be solved using an object of the underworld.system.Solver class.
The underlying element types are determined by the supporting mesh used for the ‘pressureField’.
The strong form of the given boundary value problem, for \(f\), \(q\) and \(h\) given, is
\[\begin{split}\begin{align} q_i =& \kappa \, ( u_{,i} + S_i ) & \\ q_{i,i} =& \: f & \text{ in } \Omega \\ u =& \: q & \text{ on } \Gamma_q \\ q_i n_i =& \: h & \text{ on } \Gamma_h \\ \end{align}\end{split}\]where,
 \(\kappa\) is the diffusivity, \(u\) is the pressure,
 \(S\) is a flow bodysource, for example due to gravity,
 \(f\) is a source term, \(q\) is the Dirichlet condition, and
 \(h\) is a Neumann condition.
The problem boundary, \(\Gamma\), admits the decomposition \(\Gamma=\Gamma_q\cup\Gamma_h\) where \(\emptyset=\Gamma_q\cap\Gamma_h\). The equivalent weak form is:
\[\int_{\Omega} w_{,i} \, q_i \, d \Omega = \int_{\Omega} w \, f \, d\Omega + \int_{\Gamma_h} w \, h \, d \Gamma\]where we must find \(u\) which satisfies the above for all \(w\) in some variational space.
Parameters:  pressureField (underworld.mesh.MeshVariable) – The solution field for pressure.
 fn_diffusivity (underworld.function.Function) – The function that defines the diffusivity across the domain.
 fn_bodyforce (underworld.function.Function) – A function that defines the flow bodyforce across the domain, for example gravity. Must be a vector. Optional.
 voronoi_swarm (underworld.swarm.Swarm) – A swarm with just one particle within each cell should be provided. This avoids the evaluation of the velocity on nodes and inaccuracies arising from diffusivity changes within cells. If a swarm is provided, voronoi type numerical integration is utilised. The provided swarm is used as the basis for the voronoi integration. If no voronoi_swarm is provided, Gauss integration is used.
 conditions (underworld.conditions.SystemCondition) – Numerical conditions to impose on the system. This should be supplied as the condition itself, or a list object containing the conditions.
 gauss_swarm (underworld.swarm.GaussIntegrationSwarm) – If provided this gauss_swarm will be used for (gaussian) numerical integration rather than a default gauss integration swarm that is automatically generated and dependent on the element order of the mesh. NB: if a voronoi_swarm is defined it OVERRIDES this gauss_swarm as the preferred integration swarm (quadrature method).
 velocityField (underworld.mesh.MeshVariable) – Solution field for darcy flow velocity. Optional.
 swarmVarVelocity (undeworld.swarm.SwarmVariable) – If a swarm variable is provided, the velocity calculated on the swarm will be stored. This is the most representative velocity data object, as the velocity calculation occurs on the swarm, away from mesh nodes. Optional.
Notes
Constructor must be called collectively by all processes.

fn_bodyforce
¶ The heating function. You may change this function directly via this property.

fn_diffusivity
¶ The diffusivity function. You may change this function directly via this property.

class
underworld.systems.
SteadyStateHeat
(temperatureField, fn_diffusivity, fn_heating=0.0, voronoi_swarm=None, conditions=[], gauss_swarm=None, _removeBCs=True, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
This class provides functionality for a discrete representation of the steady state heat equation.
The class uses a standard Galerkin finite element method to construct a system of linear equations which may then be solved using an object of the underworld.system.Solver class.
The underlying element types are determined by the supporting mesh used for the ‘temperatureField’.
The strong form of the given boundary value problem, for \(f\), \(h\) and \(h\) given, is
\[\begin{split}\begin{align} q_i =&  \alpha \, u_{,i} & \\ q_{i,i} =& \: f & \text{ in } \Omega \\ u =& \: g & \text{ on } \Gamma_g \\ q_i n_i =& \: h & \text{ on } \Gamma_h \\ \end{align}\end{split}\]where, \(\alpha\) is the diffusivity, \(u\) is the temperature, \(f\) is a source term, \(g\) is the Dirichlet condition, and \(h\) is a Neumann condition. The problem boundary, \(\Gamma\), admits the decomposition \(\Gamma=\Gamma_g\cup\Gamma_h\) where \(\emptyset=\Gamma_g\cap\Gamma_h\). The equivalent weak form is:
\[\int_{\Omega} w_{,i} \, q_i \, d \Omega = \int_{\Omega} w \, f \, d\Omega + \int_{\Gamma_h} w \, h \, d \Gamma\]where we must find \(u\) which satisfies the above for all \(w\) in some variational space.
Parameters:  temperatureField (underworld.mesh.MeshVariable) – The solution field for temperature.
 fn_diffusivity (underworld.function.Function) – The function that defines the diffusivity across the domain.
 fn_heating (underworld.function.Function) – A function that defines the heating across the domain. Optional.
 voronoi_swarm (underworld.swarm.Swarm) – If a voronoi_swarm is provided, voronoi type numerical integration is utilised. The provided swarm is used as the basis for the voronoi integration. If no voronoi_swarm is provided, Gauss integration is used.
 gauss_swarm (underworld.swarm.GaussIntegrationSwarm) – If provided this gauss_swarm will be used for (gaussian) numerical integration rather than a default gauss integration swarm that is automatically generated and dependent on the element order of the mesh. NB: if a voronoi_swarm is defined it OVERRIDES this gauss_swarm as the preferred integration swarm (quadrature method).
 conditions (underworld.conditions.SystemCondition) – Numerical conditions to impose on the system. This should be supplied as the condition itself, or a list object containing the conditions.
Notes
Constructor must be called collectively by all processes.
Example
Setup a basic thermal system:
>>> linearMesh = uw.mesh.FeMesh_Cartesian( elementType='Q1/dQ0', elementRes=(4,4), minCoord=(0.,0.), maxCoord=(1.,1.) ) >>> tField = uw.mesh.MeshVariable( linearMesh, 1 ) >>> topNodes = linearMesh.specialSets["MaxJ_VertexSet"] >>> bottomNodes = linearMesh.specialSets["MinJ_VertexSet"] >>> tbcs = uw.conditions.DirichletCondition(tField, topNodes + bottomNodes) >>> tField.data[topNodes.data] = 0.0 >>> tField.data[bottomNodes.data] = 1.0 >>> tSystem = uw.systems.SteadyStateHeat(temperatureField=tField, fn_diffusivity=1.0, conditions=[tbcs])
Example with non diffusivity:
>>> k = tField + 1.0 >>> tSystem = uw.systems.SteadyStateHeat(temperatureField=tField, fn_diffusivity=k, conditions=[tbcs]) >>> solver = uw.systems.Solver(tSystem) >>> solver.solve() Traceback (most recent call last): ... RuntimeError: Nonlinearity detected. Diffusivity function depends on the temperature field provided to the system. Please set the 'nonLinearIterate' solve parameter to 'True' or 'False' to continue. >>> solver.solve(nonLinearIterate=True)

fn_diffusivity
¶ The diffusivity function. You may change this function directly via this property.

fn_heating
¶ The heating function. You may change this function directly via this property.

class
underworld.systems.
Stokes
(velocityField, pressureField, fn_viscosity, fn_bodyforce=None, fn_one_on_lambda=None, fn_source=None, voronoi_swarm=None, conditions=[], gauss_swarm=None, _removeBCs=True, _fn_viscosity2=None, _fn_director=None, fn_stresshistory=None, _fn_stresshistory=None, _fn_v0=None, _fn_p0=None, _fn_fssa=None, _callback_post_solve=None, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
This class provides functionality for a discrete representation of the Stokes flow equations.
Specifically, the class uses a mixed finite element method to construct a system of linear equations which may then be solved using an object of the underworld.system.Solver class.
The underlying element types are determined by the supporting mesh used for the ‘velocityField’ and ‘pressureField’ parameters.
The strong form of the given boundary value problem, for \(f\), \(g\) and \(h\) given, is
\[\begin{split}\begin{align} \sigma_{ij,j} + f_i =& \: 0 & \text{ in } \Omega \\ u_{k,k} + \frac{p}{\lambda} =& \: H & \text{ in } \Omega \\ u_i =& \: g_i & \text{ on } \Gamma_{g_i} \\ \sigma_{ij}n_j =& \: h_i & \text{ on } \Gamma_{h_i} \\ \end{align}\end{split}\]where,
 \(\sigma_{i,j}\) is the stress tensor
 \(u_i\) is the velocity,
 \(p\) is the pressure,
 \(f_i\) is a body force,
 \(\lambda\) is pseudo compressibility factor,
 \(H\) is the compressible equation source term,
 \(g_i\) are the velocity boundary conditions (DirichletCondition)
 \(h_i\) are the traction boundary conditions (NeumannCondition).
The problem boundary, \(\Gamma\), admits the decompositions \(\Gamma=\Gamma_{g_i}\cup\Gamma_{h_i}\) where \(\emptyset=\Gamma_{g_i}\cap\Gamma_{h_i}\). The equivalent weak form is:
\[\int_{\Omega} w_{(i,j)} \sigma_{ij} \, d \Omega = \int_{\Omega} w_i \, f_i \, d\Omega + \sum_{j=1}^{n_{sd}} \int_{\Gamma_{h_j}} w_i \, h_i \, d \Gamma\]where we must find \(u\) which satisfies the above for all \(w\) in some variational space.
Parameters:  velocityField (underworld.mesh.MeshVariable) – Variable used to record system velocity.
 pressureField (underworld.mesh.MeshVariable) – Variable used to record system pressure.
 fn_viscosity (underworld.function.Function) – Function which reports a viscosity value. Function must return scalar float values.
 fn_bodyforce (underworld.function.Function, Default = None) – Function which reports a body force for the system. Function must return float values of identical dimensionality to the provided velocity variable.
 fn_one_on_lambda (underworld.function.Function, Default = None) – Pseudocompressibility factor. Note that nonzero values are incompatible with the ‘penalty’ stokes solver. Ensure a ‘penalty’ equal to 0 is used if this function is nonzero. By default this is the case.
 fn_source (underworld.function.Function, Default = None) – Mass source term. Check fn_one_on_lambda for usage caveats.
 fn_stresshistory (underworld.function.Function, Default = None) – Function which defines the stress history term used for viscoelasticity. Function is a vector of size 3 (dim=2) or 6 (dim=3) representing a symetric tensor.
 voronoi_swarm (underworld.swarm.Swarm) – If a voronoi_swarm is provided, voronoi type numerical integration is utilised. The provided swarm is used as the basis for the voronoi integration. If no voronoi_swarm is provided, Gauss integration is used.
 gauss_swarm (underworld.swarm.GaussIntegrationSwarm) – If provided this gauss_swarm will be used for (gaussian) numerical integration rather than a default gauss integration swarm that is automatically generated and dependent on the element order of the mesh. NB: if a voronoi_swarm is defined it OVERRIDES this gauss_swarm as the preferred integration swarm (quadrature method).
 conditions (underworld.conditions.SystemCondition) – Numerical conditions to impose on the system. This should be supplied as the condition itself, or a list object containing the conditions.
Notes
Constructor must be called by collectively all processes.

eqResiduals
¶ Returns the stokes flow equations’ residuals from the latest solve. Residual calculations use the matrices and vectors of the discretised problem. The residuals correspond to the momentum equation and the continuity equation.
Returns: r1 is the momentum equation residual r2 is the continuity equation residual Return type: (r1, r2)  2 tuple of doubles Notes
This method must be called collectively by all processes.

fn_bodyforce
¶ The body force function. You may change this function directly via this property.

fn_one_on_lambda
¶ A bulk viscosity parameter

fn_source
¶ The volumetric source term function. You may change this function directly via this property.

fn_viscosity
¶ The viscosity function. You may change this function directly via this property.

stokes_callback
¶ Return the callback function used by this system

class
underworld.systems.
StokesSolver
(stokesSLE, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
The Block Stokes Schur Complement Solver: This solves the saddlepoint system
\[\begin{split}\begin{bmatrix} K & G \\ G^T & C \end{bmatrix} \begin{bmatrix} u \\ p \end{bmatrix} = \begin{bmatrix}f \\ h \end{bmatrix}\end{split}\]via a Schur complement method.
We first solve:
(1)¶\[S p= G^T K^{1} f  h,\]where \(S = G^T K^{1} GC\)
Then we backsolve for the velocity:
(2)¶\[K u = f  G p.\]The effect of \(K^{1}\) in (1) is obtained via a KSPSolve in PETSc. This has the prefix ‘A11’ (often called the ‘inner’ solve)
The solve in (1) for the pressure has prefix ‘scr’.
Assuming the returned solver is called ‘solver’, it is possible to configure these solves individually via the solver.options.A11 and solver.options.scr dictionaries.
Try help(solver.options.A11) for some details.
Common configurations are provided via the set_inner_method() method.
help(solver.set_inner_method) for more.
For more advanced configurations use the solver.options.A11/scr dictionaries directly.
help(solver.options) to see more.

set_inner_method
(solve_type='mg')[source]¶ Configure velocity/inner solver (A11 PETSc prefix).
Available options below. Note that associated solver software (for example mumps) must be installed.
 mg : Geometric multigrid (default).
 nomg : Disables multigrid.
 lu : LU direct solver (serial only).
 mumps : MUMPS parallel direct solver.
 superludist : SuperLU parallel direct solver.
 superlu : SuperLU direct solver (serial only).

set_mg_levels
(levels)[source]¶ Set the number of multigrid levels manually. It is set automatically by default.

set_penalty
(penalty)[source]¶ By setting the penalty, the Augmented Lagrangian Method is used as the solve. This method is not recommended for normal use as there is additional memory and cpu overhead. This method can often help improve convergence issues for problems with large viscosity contrasts that are having trouble converging.
A penalty of roughly 0.1 of the maximum viscosity contrast is not a bad place to start as a rule of thumb. (check notes/paper)

solve
(nonLinearIterate=None, nonLinearTolerance=0.01, nonLinearKillNonConvergent=False, nonLinearMinIterations=1, nonLinearMaxIterations=500, callback_post_solve=None, print_stats=False, reinitialise=True, **kwargs)[source]¶ Solve the stokes system
Parameters:  nonLinearIterate (bool) – True will perform non linear iterations iterations, False (or 0) will not
 nonLinearTolerance (float, Default=1.0e2) – Relative tolerance criterion for the change in the velocity field
 nonLinearMaxIterations (int, Default=500) – Maximum number of non linear iteration to perform
 callback_post_sovle (func, Default=None) – Optional callback function to be performed at the end of a linear solve iteration. Commonly this will be used to perform operations between non linear iterations, for example, calibrating the pressure solution or removing the system null space.
 print_stats (bool, Default=False) – Print out solver iteration and timing counts per solver
 reinitialise (bool, Default=True,) – Rebuild the system discretisation storage (location matrix/petsc mats & vecs) and repopulate, if available, the stokes voronio swarm before the system is solved.


class
underworld.systems.
SwarmAdvector
(velocityField, swarm, order=2, **kwargs)[source]¶ Bases:
underworld.systems._timeintegration.TimeIntegration
Objects of this class advect a swarm through time using the provided velocity field.
Parameters:  velocityField (underworld.mesh.MeshVariable) – The MeshVariable field used for evaluating the velocity field that advects the swarm particles
 swarm (underworld.swarm.Swarm) – Particle swarm that will be advected by the given velocity field

integrate
(dt, update_owners=True)[source]¶ Integrate the associated swarm in time, by dt, using the velocityfield that is associated with this class
Parameters:  dt (double) – The timestep to use in the intergration
 update_owners (bool) – If this is set to False, particle ownership (which element owns a particular particle) is not updated after advection. This is often necessary when both the mesh and particles are advecting simutaneously.
Example
>>> import underworld as uw >>> import numpy as np >>> from underworld import function as fn >>> dim=2; >>> elementMesh = uw.mesh.FeMesh_Cartesian(elementType="Q1/dQ0", elementRes=(9,9), minCoord=(1.,1.), maxCoord=(1.,1.)) >>> velocityField = uw.mesh.MeshVariable( mesh=elementMesh, nodeDofCount=dim ) >>> swarm = uw.swarm.Swarm(mesh=elementMesh) >>> particle = np.zeros((1,2)) >>> particle[0] = [0.2,0.2] >>> swarm.add_particles_with_coordinates(particle) array([0], dtype=int32) >>> velocityField.data[:]=[1.0,1.0] >>> swarmAdvector = uw.systems.SwarmAdvector(velocityField=velocityField, swarm=swarm, order=2) >>> dt=swarmAdvector.get_max_dt() >>> swarmAdvector.integrate(dt) >>> np.allclose(swarm.particleCoordinates.data[0], [ 0.27856742, 0.12143258], rtol=1e4) True

class
underworld.systems.
TimeIntegration
(order, **kwargs)[source]¶ Bases:
underworld._stgermain.StgCompoundComponent
Abstract class for integrating numerical objects (fields, swarms, etc.) in time.
The integration algorithm is a modified Runge Kutta method that only evaluates midpoint information varying in space  using only the present timestep solution. The order of the integration used can be 1,2,4
Parameters: order (int {1,2,4}) – Defines the numerical order ‘in space’ of the Runge Kutta like integration scheme. 
dt
¶ Time integrator timestep size.

time
¶ Time integrator time value.
