- adds a function to create the initial model setup (also in parallel)

- simplified GMG routines
- uses `files` as default
- add plotting routines that are loaded when GLMakie is loaded

src/LaMEM.jl

@@ -2,6 +2,7 @@ module LaMEM
 
 using GeoParams
 using .GeoParams
+using Requires
 export NO_units, GEO_units, SI_units, km, m, Pa, Pas, kg, cm, yr
 #export GeoParams
 
@@ -18,19 +19,28 @@ using .Run
 export run_lamem, run_lamem_save_grid, mpiexec
 export remove_popup_messages_mac, show_paths_LaMEM
 
+# Functions that help running LaMEM directly from julia
 include("LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl")
 using .LaMEM_Model
-export  Model, Write_LaMEM_InputFile, 
+export  Model, Write_LaMEM_InputFile, create_initialsetup,
         Scaling, Grid, Time, FreeSurface, BoundaryConditions, SolutionParams,
         Solver, ModelSetup, 
         geom_Sphere,
         Output,
         Materials, Phase, Softening, PhaseTransition, Dike, 
-        add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!,
+        add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!,
         add_phasetransition!, add_dike!, add_geom!     
 
+
 using .Run.LaMEM_jll
 export LaMEM_jll        # export LaMEM_jll as well & directories
 
 
+function __init__()
+    @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin 
+        @eval include("MakieExt.jl")
+    end
+end
+
+
 end # module

src/LaMEM_ModelGeneration/ErrorChecking.jl

@@ -17,6 +17,12 @@ function Check_LaMEM_Model(m::Model)
         error("If you use internal geometries to set phases, you need to at least specify one internal geometry object. 
                Example: add_geom!(model, geom_Sphere())")
     end
+
+    if (m.ModelSetup.msetup=="files") && diff([extrema(m.Grid.Phases)...])[1]==0 && diff([extrema(m.Grid.Temp)...])[1]==0
+        error("Your initial `Temp` grid is constant, as is your initial `Phases` grid. 
+        You do need to set some variability to see action, for example with the GMG function `AddSphere!(model,cen=(0.0,0.0,0.0), radius=(0.15, ))` ")
+    end
+
 
     return nothing
 end

src/LaMEM_ModelGeneration/GMG_interface.jl

@@ -8,19 +8,13 @@ import GeophysicalModelGenerator: AddBox!, AddSphere!, AddEllipsoid!, AddCylinde
     AddBox!(model::Model; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
             Origin=nothing, StrikeAngle=0, DipAngle=0,
             phase = ConstantPhase(1),
-        T=nothing )
+            T=nothing )
 
 Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models
 See the documentation of the GMG routine
 
 """
-function  AddBox!(model::Model; xlim=Tuple{2}, ylim=Tuple{2}, zlim=Tuple{2}, 
-            Origin=nothing, StrikeAngle=0, DipAngle=0,  phase = ConstantPhase(1), T=nothing )
-
-return AddBox!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; xlim=xlim, ylim=ylim, zlim=zlim, 
-            Origin=Origin, StrikeAngle=StrikeAngle, DipAngle=DipAngle,  phase = phase, T=T )
-
-end
+AddBox!(model::Model; kwargs...) = AddBox!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 
 
 
 """
@@ -29,58 +23,29 @@ end
 See the documentation of the GMG routine
 
 """
-function  AddSphere!(model::Model; cen=Tuple{3}, radius=Tuple{1}, phase = ConstantPhase(1), T=nothing)
-
-    AddSphere!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, cen=cen, radius=radius, phase=phase, T=T)
-
-    return nothing
-end
-
-
+AddSphere!(model::Model; kwargs...) = AddSphere!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 
 
 """
-    AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid;      # required input
-                    base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},       # center and radius of the sphere
-                    phase = ConstantPhase(1),                           # Sets the phase number(s) in the sphere
-                    T=nothing )                                         # Sets the thermal structure (various fucntions are available)
+    AddCylinder!(model::Model;                                      # required input
+                    base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},   # center and radius of the sphere
+                    phase = ConstantPhase(1),                       # Sets the phase number(s) in the sphere
+                    T=nothing )                                     # Sets the thermal structure (various fucntions are available)
 
 
 See the documentation of the GMG routine
 
 """
-function  AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid;      # required input
-                base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},       # center and radius of the sphere
-                phase = ConstantPhase(1),                           # Sets the phase number(s) in the sphere
-                T=nothing )                                         # Sets the thermal structure (various fucntions are available)
-
-
-    AddCylinder!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, 
-                base=base, cap=cap, radius=radius,                  # center and radius of the sphere
-                phase = phase,                                      # Sets the phase number(s) in the sphere
-                T=T )
-
-    return nothing
-end
+AddCylinder!(model::Model; kwargs...) = AddCylinder!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 
+
 
 """
