Edit on GitHub

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]
source
GeophysicalModelGenerator.AddBox! — Function
AddBox!(Phase, Temp, Grid::LaMEM_grid; 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 - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
  • 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 = CartData(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 = CartData(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"
source
GeophysicalModelGenerator.HalfspaceCoolingTemp — Type
HalfspaceCooling(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]
source
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]
source
GeophysicalModelGenerator.LithosphericPhases — Type
LithosphericPhases(Layers=[10 20 30], 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 the layers from top downwards
  • Phases : the phases of the layers from top down. Note that this array
  • Tlab : Temperature of the lithosphere asthenosphere boundary. If specified, the phases at locations with T>Tlab is set to Phases[end]
source
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
source
GeophysicalModelGenerator.ReadData_PVTR — Function
Data::CartData = 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/")
CartData 
  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)
source