UWGeodynamics User Guide¶

The Jupyter notebook¶

The Jupyter notebook provides a powerful environment for the development and analysis of Underworld models. Underworld and UWGeodynamics recommend using Jupyter notebooks for the development of geodynamic models.

If you are not familiar with Jupyter notebooks, we suggest you follow a quick introduction here.

Where to find documentation?¶

Additional documentation and function specific documentation can be find in the python doctrings. You can acces them in the Jupyter notebook by prepending or appending the method, variable or function with ?.

Design principles¶

import UWGeodynamics¶

UWGeodynamics can be imported as follow:

>>> from underworld import UWGeodynamics as GEO
...

Visualization¶

We provide access to a wrapper around Lavavu , readily available from Underworld.

The visualisation module can be imported as follow

>>> from underworld import visualisation as vis

Warning

Although many plotting modules are available, we strongly encourage people to use the visualisation module. It integrates very well inside the Jupyter notebook, is parallel safe, and can take Underworld function as arguments.

Warning

We provide some basic examples. Look at the Lavavu documentation for more details.

Simple examples:¶

Plot Material Field or any field / variable defined on the swarm (e.g. plasticstrain, viscosityField, densityField etc.):

>>> from underworld import visualisation as vis
>>> Model = GEO.Model()
>>> Fig = vis.Figure(figsize=(1200,400), title="Material Field")
>>> Fig.Points(Model.swarm, Model.materialField, fn_size=3.0)
...
>>> Fig.show()
>>> Fig.save("MaterialField.png")
'MaterialField.png'

Plot Temperature Field or any field / variable defined on the mesh (e.g. temperature, pressureField, velocityField, strainRateField) as well as projected swarm field / variables (e.g. projMaterialField, projViscosityField etc.)

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> Fig = vis.Figure(figsize=(1200,400), title="Temperature")
>>> Fig.Surface(Model.mesh, GEO.dim(Model.temperature, u.degK))
...
>>> Fig.show()
>>> Fig.save("Temperature.png")
'Temperature.png'

Note

Fields can be dimensionalized using the GEO.dimensionalise function (see below)

Plot Velocity Fields or any vector field.

The example below plots a temperature field with the velocity vectors on top:

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> Fig = vis.Figure(figsize=(1200,400), title="Velocity")
>>> Fig.Surface(Model.mesh, GEO.dim(Model.temperature, u.degK))
...
>>> Fig.VectorArrows(Model.mesh, Model.velocityField)
...
>>> Fig.show()
>>> Fig.save("VelocityField.png")
'VelocityField.png'

Working with units¶

UWGeodynamics uses Pint, a Python package to define, operate and manipulate physical quantities (A numerical value with unit of measurement). Pint is a very powerful package that handles conversion and operation between units.

We recommend using SI units but other systems are also available.

Pint Unit Registry can be used as follow:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.UnitRegistry

or simply

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u

You can have a quick overview of all the units available by hitting tab after the . of the u object.

Quantities can then be defined as follow:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> length = 100. * u.kilometre
>>> width = 50. * u.kilometre
>>> gravity = 9.81 * u.metre / u.second**2

Pint offers the possibility to append a prefix to the units. 1 million years can thus be defined as follow:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> length = 1.0 * u.megayear

Note

Unit abbreviation is also possible u.km is equivalent to u.kilometer. You can refer to the Pint documentation for all abbreviations available.

Model Scaling¶

Model can be scaled using a series of scaling coefficients

>>> from underworld import UWGeodynamics as GEO
>>> GEO.scaling_coefficients
...

The default scaling coefficients are defined as follow:

Dimension	value
[mass]	1.0 kilogram
[length]	1.0 metre
[temperature]	1.0 kelvin
[time]	1.0 second
[substance]	1.0 mole

The scaling value can be changed by accessing each scaling coefficient as follow

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u

>>> GEO.scaling_coefficients["[length]"] = 3. * u.kilometre
>>> GEO.scaling_coefficients["[mass]"] = 4. * u.kilogram
>>> GEO.scaling_coefficients["[temperature]"] = 273.15 * u.degK
>>> GEO.scaling_coefficients["[time]"] = 300. * u.years

The unit entered are checked internally and an error is raised if the units are incompatible. The value is automatically converted to the base units (metre, second, degree, etc).

To scale a model, the user must define a series of characteristic physical values and assign them to the scaling object.

Arguments with units will be scaled by the UWGeodynamics functions.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u

>>> KL = 100 * u.kilometre
>>> Kt = 1. * u.year
>>> KM = 3000. * u.kilogram
>>> KT = 1200. * u.degK

>>> GEO.scaling_coefficients["[length]"] = KL
>>> GEO.scaling_coefficients["[time]"] = Kt
>>> GEO.scaling_coefficients["[mass]"]= KM
>>> GEO.scaling_coefficients["[temperature]"] = KT

dimensionalise / non-dimensionalise¶

We provide 2 functions GEO.non_dimensionalise and GEO.dimensionalise to convert between non-dimensional and dimensional values. The function are also available respectively as GEO.nd and GEO.dim.

Example:

define a length of 300 kilometres.
use the GEO.nd function to scale it.
convert the value back to SI units.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u

>>> GEO.scaling_coefficients["[length]"] = 300. * u.kilometre

>>> length = 300. * u.kilometre
>>> scaled_length = GEO.nd(length)
>>> print(scaled_length)
1.0
>>> length_metres = GEO.dimensionalise(scaled_length, u.metre)
>>> print(length_metres)
300000.0 meter

The Model object¶

The central element or “object” of the UWGeodynamics module is the Model object.

It has several uses:

It defines the extent and the outside geometry of your problem.
It works as a container for the field variables.

It basically defines the universe on which you are going to apply physical rules (Gravity field, boundary condition, composition, temperature etc.) It is the equivalent of the box in which you would put the sand and silicon if you were to build an analog experiment in a lab. One important difference is that the “box” his not empty, it is populated with particles that have already some properties. The properties are changed by defining new materials.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model(elementRes=(64, 64),
...                   minCoord=(0. * u.kilometre, 0. * u.kilometre),
...                   maxCoord=(64. * u.kilometre, 64. * u.kilometre))

The Material object¶

The UWGeodynamics module is designed around the idea of materials, which are essentially a way to define physical properties across the Model domain.

Predefined Material objects¶

A library of predefined material is available through the MaterialRegistry object:

from underworld import UWGeodynamics as GEO

materials_database = GEO.MaterialRegistry()

Note

The MaterialRegistry object can import a database of materials from a json file by passing its path as argument. The `default json`__ file can be used as an example.