-    AddEllipsoid!(model::Model;      # required input
-                    cen=Tuple{3}, axes=Tuple{3},                           # center and semi-axes of the ellpsoid
-                    Origin=nothing, StrikeAngle=0, DipAngle=0,             # origin & dip/strike
-                    phase = ConstantPhase(1),                              # Sets the phase number(s) in the box
+    AddEllipsoid!(model::Model;                                 # required input
+                    cen=Tuple{3}, axes=Tuple{3},                # center and semi-axes of the ellpsoid
+                    Origin=nothing, StrikeAngle=0, DipAngle=0,  # origin & dip/strike
+                    phase = ConstantPhase(1),                   # Sets the phase number(s) in the box
                     T=nothing ) 
 
 See the documentation of the GMG routine
 
 """
-function  AddEllipsoid!(model::Model;      # required input
-                    cen=Tuple{3}, axes=Tuple{3},                           # center and semi-axes of the ellpsoid
-                    Origin=nothing, StrikeAngle=0, DipAngle=0,             # origin & dip/strike
-                    phase = ConstantPhase(1),                              # Sets the phase number(s) in the box
-                    T=nothing )                                            # Sets the thermal structure (various fucntions are available)
-
-    AddEllipsoid!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, 
-                    cen=cen, axes=axes, Origin=Origin, StrikeAngle=StrikeAngle, DipAngle=DipAngle,
-                    phase=phase, T=T)
-
-    return nothing
-end
+AddEllipsoid!(model::Model; kwargs...) = AddEllipsoid!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 

src/LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl

@@ -83,10 +83,15 @@ include("Utils.jl")
 export  add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!,
         add_phasetransition!, add_dike!, add_geom!        
 
-
 include("ErrorChecking.jl")
 export Check_LaMEM_Model
 
+# optional dependency
+#if !isdefined(Base, :get_extension)
+#    include("MakieExt.jl")
+#end
+
+
 end
 
 

src/LaMEM_ModelGeneration/Model.jl

@@ -1,10 +1,11 @@
 # This is the main LaMEM Model struct
 using GeophysicalModelGenerator.GeoParams
-import LaMEM.Run: run_lamem
+import LaMEM.Run: run_lamem, run_lamem_save_grid
+using LaMEM.Run.LaMEM_jll
 
-export Model, Write_LaMEM_InputFiles
+export Model, Write_LaMEM_InputFiles, create_initialsetup
 
-"""
+""";
     Model
 
 Structure that holds all the information to create a LaMEM input file
@@ -176,15 +177,49 @@ Performs a LaMEM run for the parameters that are specified in `model`
 """
 function run_lamem(model::Model, cores::Int64=1, args::String=""; wait=true)
 
+    create_initialsetup(model, cores, args);    
+
+    run_lamem(model.Output.param_file_name, cores, args; wait=wait)
+
+    return nothing
+end
+
+"""
+    create_initialsetup(model::Model, cores::Int64=1, args::String="")
+
+Creates the initial model setup of LaMEM from `model`, which includes:
+- Writing the LaMEM (*.dat) input file
+- Write the VTK file (if requested when `model.Output.write_VTK_setup=true`)
+- Write the marker files to disk (if `model.ModelSetup.msetup="files"`)
+
+"""
+function create_initialsetup(model::Model, cores::Int64=1, args::String="")
+
+    # Move to the working directory
+    cur_dir = pwd()
+    if !isempty(model.Output.out_dir)
+        if !exist(model.Output.out_dir);  mkdir(model.Output.out_dir); end # create directory if needed
+        cd(model.Output.out_dir)
+    end
+
     Write_LaMEM_InputFile(model, model.Output.param_file_name)
 
     if model.ModelSetup.msetup=="files"
         # write marker files to disk before running LaMEM
+        Model3D = CartData(model.Grid.Grid, (Phases=model.Grid.Phases,Temp=model.Grid.Temp));
 
-    end
+        if cores>1
+           # PartFile = CreatePartitioningFile(model.Output.param_file_name, cores, LaMEM_options=args)
+
+            PartFile = run_lamem_save_grid(model.Output.param_file_name, cores)
 
-    run_lamem(model.Output.param_file_name, cores, args; wait=wait)
+            Save_LaMEMMarkersParallel(Model3D, PartitioningFile=PartFile)
+        else
+            Save_LaMEMMarkersParallel(Model3D)
+        end
+    end
 
+    cd(cur_dir)
     return nothing
 end
 

src/LaMEM_ModelGeneration/ModelSetup.jl

@@ -9,7 +9,7 @@ export ModelSetup, Write_LaMEM_InputFile, geom_Sphere, set_geom!
 """
 Base.@kwdef mutable struct ModelSetup
     "Setup type - can be `geom` (phases are assigned from geometric primitives), `files` (from julia input), `polygons` (from geomIO input, which requires `poly_file` to be specified) "
-    msetup::String          = "geom"              
+    msetup::String          = "files"              
 
     "add random noise to the particle location"
     rand_noise::Int64      = 1                 

src/LaMEM_ModelGeneration/Utils.jl

@@ -1,5 +1,6 @@
 # Contains a number of useful functions
-export  add_phase!, rm_phase!, rm_last_phase!, add_softening!, add_phasetransition!, 
+export  add_phase!, rm_phase!, rm_last_phase!, replace_phase!,
+        add_softening!, add_phasetransition!, 
         add_dike!, add_geom! 
 
 
@@ -12,6 +13,17 @@ function add_phase!(model::Model, phase::Phase)
     return nothing
 end
 
+"""
+    add_phase!(model::Model, phases...) 
+Add several phases @ once.
+"""
+function add_phase!(model::Model, phases...) 
+    for phase in phases
+        push!(model.Materials.Phases, phase);
+    end
+end
+
+
 """
     rm_last_phase!(model::Model, phase::Phase)
 This removes the last added `phase` from `model`
@@ -35,6 +47,26 @@ function rm_phase!(model::Model, ID::Int64)
 end
 
 
+"""
+    replace_phase!(model::Model, phase_new::Phase; ID::Int64=nothing, Name::String=nothing)
+
+This replaces a `phase` within a LaMEM Model Setup `model` with `phase_new` either based on its `Name` or `ID`. 
+Note that it is expected that only one such phase is present in the current setup.
+"""
+function replace_phase!(model::Model, phase_new::Phase; ID::Union{Nothing,Int64}=nothing, Name::Union{Nothing,String}=nothing) 
+    id_vec   = [phase.ID for phase in model.Materials.Phases]
+    name_vec = [phase.Name for phase in model.Materials.Phases]
+    if !isnothing(ID)
+        id = findfirst(id_vec .== ID)
+    elseif !isnothing(Name)
+        id = findfirst(name_vec .== Name)
+    end
+    model.Materials.Phases[id] = phase_new
+
+    return nothing
+end
+
+
 """
     add_petsc!(model::Model, option::String) 
 

src/MakieExt.jl

@@ -0,0 +1,65 @@
+# Plotting extensions, that are only loaded when GLMakie is available
+println("adding GLMakie plotting extensions for LaMEM")
+
+using LaMEM, GeophysicalModelGenerator
+using .GLMakie
+
+export plot_initial_setup
+
+"""
+    plot_initial_setup(model::LaMEM.Model; x=nothing, y=nothing, z=nothing, phases=true)
+
+Plots a 2D setup through the LaMEM model setup `model`. if `phases=true` we plot phases; otherwise temperature
+"""
+function plot_initial_setup(model::LaMEM.Model; x=nothing, y=nothing, z=nothing, phases=true)
+
+    Model3D = CartData(model.Grid.Grid, (Phases=model.Grid.Phases,Temp=model.Grid.Temp));
+
+    if !isnothing(z); Depth_level = z*km; else Depth_level=nothing; end
+    if !isnothing(x); Lon_level = x; else Lon_level=nothing; end
+    if !isnothing(y); Lat_level = y; else Lat_level=nothing; end
+
+    Cross = CrossSectionVolume(Model3D, Depth_level=Depth_level, Lat_level=Lat_level, Lon_level=Lon_level)
+
+    # Flatten arrays + coordinates
+    if      !isnothing(x)
+        Phase2D, Temp2D =   dropdims(Cross.fields.Phases,dims=1), dropdims(Cross.fields.Temp,dims=1)
+        X2D, Z2D = dropdims(Cross.y.val,dims=1), dropdims(Cross.z.val,dims=1)
+        x_vec, z_vec = X2D[:,1],  Z2D[1,:]
+        title_str = "x = $x"
+        x_str,z_str = "y","z"
+    elseif !isnothing(y)
+        Phase2D, Temp2D =   dropdims(Cross.fields.Phases,dims=2), dropdims(Cross.fields.Temp,dims=2)
+        X2D, Z2D = dropdims(Cross.x.val,dims=2), dropdims(Cross.z.val,dims=2)
+        x_vec, z_vec = X2D[:,1],  Z2D[1,:]
+        title_str = "y = $y"
+        x_str,z_str = "x","z"
+    elseif !isnothing(z)
+        Phase2D, Temp2D =   dropdims(Cross.fields.Phases,dims=3), dropdims(Cross.fields.Temp,dims=3)
+        X2D, Z2D = dropdims(Cross.x.val,dims=3), dropdims(Cross.y.val,dims=3)
+
+        x_vec, z_vec = X2D[:,1],  Z2D[1,:]
+        title_str = "z = $z"
+        x_str,z_str = "x","y"
+    end
+
+
+    # Create figure
+    if phases
+        data = Phase2D
+        title_str = "Phases "*title_str
+    else
+        data = Temp2D
+        title_str = "Temperature "*title_str
+    end
+
+    fig = Figure()
+    ax = Axis(fig[1, 1], aspect = DataAspect(), title = title_str, xlabel=x_str, ylabel=z_str)
+    hm = heatmap!(ax, x_vec, z_vec, data)
+
+    Colorbar(fig[:, end+1], hm)
+
+    display(fig)
+
+    return nothing
+end