LaMEM

In order to generate geodynamic simulations from setups created with GeophysicalModelGenerator.jl, we provide a few routines that directly create marker input files for the 3D geodynamic modelling software LaMEM, which is an open-source cartesian code that is well-suited to perform crustal and lithospheric-scale simulations. If you want to learn how to run LaMEM simulations, please have a look at the wiki page.

The routines provided here have the following functionality:

Read LaMEM *.dat files (to get the size of the domain)
Read LaMEM processor partitioning file
Add lithospheric boxes to a setup, that may have a layered structure and various thermal structures
Save LaMEM marker files in serial or in parallel
Read a LaMEM timestep

GeophysicalModelGenerator.ReadLaMEM_InputFile — Function

Grid::LaMEM_grid = ReadLaMEM_InputFile(file)

Parses a LaMEM input file and stores grid information in the Grid structure.

Example

julia> Grid = ReadLaMEM_InputFile("SaltModels.dat") 
LaMEM Grid: 
nel         : (32, 32, 32)
marker/cell : (3, 3, 3)
markers     : (96, 96, 96)
x           ϵ [-3.0 : 3.0]
y           ϵ [-2.0 : 2.0]
z           ϵ [-2.0 : 0.0]

GeophysicalModelGenerator.GetProcessorPartitioning — Function

nProcX,nProcY,nProcZ, xc,yc,zc, nNodeX,nNodeY,nNodeZ = GetProcessorPartitioning(filename)

Reads a LaMEM processor partitioning file, used to create marker files, and returns the parallel layout

GeophysicalModelGenerator.Save_LaMEMTopography — Function

Save_LaMEMTopography(Topo::CartData, filename::String)

This writes a topography file Topo for use in LaMEM, which should have size (nx,ny,1) and contain the field :Topography

GeophysicalModelGenerator.Save_LaMEMMarkersParallel — Function

Save_LaMEMMarkersParallel(Grid::CartData; PartitioningFile=empty, directory="./markers", verbose=true)

Saves a LaMEM marker file from the CartData structure Grid. It must have a field called Phases, holding phase information (as integers) and optionally a field Temp with temperature info. It is possible to provide a LaMEM partitioning file PartitioningFile. If not, output is assumed to be for one processor.

The size of Grid should be consistent with what is provided in the LaMEM input file. In practice, the size of the mesh can be retrieved from a LaMEM input file using ReadLaMEM_InputFile.

Example

julia> Grid    = ReadLaMEM_InputFile("LaMEM_input_file.dat")
julia> Phases  = zeros(Int32,size(Grid.X));
julia> Temp    = ones(Float64,size(Grid.X));
julia> Model3D = CartData(Grid, (Phases=Phases,Temp=Temp))
julia> Save_LaMEMMarkersParallel(Model3D)
Writing LaMEM marker file -> ./markers/mdb.00000000.dat

If you want to create a LaMEM input file for multiple processors:

julia> Save_LaMEMMarkersParallel(Model3D, PartitioningFile="ProcessorPartitioning_4cpu_1.2.2.bin")
Writing LaMEM marker file -> ./markers/mdb.00000000.dat
Writing LaMEM marker file -> ./markers/mdb.00000001.dat
Writing LaMEM marker file -> ./markers/mdb.00000002.dat
Writing LaMEM marker file -> ./markers/mdb.00000003.dat

GeophysicalModelGenerator.ReadData_PVTR — Function

Data::ParaviewData = ReadData_PVTR(fname, dir)

Reads a parallel, rectilinear, *.vts file with the name fname and located in dir and create a 3D Data struct from it.

Example

julia> Data = ReadData_PVTR("Haaksbergen.pvtr", "./Timestep_00000005_3.35780500e-01/")
ParaviewData 
  size  : (33, 33, 33)
  x     ϵ [ -3.0 : 3.0]
  y     ϵ [ -2.0 : 2.0]
  z     ϵ [ -2.0 : 0.0]
  fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res)

GeophysicalModelGenerator.AddBox! — Function

AddBox!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1),
        T=nothing )

Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - grid structure (usually obtained with ReadLaMEM_InputFile, but can also be other grid types)
xlim - left/right coordinates of box
ylim - front/back coordinates of box [optional; if not specified we use the whole box]
zlim - bottom/top coordinates of box
Origin - the origin, used to rotate the box around. Default is the left-front-top corner
StrikeAngle - strike angle of slab
DipAngle - dip angle of slab
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Examples

Example 1) Box with constant phase and temperature & a dip angle of 10 degrees:

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid: 
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview 
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

Example 2) Box with halfspace cooling profile

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview 
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

GeophysicalModelGenerator.AddSphere! — Function

AddSphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, radius=Tuple{1},
        phase = ConstantPhase(1).
        T=nothing )

Adds a sphere with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
cen - center coordinates of sphere
radius - radius of sphere
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Example

Sphere with constant phase and temperature:

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid: 
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddSphere!(Phases,Temp,Grid, cen=(0,0,-1), radius=0.5, phase=ConstantPhase(2), T=ConstantTemp(800))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview 
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

GeophysicalModelGenerator.AddEllipsoid! — Function

AddEllipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, axes=Tuple{3},
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1).
        T=nothing )

Adds an Ellipsoid with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
cen - center coordinates of sphere
axes - semi-axes of ellipsoid in X,Y,Z
Origin - the origin, used to rotate the box around. Default is the left-front-top corner
StrikeAngle - strike angle of slab
DipAngle - dip angle of slab
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Example

Ellipsoid with constant phase and temperature, rotated 90 degrees and tilted by 45 degrees:

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid: 
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddEllipsoid!(Phases,Temp,Grid, cen=(-1,-1,-1), axes=(0.2,0.1,0.5), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(3), T=ConstantTemp(600))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview 
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

