Example: simple model, linear viscous rheology

This example shows how to build a very simple model using default linear viscous rheology and standard Dirichlet type boundary conditions for a velocity field imposing shortening in \(z\) direction as shown by the figure below. Because the viscosity chosen in this example does not depend on temperature, the thermal part of the model is not included.

../_images/compression.png

1. Create a domain

We define a 3D domain \(\Omega = [0,600]\times[-250,0]\times[0,300]\) km3 \(\in \mathbb R^3\) discretized by a regular grid of 9x9x9 nodes.

import os
import numpy as np
import genepy as gp

# 3D domain
dimensions = 3
O = np.array([0,-250e3,0],    dtype=np.float64) # Origin
L = np.array([600e3,0,300e3], dtype=np.float64) # Length
n = np.array([9,9,9],         dtype=np.int32)   # Number of Q1 nodes i.e. elements + 1
# Create Domain class instance
Domain = gp.Domain(dimensions,O,L,n)

2. Velocity function

We define a simple orthogonal shortening velocity in the \(z\) direction.

# velocity
cma2ms  = 1e-2 / (3600.0 * 24.0 * 365.0) # cm/a to m/s conversion
u_norm  = 1.0 * cma2ms                   # horizontal velocity norm
u_dir   = "z"                            # direction in which velocity varies
u_type  = "compression"                  # extension or compression
# Create Velocity class instance
BCs = gp.VelocityLinear(Domain,u_norm,u_dir,u_type)

# Access the symbolic velocity function
u = BCs.u

Note

In this example, the derivatives of the velocity are not used.

3. Initial conditions

In this example we do not impose any initial plastic strain value nor mesh refinement. Therefore the initial conditions are only the Domain and the velocity function. They will be used to generate the options for pTatin3d model.

# Initial conditions
model_ics = gp.InitialConditions(Domain,u)

4. Boundary conditions

Because the imposed velocity is orthogonal to the boundary we can define the velocity boundary conditions using Dirichlet type boundary conditions.

Note

In the following example a path to the mesh files describing the boundaries is provided. These mesh files are located in "ptatin-gene/src/models/gene3d/examples". You can modify the root variable to match the location of the mesh files on your system or remove that part of the code if you do not have access to these files. Note however that pTatin3d requires mesh files to define the boundaries.

Details on the methods used to define the boundary conditions can be found in the boundary conditions section.

# boundary conditions
# path to mesh files (system dependent, change accordingly)
root = os.path.join(os.environ['PTATIN'],"ptatin-gene/src/models/gene3d/examples")
# Velocity boundary conditions
u_bcs = [
  gp.Dirichlet(23,"Zmax",["z"],u, mesh_file=os.path.join(root,"box_ptatin_facet_23_mesh.bin")), # orthogonal shortening
  gp.Dirichlet(37,"Zmin",["z"],u, mesh_file=os.path.join(root,"box_ptatin_facet_37_mesh.bin")), # orthogonal shortening
  gp.Dirichlet(32,"Xmax",["x"],u, mesh_file=os.path.join(root,"box_ptatin_facet_23_mesh.bin")), # free-slip
  gp.Dirichlet(14,"Xmin",["x"],u, mesh_file=os.path.join(root,"box_ptatin_facet_37_mesh.bin")), # free-slip
  gp.DirichletUdotN(33,"Bottom",  mesh_file=os.path.join(root,"box_ptatin_facet_33_mesh.bin")), # basal outflow
]
# collect all boundary conditions
model_bcs = gp.ModelBCs(u_bcs)

5. Material parameters

Next we define the material properties of each Region and gather them all in a ModelRegions class instance. In this example we use the default values for all regions:

regions = [
  # Upper crust
  gp.Region(38),
  # Lower crust
  gp.Region(39),
  # Lithosphere mantle
  gp.Region(40),
  # Asthenosphere
  gp.Region(41)
]
model_regions = gp.ModelRegions(regions,
                                mesh_file=os.path.join(root,"box_ptatin_md.bin"),
                                region_file=os.path.join(root,"box_ptatin_region_cell.bin"))

6. Create the model and generate options

Finally, we create the model by gathering all the information defined previously and we save the options to a file named simple_shortening_model.opts.

# create class instance
model = gp.Model(model_ics,model_regions,model_bcs)
# write the options for ptatin3d
with open("simple_shortening_model.opts","w") as f:
  f.write(model.options)