User defined¶

Materials are defined using the Material object as follow:

>>> from underworld import UWGeodynamics as GEO

>>> crust = GEO.Material(name="Crust")

Typing the name of the material in an empty cell will return a table which summarizes the property of the material:

As you can see, most of the property are undefined.

They are several ways to define the physical parametres of our Material.

The first one is to add them directly when creating the object itself:

>>> from underworld import UWGeodynamics as GEO

>>> u = GEO.u
>>> crust = GEO.Material(name="Crust", density=3000*u.kilogram/u.metre**3)

The second option is to change the property after creating the Material:

>>> from underworld import UWGeodynamics as GEO

>>> u = GEO.u
>>> crust = GEO.Material(name="Crust")
>>> crust.density = 3000. * u.kilogram / u.metre **3

The second option is often easier to read.

Warning

UWGeodynamics contains some basic dimensionality checks. Entering wrong units will raise an error

Material can be added to a model as follow:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()
>>> crust = Model.add_material(name="Crust")

Although optional, it is a good idea to give a name to the material. The Model.add_material method will return a Material object. That object is a python object that will then be used to define the property of the material.

Material Attributes¶

The Material object comes with a series of attribute that can be used to define its physical behavior.

Materials attributes¶
Name	Description
shape	Initial Geometrical Representation
density	Density
diffusivity	Thermal Diffusivity
capacity	Thermal Capacity
radiogenicHeatProd	Radiogenic Heat Production
viscosity	Viscous behavior
plasticity	Plastic behavior
elasticity	Elastic behavior
minViscosity	Minimum Viscosity allowed
maxViscosity	Maximum Viscosity allowed
stressLimiter	Maximum sustainable stress
healingRate	Plastic Strain Healing Rate
solidus	Solidus
liquidus	Liquidus
latentHeatFusion	Latent Heat Fusion (Enthalpy of Fusion)
meltExpansion	Melt Expansion
meltFraction	Initial Melt Fraction
meltFractionLimit	Maximum Fraction of Melt
viscosityChange	Change in Viscosity over Melt Fraction range
viscosityChangeX1	Melt Fraction Range begin
viscosityChangeX2	Melt Fraction Range end

Examples

>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.density = 200. * u.kg / u.m**3
>>> myMaterial = GEO.Material(name="My Material")
>>> myMaterial.density = 3000 * u.kilogram / u.metre**3
>>> myMaterial.viscosity = 1e19 * u.pascal * u.second
>>> myMaterial.radiogenicHeatProd = 0.7 * u.microwatt / u.metre**3
>>> myMaterial.diffusivity = 1.0e-6 * u.metre**2 / u.second

Global properties¶

The user can define attributes on the Model itself. The values will be used as global values for materials with undefined attributes

Example

>>> u = GEO.u
>>> Model = GEO.Model()
>>> Model.density = 200. * u.kg / u.m**3
>>> myMaterial = GEO.Material(name="My Material")

The density of myMaterial will default to 200. kilogram / cubic metre unless its density attribute is explicitly specified.

Material shape¶

The shape attribute essentially describes the initial location of a material. It is used to build the initial geometry of the model.

There are a range of available/pre-defined shapes

Layer (2D/3D)
Polygon (2D)
Box (2D)
Disk (2D)
Spheres (3D)
Annulus (2D)
CombinedShape (Combination of any of the above) (2D)
HalfSpace (3D)

Layer

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> shape = GEO.shapes.Layer(top=30.*u.kilometre, bottom=0.*u.kilometre)
>>> material = Model.add_material(name="Material", shape=shape)

>>> Fig = vis.Figure(figsize=(1200,400))
>>> Fig.Points(Model.swarm, Model.materialField)
...
>>> Fig.show()

Polygon

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> polygon = GEO.shapes.Polygon(vertices=[(10.* u.kilometre, 10.*u.kilometre),
...                                        (20.* u.kilometre, 35.*u.kilometre),
...                                        (35.* u.kilometre, 5.*u.kilometre)])
>>> material = Model.add_material(name="Material", shape=polygon)