GeophysicalModelGenerator.AddCylinder! — Function

AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},
        phase = ConstantPhase(1).
        T=nothing )

Adds a cylinder with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

Phase - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - Grid structure (usually obtained with ReadLaMEM_InputFile)
base - center coordinate of bottom of cylinder
cap - center coordinate of top of cylinder
radius - radius of the cylinder
phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Example

Cylinder with constant phase and temperature:

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid: 
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddCylinder!(Phases,Temp,Grid, base=(-1,-1,-1.5), cap=(1,1,-0.5), radius=0.25, phase=ConstantPhase(4), T=ConstantTemp(400))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview 
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

GeophysicalModelGenerator.makeVolcTopo — Function

makeVolcTopo(Grid::LaMEM_grid; center::Array{Float64, 1}, height::Float64, radius::Float64, crater::Float64, base=0.0m, background=nothing)

Creates a generic volcano topography (cones and truncated cones)

Parameters

Grid - LaMEM grid (created by ReadLaMEM_InputFile)
center - x- and -coordinates of center of volcano
height - height of volcano
radius - radius of volcano

Optional Parameters

crater - this will create a truncated cone and the option defines the radius of the flat top
base - this sets the flat topography around the volcano
background - this allows loading in a topography and only adding the volcano on top (also allows stacking of several cones to get a volcano with different slopes)

Example

Cylinder with constant phase and temperature:

julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid: 
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Topo = makeVolcTopo(Grid, center=[0.0,0.0], height=0.4, radius=1.5, crater=0.5, base=0.1)
CartData 
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.4]
    fields  : (:Topography,)
  attributes: ["note"]
julia> Topo = makeVolcTopo(Grid, center=[0.0,0.0], height=0.8, radius=0.5, crater=0.0, base=0.4, background=Topo.fields.Topography)
CartData 
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.8]
    fields  : (:Topography,)
  attributes: ["note"]
julia> Write_Paraview(Topo,"VolcanoTopo")           # Save topography to paraview 
Saved file: VolcanoTopo.vts

GeophysicalModelGenerator.ConstantTemp — Type

ConstantTemp(T=1000)

Sets a constant temperature inside the box

Parameters

T : the value

GeophysicalModelGenerator.LinearTemp — Type

LinearTemp(Ttop=0, Tbot=1000)

Set a linear temperature structure from top to bottom

Parameters

Ttop : the value @ the top
Tbot : the value @ the bottom

GeophysicalModelGenerator.HalfspaceCoolingTemp — Type

HalfspaceCoolingTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)

Sets a halfspace temperature structure in plate

Parameters

Tsurface : surface temperature [C]
Tmantle : mantle temperature [C]
Age : Thermal Age of plate [Myrs]
Adiabat : Mantle Adiabat [K/km]

GeophysicalModelGenerator.SpreadingRateTemp — Type

SpreadingRateTemp(Tsurface=0, Tmantle=1350, Adiabat=0, MORside="left",SpreadingVel=3, AgeRidge=0, maxAge=80)

Sets a halfspace temperature structure within the box, combined with a spreading rate (which implies that the plate age varies)

Parameters

Tsurface : surface temperature [C]
Tmantle : mantle temperature [C]
Adiabat : Mantle Adiabat [K/km]
MORside : side of the box where the MOR is located ["left","right","front","back"]
SpreadingVel : spreading velocity [cm/yr]
AgeRidge : thermal age of the ridge [Myrs]
maxAge : maximum thermal Age of plate [Myrs]

GeophysicalModelGenerator.ConstantPhase — Type

ConstantPhase(phase=1)

Sets a constant phase inside the box

Parameters

phase : the value

GeophysicalModelGenerator.Compute_Phase — Function

Phase = Compute_Phase(Phase, Temp, X, Y, Z, s::LithosphericPhases, Ztop)

or

Phase = Compute_Phase(Phase, Temp, Grid::AbstractGeneralGrid, s::LithosphericPhases)

This copies the layered lithosphere onto the Phase matrix.

Parameters

Phase - Phase array
Temp - Temperature array
X - x-coordinate array (consistent with Phase and Temp)
Y - y-coordinate array (consistent with Phase and Temp)
Z - Vertical coordinate array (consistent with Phase and Temp)
s - LithosphericPhases
Ztop - Vertical coordinate of top of model box
Grid - Grid structure (usually obtained with ReadLaMEM_InputFile)

GeophysicalModelGenerator.LithosphericPhases — Type

LithosphericPhases(Layers=[10 20 15], Phases=[1 2 3 4], Tlab=nothing )

This allows defining a layered lithosphere. Layering is defined from the top downwards.

Parameters

Layers : The thickness of each layer, ordered from top to bottom. The thickness of the last layer does not have to be specified.
Phases : The phases of the layers, ordered from top to bottom.
Tlab : Temperature of the lithosphere asthenosphere boundary. If specified, the phases at locations with T>Tlab are set to Phases[end].

GeophysicalModelGenerator.LaMEM_grid — Type

Structure that holds information about the LaMEM grid (usually read from an input file).

GeophysicalModelGenerator.CreatePartitioningFile — Function

CreatePartitioningFile(LaMEM_input::String, NumProc::Int64; LaMEM_dir::String=pwd(), LaMEM_options::String="", MPI_dir="", verbose=true)

This executes LaMEM for the input file LaMEM_input & creates a parallel partitioning file for NumProc processors. The directory where the LaMEM binary is can be specified; if not it is assumed to be in the current directory. Likewise for the mpiexec directory (if not specified it is assumed to be available on the command line).