>>> Fig = vis.Figure(figsize=(1200,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

Box

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> box = GEO.shapes.Box(top=10.* u.kilometre, bottom=5*u.kilometre,
...                      minX=10.*u.kilometre, maxX=15*u.kilometre)
>>> material = Model.add_material(name="Material", shape=box)

>>> Fig = vis.Figure(figsize=(1200,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

Disk

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> disk = GEO.shapes.Disk(center=(32. * u.kilometre, 32. * u.kilometre),
...                        radius=10.*u.kilometre)
>>> material = Model.add_material(name="Material", shape=disk)

>>> Fig = vis.Figure(figsize=(1200,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

Sphere (3D)

>>> from underworld import UWGeodynamics as GEO

>>> u = GEO.u
>>> Model = GEO.Model(elementRes=(16, 16, 16),
...                   minCoord=(-1. * u.m, -1. * u.m, -50. * u.cm),
...                   maxCoord=(1. * u.m, 1. * u.m, 50. * u.cm))

>>> sphereShape = GEO.shapes.Sphere(center=(0., 0., 20.*u.centimetre),
                                    radius=20. * u.centimetre))

Annulus

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> annulus = GEO.shapes.Annulus(center=(35.*u.kilometre, 50.*u.kilometre),
...                              r1=5.*u.kilometre,
...                              r2=10.*u.kilometre)
>>> material = Model.add_material(name="Material", shape=annulus)

>>> Fig = vis.Figure(figsize=(400,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

CombinedShape

Several shapes can be combined to form a material shape:

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> disk1 = GEO.shapes.Disk(center=(10. * u.kilometre, 10. * u.kilometre),
...                         radius=10.*u.kilometre)
>>> disk2 = GEO.shapes.Disk(center=(20. * u.kilometre, 20. * u.kilometre),
...                         radius=5.*u.kilometre)

>>> shape = disk1 | disk2
>>> material = Model.add_material(name="Material", shape=shape)

>>> Fig = vis.Figure(figsize=(400,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

You can also take the intersection of some shapes:

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()
>>> disk1 = GEO.shapes.Disk(center=(32. * u.kilometre, 32. * u.kilometre),
...                         radius=10.*u.kilometre)
>>> disk2 = GEO.shapes.Disk(center=(32. * u.kilometre, 22. * u.kilometre),
...                         radius=10.*u.kilometre)

>>> shape = disk1 & disk2
>>> material = Model.add_material(name="Material", shape=shape)

>>> Fig = vis.Figure(figsize=(400,400))
>>> Fig.Points(Model.swarm, Model.materialField)
>>> Fig.show()

HalfSpace

HalfSpaces can be used to divide space in 2 domains. The divide is a plan defined by its normal vector. The convention is to keep the domain opposite to direction defined by the normal vector.

Note

HalfSpaces can be combined to define 3D shapes / volumes.

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.UnitRegistry

>>> Model = GEO.Model(elementRes=(34, 34, 12),
...                   gravity=(0., 0., -9.81 * u.m / u.s**2),
...                   minCoord=(0. * u.km, 0. * u.km, -2880. * u.km),
...                   maxCoord=(9000. * u.km, 2000. * u.km, 20. * u.km))

>>> halfspace1 = GEO.shapes.HalfSpace(normal=(-1.,0.,1.), origin=(4000. * u.km, 0. * u.km, -1000. * u.km))
>>> halfspace2 = GEO.shapes.HalfSpace(normal=(0.,0.,1.), origin=(7000. * u.km, 1000. * u.km, 0. * u.km))
>>> halfspace3 = GEO.shapes.HalfSpace(normal=(1.,0.,0.), origin=(9000. * u.km, 1000. * u.km, -500. * u.km))
>>> halfspace4 = GEO.shapes.HalfSpace(normal=(0.,0.,-1.), origin=(6500. * u.km, 1000. * u.km, -1000. * u.km))

>>> compositeShape = halfspace1 & halfspace2 & halfspace3 & halfspace4
>>> polygon= Model.add_material(name="polygon", shape=compositeShape)

>>> Fig = vis.Figure()
>>> Fig.Points(Model.swarm, Model.materialField, cullface=False, opacity=1.)
>>> Fig.Mesh(Model.mesh)
>>> viewer = Fig.viewer(resolution=(1200,600))
>>> viewer = Fig.viewer(axis=True)
>>> viewer.rotatex(-70)
>>> viewer.rotatey(-10)
>>> viewer.window()

Multiple materials

You can add as many materials as needed:

>>> from underworld import UWGeodynamics as GEO
>>> from underworld import visualisation as vis

>>> u = GEO.u
>>> Model = GEO.Model()

>>> shape = GEO.shapes.Layer(top=30.*u.kilometre, bottom=0.*u.kilometre)
>>> material1 = Model.add_material(name="Material 1", shape=shape)

>>> polygon = GEO.shapes.Polygon(vertices=[(10.* u.kilometre, 10.*u.kilometre),
...                                        (20.* u.kilometre, 35.*u.kilometre),
...                                        (35.* u.kilometre, 5.*u.kilometre)])

>>> material2 = Model.add_material(name="Material 2", shape=polygon)

>>> Fig = vis.Figure(figsize=(400,400))
>>> Fig.Points(Model.swarm, Model.materialField, fn_size=3.)
...
>>> Fig.show()
>>> Fig.save("multiple_materials.png")
'multiple_materials.png'

Temperature and Pressure dependent densities¶

Temperature and Pressure dependent densities can be assigned to a Material using the GEO.LinearDensity function which calculates:

(1)¶\[`\rho = \rho_0 (1 + (\beta * \delta P) - (\alpha * \delta T))`\]

where $\rho$ is the reference density, $\beta$ a factor, $\delta P$ the difference between the pressure and the reference pressure, $\alpha$ is the thermal expansivity and $\delta T$ is the difference between the temperature and the reference temperature.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> material1 = Model.add_material(name="Material 1", shape=shape)
>>> material1.density = GEO.LinearDensity(reference_density=3370. * u.kilogram / u.metre**3,
...                                       thermalExpansivity= 2.8e-5 * u.kelvin**-1,
...                                       beta=1.0)

Rheologies¶

Newtonian Rheology¶

A newtonian rheology can be applied by assigning a viscosity to a already defined material

>>> from underworld import UWGeodynamics as GEO

>>> myMaterial = GEO.Material(name="Newtonian Material")
>>> myMaterial.viscosity = 1e19 * u.pascal * u.second

Non-Newtonian Rheology¶

Deformation of materials on long timescale is predominantly achieved through viscous diffusion and dislocation creep. Those processes can be expressed using the following equation:

\[\eta_{\text{eff}}^{vcreep} = \frac{1}{2}A^{\frac{-1}{n}} \dot{\epsilon}^{\frac{(1-n)}{n}}d^{\frac{m}{n}}\exp{\left(\frac{E + PV}{nRT}\right)}\]

with A the prefactor, $\dot{\epsilon}$ the square root of the second invariant of the deviatoric strain rate tensor, d the grain size, p the grain size exponent, E the activation energy, P the pressure, V the activation volume, n the stress exponent, R the Gas Constant and T the temperature.

UWGeodynamics provides a library of commonly used Viscous Creep Flow Laws. These can be accessed using the GEO.ViscousCreepRegistry:

Note

The ViscousCreepRegistry object can import a database of rheologies from a json file by passing its path as argument. The `default json`__ file can be find here and can be used as an example.

Example:

>>> from underworld import UWGeodynamics as GEO
>>> material = GEO.Material(name="Material")

>>> rh = GEO.ViscousCreepRegistry()
>>> material.viscosity = rh.Wet_Quartz_Dislocation_Gleason_and_Tullis_1995

You can scale viscosity by using a multiplier. For example to make the Gleason and Tullis, 1995 rheology 30 times stronger you can do:

>>> from underworld import UWGeodynamics as GEO
>>> material = GEO.Material(name="Material")

>>> rh = GEO.ViscousCreepRegistry()
>>> material.viscosity = 30 * rh.Wet_Quartz_Dislocation_Gleason_and_Tullis_1995

The user can of course define their own ViscousCreep rheology.

>>> from underworld import UWGeodynamics as GEO
>>> viscosity = GEO.ViscousCreep(preExponentialFactor=1.0,
...                              stressExponent=1.0,
...                              activationVolume=0.,
...                              activationEnergy=200 * u.kilojoules,
...                              waterFugacity=0.0,
...                              grainSize=0.0,
...                              meltFraction=0.,
...                              grainSizeExponent=0.,
...                              waterFugacityExponent=0.,
...                              meltFractionFactor=0.0,
...                              f=1.0)

Single parametres can then be modified.

viscosity.activationEnergy = 300. * u.kilojoule

Composite Viscosity¶

Material viscosity can be assigned a combination of viscosities. The effective viscosity is calculated as the harmonic mean of all viscosities.

This is useful to combine diffusion and dislocation creep:

\[\eta_{\text{eff}}^{vcreep} = \left(\frac{1}{\eta_{\text{eff}}^{\text{diff}}} + \frac{1}{\eta_{\text{eff}}^{\text{disl}}}\right)\]

>>> from underworld import UWGeodynamics as GEO

>>> viscosity1 = GEO.ViscousCreep(preExponentialFactor=1.0,
...                              stressExponent=1.0,
...                              activationVolume=0.,
...                              activationEnergy=200 * u.kilojoules,
...                              waterFugacity=0.0,
...                              grainSize=0.0,
...                              meltFraction=0.,
...                              grainSizeExponent=0.,
...                              waterFugacityExponent=0.,
...                              meltFractionFactor=0.0,
...                              f=1.0)


>>> viscosity2 = GEO.ViscousCreep(preExponentialFactor=1.0,
...                              stressExponent=1.0,
...                              activationVolume=0.,
...                              activationEnergy=200 * u.kilojoules,
...                              waterFugacity=0.0,
...                              grainSize=0.0,
...                              meltFraction=0.,
...                              grainSizeExponent=0.,
...                              waterFugacityExponent=0.,
...                              meltFractionFactor=0.0,
...                              f=1.0)
>>> combined_viscosity = GEO.CompositeViscosity([viscosity1, viscosity2])

Plastic Behavior (Yield)¶

Plastic yielding can be added and will result in rescaling the effective viscosity for a stress limited to the yield stress of the material.

The effective plastic viscosity is given by:

\[\eta_{\text{eff}} = \frac{\sigma_y}{2\dot{\epsilon}}\]

Where $\dot{\epsilon}$ is the second invariant of the strain rate tensor defined as $\dot{\epsilon}=\sqrt{\frac{1}{2}\dot{\epsilon}_{ij}\dot{\epsilon}_{ij}}$ The yield value $\sigma_y$ is defined using a Drucker-Prager yield-criterion:

(2)¶\[ \begin{align}\begin{aligned}\sigma_y = C \cos\phi + \sin\phi P \quad \text{(2D)}\\\sigma_y = \frac{6C\cos\phi}{\sqrt3(3-\sin\phi)} + \frac{6\sin\phi P}{\sqrt3(3-\sin\phi)} \quad \text{(3D)}\end{aligned}\end{align} \]

Setting the friction angle $\phi=0$ gives the von Mises Criterion. In 2D, equation (2) corresponds to the Mohr-Coulomb criterion, while in 3D it circumscribes the Mohr-Coulomb yield surface.

Linear cohesion and friction weakening can be added by defining their initial and final values over a range of accumulated plastic strain.

As with Viscous Creep, we also provide a registry of commmonly used plastic behaviors. They can be accessed using the GEO.PlasticityRegistry registry.

Note

The PlasticityRegistry object can import a database of plasticity from a json file by passing its path as argument. The `default json`__ file can be find here and can be used as an example.

Users can define their own parametres:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()
>>> material = Model.add_material()

>>> material.plasticity = GEO.DruckerPrager(
...     cohesion=10. * u.megapascal,
...     cohesionAfterSoftening=10. * u.megapascal,
...     frictionCoefficient = 0.3,
...     frictionAfterSoftening = 0.2,
...     epsilon1=0.5,
...     epsilon2=1.5)

Viscous Creep and Plastic yielding are combined by assuming that they act in parallel as independent processes. The effective viscoplastic viscosity is calculated as:

\[\eta_{\text{eff}^{\text{vp}}} = \min{(\eta_{\text{eff}}^{\text{vcreep}}, \eta_{\text{eff}}^{\text{pl}})}\]

Elasticity¶

Elastic behavior can be added to a material:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()
>>> material = Model.add_material()

>>> material.elasticity(shear_modulus=10e9 * u.pascal,
                        observation_time=10000 * u.year)

Simple phase change¶

One can change the property of one material to another depending on some time, tepmerature, pressure etc. criteria. This is not a phase change sensu-stricto but this allows for simple change such as transition from mantle to oceanic-crust behavior or simply air to water…

Warning

Phase changes can only occur between predefined material. If you plan to add a material during the Model run, you will have to define it beforehand.

In the following example we change air into water when the air particles move below the 0. level.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> air = Model.add_material(name="air")
>>> water = Model.add_material(name="water")
>>> air.phase_changes = GEO.PhaseChange((Model.y < 0.), water.index)

The above example essentially fills the basins with water. For such a specific purpose you can use the WaterFill class.

>>> air = Model.add_material(name="air")
>>> water = Model.add_material(name="water")
>>> air.phase_changes = GEO.WaterFill(sealevel=0., result=water)

This is easier to read but equivalent.

Melt¶

Materials can be assigned a Solidus and a Liquidus defined as polynomial function of temperature. This allows to calculate the fraction of melt present in the material.

Warning

There is no seggregation of the melt from its source.

A registry of Solidii and Liquidii are available:

>>> from underworld import UWGeodynamics as GEO

>>> solidii = GEO.SolidusRegistry()
>>> crust_solidus = solidii.Crustal_Solidus

>>> liquidii = GEO.LiquidusRegistry()
>>> crust_liquidus = liquidii.Crustal_Liquidus

The percentage of melt results in a linear decrease of the viscosity of a factor viscosityChange over the viscosityChangeX1 - viscosityChangeX2 melt fraction interval.

The latent heat of fusion is embedded in the energy equation and affects the temperature field of the Model.

The meltExpansion factor affects the density of the materials and equation (1) becomes:

\[ \begin{align}\begin{aligned}:label: linearDensityMelt\\`\rho = \rho_0 (1 + (\beta * \delta P) - (\alpha * \delta T) - \gamma F) `\end{aligned}\end{align} \]

with gamma the factor of melt expansion and F the fraction of melt.

The following example prescribes a decrease in viscosity of 3 order of magnitude over a range of 0.15 to 0.30 fraction of melt. The fraction of the melt is limited to 0.3

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> continentalcrust = Model.add_material(name="Continental Crust")
>>> continentalcrust.radiogenicHeatProd = 7.67e-7 * u.watt / u.meter**3
>>> continentalcrust.density  = 2720. * u.kilogram / u.metre**3

>>> continentalcrust.add_melt_modifier(crust_solidus, crust_liquidus,
...                                    latentHeatFusion=250.0 * u.kilojoules / u.kilogram / u.kelvin,
...                                    meltFraction=0.,
...                                    meltFractionLimit=0.3,
...                                    meltExpansion=0.13,
...                                    viscosityChangeX1 = 0.15,
...                                    viscosityChangeX2 = 0.30,
...                                    viscosityChange = 1e-3
...                                   )
...

Mechanical Boundary Conditions¶

Mechanical boundary conditions are a critical part of any geodynamic model design. In what follows, we quickly detail the options available for defining the mechanical boundary conditions in Underworld using the UWGeodynamics module.

Questions like how to define boundary conditions and to make sure that those are consistent are beyond the scope of this manual.

We will define a simple model for the sake of the example.

>>> from underworld import UWGeodynamics as GEO

>>> u = GEO.u

>>> Model = GEO.Model(elementRes=(64, 64),
...                   minCoord=(0. * u.kilometre, 0. * u.kilometre),
...                   maxCoord=(64. * u.kilometre, 64. * u.kilometre))

Kinematic boundary conditions¶

Kinematic boundary conditions are set using the set_velocityBCs method. Conditions are defined for each wall (left, right, bottom, top, back and front (3D only)). For each wall, the user must define the condition for each degree of freedom (2 in 2D (x,y), 3 in 3D (x,y,z).

if $V$ is a vector $(V_x, V_y, V_z)$ that we want to apply on the left wall, the left parametre must be defined as left=[Vx, Vy, Vz].

In the following example we set the boundary condition to be:

left wall: $V_x = -1.0 \text{cm / yr}$, $Vy = None$
right wall: $V_x = 1.0 \text{cm / yr}$, $Vy=None$
bottom wall: $V_x = None$, $V_y= 0.$ (free slip)

It is an extension model with a total rate of extension equal to 2.0 centimetre / year. No $V_x$ is prescribed at the bottom, while $V_y$ is set to $0.$ no material will be able to enter or leave the model domain from that side. The material is free to move vertically along the side walls.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_velocityBCs(left=[1.0*u.centimetre/u.year, None],
...                       right=[-1.0*u.centimetre/u.year, None],
...                       bottom=[None, 0.],
...                       top=[None,0.])
...

3D¶

Defining boundary conditions for a 3D model is no different than above. The user must define the velocity components with 3 degree of freedom instead of 2.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model(elementRes=(16, 16, 16),
...                   minCoord=(0. * u.kilometre, 0. * u.kilometre, 0. * u.kilometre),
...                   maxCoord=(64. * u.kilometre, 64. * u.kilometre, 64. * u.kilometre))

>>> Model.set_velocityBCs(left=[1.0*u.centimetre/u.year, None, 0.],
...                        right=[-1.0*u.centimetre/u.year, None, 0.],
...                        bottom=[None, None, 0.],
...                        top=[None, None, 0.],
...                        front=[None, 0., None],
...                        back=[None, 0., None])
...

Velocity varying along a wall¶

At times it is necessary to define a velocity only for a section of a wall and or varying velocities along that wall.

An Underworld function can be passed as a condition.

As an example, we will apply a velocity of $5.0\text{cm/yr}$ for the part of the left wall below 32 kilometre. Velocity is set to be $1.0\text{cm/yr}$ above.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> conditions = [(Model.y < GEO.nd(32 * u.kilometre), GEO.nd(5.0 * u.centimetre/u.year)),
                  (True, GEO.nd(1.0*u.centimetre/u.year))]

>>> function = GEO.uw.fn.branching.conditional(conditions)

>>> Model.set_velocityBCs(left=[function, None],
...                       right=[-1.0*u.centimetre/u.year, None],
...                       bottom=[None, 10.*u.megapascal],
...                       top=[None,0.])
...

Stress Conditions¶

Stress conditions can be applied to the boundaries using the set_stressBCs method:

In the following example we apply a stress of 200.0 megapascal to the bottom of our model:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_stressBCs(bottom=[None, 200. * u.megapascal])
...

Note that you will have to make sure that kinematic and stress conditions are compatible.

Frictional Boundaries¶

Frictional Boundaries can be set as follow:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_frictional_boundary(left=frictionCoeff,
...                               right=frictionCoeff,
...                               bottom=frictionCoeff,
...                               top=frictionCoeff,
...                               thickness=3)
...

Where the left, right, top, bottom parametres indicate the side to which you apply a frictional boundary condition on. frictionCoeff is the friction coefficient (tangent of the friction angle in radians). thickness is the thickness of the boundary in number of elements.

Isostasy¶

Isostasy is an important concept in geodynamics. It is essentially a consequence of the redistribution of mass within a deforming Earth. One important limitation of our geodynamic model is that we model special cases inside rectangular boxes while earth is actually a sphere. One may however need to provide a way to maintain the volume / mass inside the domain in order to mimic isostasy. There is no ideal way to model isostasy in a boxed model, it is however possible to approach isostasy using a support condition.

Options are to:

Balance flows using a kinematic condition at the base of the model.
Balance flows using a stress condition at the base of the model.
Balance flows along the sides.

Lecode Isostasy (kinematic)¶

The Lecode Isostasy submodule provides a way to model isostatic support at the base of the model. It calculates the velocity to apply at the base of each elemental column. It applies the principles of Airy isostatic model by approximating the weight of each column. The calculation is done dynamically and velocities will change from one step to the next. It is a good option to use in most cases.

The option can be used by creating a LecodeIsostasy object using the GEO.LecodeIsostasy class. The object requires the index of the material of reference (the material number). One can apply an average velocity (calculated across each column base) using the average parameter (default to False).

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_velocityBCs(left=[1.0*u.centimetre/u.year, None],
...                       right=[-1.0*u.centimetre/u.year, None],
...                       bottom=[None, GEO.LecodeIsostasy(reference_mat=mantle)],
...                       top=[None,0.])
...

Traction Condition (stress)¶

Another approach to model isostasy is to defined a certain stress at the base of the model. This is done using units of stress (derived SI units = pascal). The model will then maintain the denfined stress by adjusting the flow across the border/boundary.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_stressBCs(bottom=[None, 10.*u.megapascal])
...

Lithostatic Pressure Condition (stress)¶

The lithostatic pressure field can be passed as a boundary condition (stress)

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_stressBCs(left=[-Model.lithostatic_pressureField, None])
...

Thermal Boundary Conditions¶

Absolute temperatures¶

Setting the temperature at the top of a model to be $500 \text{kelvin}$ at the top and $1600 \text{kelvin}$ at the bottom is done as follow.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_temperatureBCs(top=500. * u.degK, bottom=1600. * u.degK)
...

You can of course define temperatures on the sidewalls:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_temperatureBCs(right=500. * u.degK, left=1600. * u.degK)
...

Fix the temperature of a Material

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> Model.set_temperatureBCs(top=500. * u.degK,
...                          bottom=-0.022 * u.milliwatt / u.metre**2,
...                          bottom_material=Model,
...                          materials=[(air, 273. * u.Kelvin)])
...

Note

Model inflow is negative, outflow is positive.

Fix the temperature of internal nodes

You can assign a temperature to a list of nodes by passing a list of node indices (global).

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> Model = GEO.Model()

>>> nodes = [0, 1, 2]
>>> Model.set_temperatureBCs(top=500. * u.degK,
...                          bottom=-0.022 * u.milliwatt / u.metre**2,
...                          bottom_material=Model,
...                          nodeSets=[(273. * u.Kelvin, nodes)])
...

Heat flux¶

Heat Flux can be assign as follow:

>>> from underworld import UWGeodynamics as GEO

>>> u = GEO.u

>>> Model = GEO.Model()
>>> Material = Model.add_material(shape=GEO.Layer(top=Model.top,
...                                               bottom=Model.bottom)
>>> Model.set_heatFlowBCs(bottom=(-0.22 * u.milliwatt / u.metre**2,
...                               Material))
...

Model initialization¶

Initialization of the pressure and temperature fields is done by using the

Model.init_model method.

The default behavior is to not initialise the pressure nor the temperature fields.

You can initialise the fields by passing an Underworld function or a Mesh variable. You can initialise the temperature field to steady-state using the “steady-state” value. Yuo can initialise the pressure field to be lithostatic using the “lithostatic” value.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u

>>> Model = GEO.Model()
>>> Model.density = 2000. * u.kilogram / u.metre**3
>>> Model.init_model(temperature="steady-state", pressure="lithostatic")
...

Warning

The lithostatic pressure calculation relies on a regular quadratic mesh. Most of the time this is fine for model initialization as models often starts on a regular mesh. However, this will not work on a deformed mesh

Running the Model¶

Once your model is set up and initialized. You can run it using the Model.run_for method.

You have 2 options:

Run the model for some given number of steps:

Model.run_for(nstep=10)

Specify an endTime

Model.run_for(endTime=1.0* u.megayears)

which is equivalent to

Model.run_for(1.0*u.megayears)

Specify a timestep¶

UWGeodynamics calculates the time step automatically based on some numerical stability criteria. You can force a specific time step or force the time step to be constant throughout

Model.run_for(1.0*u.megayears, dt=10000. * u.years)

Saving data¶

As your model is running you will need to save the results to files.

The Model.run_for method provides a series of arguments to help you save the results at some regular intervals and/or specified times. You can define:

A checkpoint_interval

Model.run_for(endTime=1.0*u.megayears,
              checkpoint_interval=0.1* u.megayears)

The value passed to the checkpoint_interval must have units of time 1. A list of checkpoint times:

Model.run_for(endTime=1.0*u.megayears,
                 checkpoint_interval=0.1* u.megayears,
                 checkpoint_times=[0.85 * u.megayears,
                                   0.21 * u.megayears])

This can be used together or without the checkpoint_interval

UWGeodynamics will save all the fields defined in the GEO.rcParams[“default.outputs”] list. You can change that list before running the model.

Checkpointing¶

By checkpointing we mean saving the data required to restart the Model. This includes the mesh, the swarm and all the associated variables.

However, as the swarm and the swarm variables can be very large and can take a lot of space on disk, the user can decide to save them only every second, third, fourth etc. checkpoint step.

This is done passing the restart_checkpoint argument to the Model.run_for function:

Model.run_for(endTime=1.0*u.megayears,
              checkpoint_interval=0.1* u.megayears,
              restart_checkpoint=2

By default, the swarm and the swarm variables are saved every time the model reaches a checkpoint time (restart_checkpoint=1).

Pre / Post-solve hook functions¶

We provide 2 access points for injection of custom functions.

def my_functionA():
    # do something
    print("Hello, I am running a pre-solve function")
    return

def my_functionB():
    # do something
    print("Hello, I am running a post-solve function")
    return

Model.pre_solve_function["A"] = my_functionA
Model.post_solve_function["B"] = my_functionB

Note that the functions are executed in the order they were defined.

Solver Callback functions¶

User can provide custom callback functions to the solver itself. The function(s) will be executed after each solve. This gives the possibility to tweak the behaviour of the non-linear iterations loop.

def my_function()
    # do something
    print("Hello, This is a solver callback")

Model.callback_function["my_function"] = my_function

Restarting the Model¶

When checkpointing a model only the mesh, swarms their associates variables are explicitely saved. Since the model state is not explicitly saved, thus the user needs to recreate the Model object before restarting it. In practice, this means the user must run all commands preceding the Model.run_for command.

The user can then restart a model using the restartStep and restartDir arguments:

restartStep is None by default. The step numbercyou want to restart from. If -1, restarts from the last available step in restartDir
restartDir is the folder where the program should look for previously saved data or checkpoints. It is set to Model.outputs by default.

from underworld import UWGeodynamics as GEO

u = GEO.u

Model = GEO.Model(elementRes=(64, 64),
                  minCoord=(0. * u.kilometre, 0. * u.kilometre),
                  maxCoord=(64. * u.kilometre, 64. * u.kilometre))

# Default (restart, restartDir are optional in this case)
Model.run_for(2.0 * u.megayears, restartStep=-1, restartDir="your_restart_directory")

# Restart from step 10
Model.run_for(2.0 * u.megayears, restartStep=10, restartDir="your_restart_directory")

# Overwrite existing outputs
Model.run_for(2.0 * u.megayears, restartStep=False)

Model outputs¶

All mesh variables / fields defined in the GEO.rcParams["default.outputs"] are saved as HDF5 files to the outputDir directory at every output times. An XMF file is provided and can be used to open the files in Paraview

All variables required for a restart are saved as HDF5 files to the outputDir directory at each checkpoint time. An XMF file is also provided.

Passive Tracers and tracked fields are also saved as HDF5 files at every output and checkpoint times. Each of then has an associated XMF file.

Parallel run¶

A Model can be run on multiple processors:

You first need to convert your jupyter notebook to a python script:

jupyter nbconvert --to python my_script.ipynb

You can then run the python script as follow:

mpirun -np 4 python my_script.py

Warning

Underworld and UWGeodynamics functions are parallel safe and can be run on multiple CPUs. This might not be the case with other python libraries you might be interested in using with your Model. For example, matplotlib plots will not work in parallel and must be processed in serial. Tutorial 1 has examples of matplotlib plots which are only done on the rank 0 CPU.

Passive Tracers¶

>>> from underworld import UWGeodynamics as GEO
>>> import numpy as np

>>> u = GEO.u

>>> Model = GEO.Model()
>>> npoints = 1000
>>> coords = np.ndarray((npoints, 2))
>>> coords[:, 0] = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), npoints)
>>> coords[:, 1] = GEO.nd(32. * u.kilometre)
>>> Model.add_passive_tracers(vertices=coords)

You can pass a list of centroids to the Model.add_passive_tracers method. In that case, the coordinates of the passive tracers are relative to the position of the centroids. The pattern is repeated around each centroid.

>>> from underworld import UWGeodynamics as GEO
>>> import numpy as np

>>> u = GEO.u
>>> Model = GEO.Model()
>>> cxpos = np.linspace(GEO.nd(20*u.kilometer), GEO.nd(40*u.kilometer),5)
>>> cypos = np.linspace(GEO.nd(20*u.kilometer), GEO.nd(40*u.kilometer),5)
>>> cxpos, cypos = np.meshgrid(cxpos, cypos)
>>> coords_centroid = np.ndarray((cxpos.size, 2))
>>> coords_centroid[:, 0] = cxpos.ravel()
>>> coords_centroid[:, 1] = cypos.ravel()
>>>
>>> coords = np.zeros((1, 2))
>>> Model.add_passive_tracers(vertices=coords, centroids=coords_centroid)

We provide a function to create circles on a grid:

>>> from underworld import UWGeodynamics as GEO

>>> coords = GEO.circles_grid(radius = 2.0 * u.kilometer,
...                           minCoord=[Model.minCoord[0], lowercrust.bottom],
...                           maxCoord=[Model.maxCoord[0], 0.*u.kilometer])

Tracking Values¶

Passive tracers can be used to track values of fields at specific location through time. Tracking projected fields (fields with the prefix proj) is discouraged as these fields can’t be restarted after checkpointing. Instead, one can track the non projected field, i.e. projViscosityField, viscosityField)

>>> from underworld import UWGeodynamics as GEO
>>> import numpy as np

>>> u = GEO.u
>>> Model = GEO.Model()

>>> npoints = 1000
>>> coords = np.ndarray((npoints, 2))
>>> coords[:, 0] = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), npoints)
>>> coords[:, 1] = GEO.nd(32. * u.kilometre)
>>> Model.add_passive_tracers(name="p1", vertices=coords)

>>> Model.p1_tracers.add_tracked_field(Model.pressureField,
                                       name="tracers_press",
                                       units=u.megapascal,
                                       dataType="float")
>>> Model.p1_tracers.add_tracked_field(Model.strainRateField,
                                       name="tracers_strainRate",
                                       units=1.0/u.second,
                                       dataType="float")

Surface Processes¶

A range of basic surface processes function are available from the surfaceProcesses submodule. Surface processes are turned on once you have passed a valid surface processes function to the surfaceProcesses method of the Model object.

Example:

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> air = GEO.Material()
>>> sediment = GEO.Material()
>>> Model.surfaceProcesses = GEO.surfaceProcesses.SedimentationThreshold(
...     air=[air], sediment=[sediment], threshold=0. * u.metre)

Three simple function are available:

Total Erosion Above Threshold (ErosionThreshold).
Total Sedimentation Below Threshold (SedimentationThreshold)
Combination of the 2 above. (ErosionAndSedimentationThreshold)

Erosion and sedimentation rate¶

Adds an erosion and sedimentation rate to the surface. A pre-defined vertical co-ordinate (surfaceElevation) needs to be defined to stop erodion below that level and sedimenation above it.

Example:

  >>> from underworld import UWGeodynamics as GEO
  >>> u = GEO.u
  >>> air = GEO.Material()
  >>> sediment = GEO.Material()
  >>> Model.surfaceProcesses = GEO.surfaceProcesses.velocitySurface_2D(
  ...     airIndex=air.index, sedimentIndex=sediment.index,
  ...     sedimentationRate= 2.*u.millimeter / u.year, erosionRate= 2.*u.millimeter / u.year,
  ...     surfaceElevation=0.*u.kilometer,
  ...     surfaceArray = coords)

Diffusive surface
~~~~~~~~~~~~~~~~~~~~~~
Adds a linear diffusive surface to the model.

Example:

.. code:: python

   >>> from underworld import UWGeodynamics as GEO
   >>> u = GEO.u
   >>> air = GEO.Material()
   >>> sediment = GEO.Material()
   >>> Model.surfaceProcesses = GEO.surfaceProcesses.diffusiveSurface_2D(
   ...     airIndex=air.index, sedimentIndex=sediment.index,
   ...     D= 1000.0*u.meter**2/u.year,
   ...     surfaceArray = coords)

Coupling with Badlands¶

UWGeodynamics provides a way to couple an Underworld model to Badlands.

>>> from underworld import UWGeodynamics as GEO
>>> u = GEO.u
>>> air = GEO.Material()
>>> sediment = GEO.Material()
>>> Model.surfaceProcesses = GEO.surfaceProcesses.Badlands(
...     airIndex=[air.index], sedimentIndex=sediment.index,
...     XML="ressources/badlands.xml", resolution=1. * u.kilometre,
...     checkpoint_interval=0.01 * u.megayears)

This will allow communication between the UWGeodynamics model and the Badlands surface processes model. Badlands input parameters must be defined inside an XML file as described in the module documentation. We provide an XML example. The resulting Model is a 2-way coupled thermo-mechanical model with surface processes, where the velocity field retrieved from the thermo-mechanical model is used to advect the surface in the Surface Processes Model. The surface in subjected to erosion and deposition. The distribution of materials in the thermomechanical model is then updated.

Users must define a list of material describing the air layers (usually, air and sticky air). It is also require to define an UWGeodynamics.Material object describing the sediment that will be deposited. The index of the Material is passed to the surfaceProcesses function. Users can also provide an Underworld function returning an index of an existing UWGeodynamics.Material.

It is recommended to use a higher spatial resolution in the surface processes model than in the thermo-mechanical model.

Note

When the Thermomechanical model is 2D, the velocity field at the surface is extrapolated in the 3D dimension and the resulting model is a T or 2.5D model (symmetric regional uplift). If the thermomechanical model is 3D the coupling is done in 3D.

Deforming Mesh¶

Uniaxial deformation can be turned on using the Model.mesh_advector() method. The method takes an axis argument which defines the direction of deformation (x=0, y=1, z=2)

>>> Model.mesh_advector(axis=0)

Element are stretched or compressed uniformly across the model. This will result in a change in resolution with time.

Top Free surface¶

Free surface can be turned on using the Model.freesurface switch.

>>> Model.freesurface = True

Warning

No stabilization algorithm has been implemented yet.

Dynamic rc settings¶

You can dynamically change the default rc settings in a python script or interactively from the python shell. All of the rc settings are stored in a dictionary-like variable called UWGeodynamics.rcParams, which is global to the UWGeodynamics package. rcParams can be modified directly, for example:

>>> from underworld import UWGeodynamics as GEO
>>> GEO.rcParams['solver'] = "mumps"
>>> GEO.rcParams['penalty'] = 1e6

The UWGeodynamics.rcdefaults command will restore the standard UWGeodynamics default settings.

There is some degree of validation when setting the values of rcParams, see UWGeodynamics.rcsetup for details.

rcParams¶
Name	Function	Default value
CFL	Set CFL Factor	0.5
solver	Set Solver	“mg” (multigrid), options are “mumps”, “lu”
penalty	Set penalty value	0.0 or None
initial.nonlinear.tolerance	Set nonlinear tolerance for Stokes first solve	1e-2
nonlinear.tolerance	Set nonlinear tolerance for solves	1e-2
initial.nonlinear.min.iterations	Set minimal number of Picard iterations (first solve)	2
initial.nonlinear.max.iterations	Set maximal number of Picard iterations (first solve)	500
nonlinear.min.iterations	Set minimal number of Picard iterations	2
nonlinear.max.iterations	Set maximal number of Picard iterations	500
default.outputs	List of fields to be saved at checkpoint	[“temperature”, “pressureField”, “strainRateField”, “velocityField”, “projStressField”, “projTimeField”, “projMaterialField”, “projViscosityField”, “projPlasticStrain”, “projDensityField”]
swarm.particles.per.cell.2D	Initial number of particles per cell for 2D models	40
swarm.particles.per.cell.3D	Initial number of particles per cell for 3D models	120
popcontrol.particles.per.cell.2D	Minimum number of particles per cell	40
popcontrol.particles.per.cell.3D	Maximum number of particles per cell	120
popcontrol.aggressive	Turn on Aggressive population control	True
popcontrol.split.threshold	Population control split threshold	0.15
popcontrol.max.splits	Population control maximum number of splits	10
shear.heating	Turn shear heating on / off	False
surface.pressure.normalization	Turn surface pressure normalization on / off	True
pressure.smoothing	Turn pressure smoothing after solve on / off	True
advection.diffusion.method	Advection Diffusion solve method	“SUPG”
rheologies.combine.method	Visco-plastic rheology combination	“Minimum”, options are “Minimum”, “Harmonic”
averaging.method	Multi-material element averaging method	“arithmetic” options are “arithmetic”, “geometric”, “harmonic”, “maximum”, “minimum”, “root mean square”
time.SIunits	Default output units for time field	u.year
viscosityField.SIunits	Default output units for viscosity field	u.pascal * u.second
densityField.SIunits	Default output units for density field	u.kilogram / u.metre**3
velocityField.SIunits	Default output units for velocity field	u.centimetre / u.year
temperature.SIunits	Default output units for temperature field	u.degK
pressureField.SIunits	Default output units for pressure field	u.pascal
strainRateFieldSIunits	Default output units for strain rate field	u.pascal
projStressTensor.SIunits	Default output units for mesh projected stress tensor field	u.pascal
projStressField.SIunits	Default output units for mesh projected stress field	u.pascal
projViscosityFIeld.SIunits	Default output units for mesh projected viscosities	u.pascal * u.second
projTimeField.SIunits	Default output units for mesh projected time field.	u.year

The `uwgeodynamicsrc` file¶

UWGeodynamics uses uwgeodynamicsrc configuration files to customize all kinds of properties, which we call rc settings or rc parameters. For now, you can control the defaults of a limited set of properties.

UWGeodynamics looks for uwgeodynamicsrc in four locations, in the following order:

uwgeodynamicsrc in the current working directory, usually used for specific customizations for a particular model setup that you do not want to apply elsewhere.
$UWGEODYNAMICSRC if it is a file, else $UWGEODYNAMICSRC/uwgeodynamicsrc.
It next looks in a user-specific place, depending on your platform:
- On Linux, it looks in .config/uwgeodynamics/uwgeodynamicsrc (or $XDG_CONFIG_HOME/uwgeodynamics/uwgeodynamicsrc) if you’ve customized your environment.
- On other platforms, it looks in .uwgeodynamics/uwgeodynamicsrc.
{INSTALL}/UWGeodynamics/uwgeo-data/uwgeodynamicsrc, where {INSTALL} is something like /usr/lib/python2.7/site-packages on Linux, and maybe C:\\Python27\\Lib\\site-packages on Windows. Every time you install UWgeodynamics, this file will be overwritten, so if you want your customizations to be saved, please move this file to your user-specific directory.

To display where the currently active uwgeodynamicsrc file was loaded from, one can do the following:

>>> from underworld import UWGeodynamics as GEO
>>> GEO.uwgeodynamics_fname()
'/workspace/user_data/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc'

See below for a sample.

UWGeodynamics User Guide¶

The Jupyter notebook¶

Where to find documentation?¶

Design principles¶

import UWGeodynamics¶

Visualization¶

Simple examples:¶

Working with units¶

Model Scaling¶

dimensionalise / non-dimensionalise¶

The Model object¶

The Material object¶

Predefined Material objects¶

User defined¶

Material Attributes¶

Global properties¶

Material shape¶

Temperature and Pressure dependent densities¶

Rheologies¶

Newtonian Rheology¶

Non-Newtonian Rheology¶

Composite Viscosity¶

Plastic Behavior (Yield)¶

Elasticity¶

Simple phase change¶

Melt¶

Mechanical Boundary Conditions¶

Kinematic boundary conditions¶

3D¶

Velocity varying along a wall¶

Stress Conditions¶

Frictional Boundaries¶

Isostasy¶

Lecode Isostasy (kinematic)¶

Traction Condition (stress)¶

Lithostatic Pressure Condition (stress)¶

Thermal Boundary Conditions¶

Absolute temperatures¶

Heat flux¶

Model initialization¶

Running the Model¶

Specify a timestep¶

Saving data¶

Checkpointing¶

Pre / Post-solve hook functions¶

Solver Callback functions¶

Restarting the Model¶

Model outputs¶

Parallel run¶

Passive Tracers¶

Tracking Values¶

Surface Processes¶

Erosion and sedimentation rate¶

Coupling with Badlands¶

Deforming Mesh¶

Top Free surface¶

Dynamic rc settings¶

The uwgeodynamicsrc file¶

_uwgeodynamicsrc-sample:¶

The `uwgeodynamicsrc` file¶