Better default parameters

src/LaMEM.jl

@@ -36,7 +36,7 @@ export  LaMEM_Model, Model, Write_LaMEM_InputFile, create_initialsetup,
         add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!, add_phaseaggregate!,
         add_phasetransition!, add_dike!, add_geom!, set_air, copy_phase,
         add_topography!, AboveSurface!, BelowSurface!,
-        prepare_lamem    
+        prepare_lamem, isdefault, hasplasticity    
 
 
 using .Run.LaMEM_jll

src/LaMEM_ModelGeneration/DefaultParams.jl

@@ -35,6 +35,57 @@ function UpdateDefaultParameters(model::Model)
         model.SolutionParams.act_p_shift = 1
     end
 
+    # Output some fields @ all times
+    model.Output.out_density          = 1
+    model.Output.out_j2_dev_stress    = 1 
+    model.Output.out_j2_strain_rate   = 1 
+    model.Output.out_pressure         = 1 
+    model.Output.out_temperature      = 1  
+
+    if isdefault(model.Solver,Solver())
+        # the LaMEM default (superlu_dist) currently doesn't work with PETSc_jll
+        model.Solver.DirectSolver = "mumps";    
+    end
+
+    # if using multigrid in a 2D setup
+    if model.Solver.SolverType=="multigrid" &&  model.Grid.nel_y[1]==1
+        push!(model.Solver.PETSc_options,"-da_refine_y 1") 
+    end
+
+    # if we have a free surface, you'll generally want output  
+    if  model.FreeSurface.surf_use==1
+        model.Output.out_surf=1
+        model.Output.out_surf_pvd         = 1 
+        model.Output.out_surf_velocity    = 1 
+        model.Output.out_surf_topography  = 1 
+        model.Output.out_surf_amplitude   = 1  
+    end
+
+    # Scaling: if we use default values, employ smarter default values based on the model setup
+    if isdefault(model.Scaling,Scaling())
+        le = abs(diff(model.Grid.coord_z)[1])*km    # length
+        η  = model.SolutionParams.eta_ref*Pas       # viscosity
+        model.Scaling = Scaling(GEO_units(length=le, viscosity = η))
+    end
+
+    # exx_strain_rates: no need to specify exx_num_periods
+    if length(model.BoundaryConditions.exx_strain_rates)==1 && model.BoundaryConditions.exx_num_periods>1
+        # we only specified a single strainrate so it'll be constant with time
+        model.BoundaryConditions.exx_num_periods = 1
+        model.BoundaryConditions.exx_time_delims = [1e20];
+    end
+
+    if  model.FreeSurface.surf_use==1
+        # if we use a free surface & surf_level==nothing, set it to zero
+        if isnothing(model.FreeSurface.surf_level)
+            model.FreeSurface.surf_level = 0.0;
+        end
+        if isnothing(model.FreeSurface.surf_air_phase)
+            model.FreeSurface.surf_air_phase = 0
+        end
+
+
+    end
 
     return model
 end

src/LaMEM_ModelGeneration/FreeSurface.jl

@@ -129,7 +129,8 @@ function Write_LaMEM_InputFile(io, d::FreeSurface)
 
     if surf_use==1
         for f in fields
-            if getfield(d,f) != getfield(Reference,f) && (f != :Topography)
+            if (getfield(d,f) != getfield(Reference,f) && (f != :Topography)) 
+
                 # only print if value differs from reference value
                 name = rpad(String(f),15)
                 comment = get_doc(FreeSurface, f)

src/LaMEM_ModelGeneration/GMG_interface.jl

@@ -2,7 +2,7 @@
 #
 # Some wrappers around GMG routines
 
-import GeophysicalModelGenerator: AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AboveSurface, BelowSurface
+import GeophysicalModelGenerator: AddBox!, AddLayer!, AddSphere!, AddEllipsoid!, AddCylinder!, AboveSurface, BelowSurface
 export AboveSurface!, BelowSurface!
 
 """
@@ -12,11 +12,22 @@ export AboveSurface!, BelowSurface!
             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
+See the documentation of the GMG routine for the full options.
 
 """
 AddBox!(model::Model; kwargs...) = AddBox!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 
 
+"""
+    AddLayer!(model::Model; xlim, ylim, zlim=Tuple{2},
+            phase = ConstantPhase(1),
+            T=nothing )
+
+Adds a layer with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models
+See the documentation of the GMG routine for the full options.
+
+"""
+AddLayer!(model::Model; kwargs...) = AddLayer!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) 
+
 
 """
     AddSphere!(model::Model; cen=Tuple{3}, radius=Tuple{1}, phase = ConstantPhase(1), T=nothing)

src/LaMEM_ModelGeneration/Materials.jl

@@ -120,6 +120,9 @@ Base.@kwdef mutable struct Phase
     "stabilization viscosity (default is eta_min)"
     eta_st::Union{Nothing,Float64}       = nothing   
 
+    "viscoplastic plasticity regularisation viscosity"
+    eta_vp::Union{Nothing,Float64}       = nothing   
+
     "pore-pressure ratio"
     rp::Union{Nothing,Float64}       = nothing   
 
@@ -221,7 +224,7 @@ Base.@kwdef mutable struct Softening
     A::Float64           =   0.7     
 
     "Material length scale (in selected units, e.g. km in geo)"
-    Lm::Float64          =   0.2    
+    Lm::Union{Nothing,Float64}          =   nothing   
 
     # healing parameters
     "APS when healTau2 activates"

src/LaMEM_ModelGeneration/Model.jl

@@ -160,7 +160,13 @@ function Write_LaMEM_InputFile(d::Model, fname::String="input.dat"; dir=pwd())
         # If we want to write an input file 
         Write_Paraview(CartData(d.Grid.Grid, (Phases=d.Grid.Phases,Temp=d.Grid.Temp)),"Model3D")
     end
-
+
+    if any(hasplasticity.(d.Materials.Phases))
+        # We have plasticity, so we likely want to see that
+        d.Output.out_plast_strain    = 1     # accumulated plastic strain
+        d.Output.out_plast_dissip    = 1      # plastic dissipation
+    end
+
     io = open(fname,"w")
 
     Write_LaMEM_InputFile(io, d.Scaling)
@@ -188,13 +194,15 @@ Performs a LaMEM run for the parameters that are specified in `model`
 """
 function run_lamem(model::Model, cores::Int64=1, args::String=""; wait=true)
 
+    cur_dir = pwd(); 
+
+    #if !isdir(model.Output.out_dir); mkdir(model.Output.out_dir); end # create directory if needed
     create_initialsetup(model, cores, args);    
 
-    cur_dir = pwd(); 
     if !isempty(model.Output.out_dir)
         cd(model.Output.out_dir)
     end
-    
+
     run_lamem(model.Output.param_file_name, cores, args; wait=wait)
 
     cd(cur_dir)
@@ -219,9 +227,7 @@ function prepare_lamem(model::Model, cores::Int64=1, args::String=""; verbose=fa
 
     println("Creating LaMEM input files in the directory: $(model.Output.out_dir)")
     cur_dir = pwd(); 
-    if !isempty(model.Output.out_dir)
-        cd(model.Output.out_dir)
-    end
+
     create_initialsetup(model, cores, args,  verbose=verbose);    
 
     cd(cur_dir)

src/LaMEM_ModelGeneration/ModelSetup.jl

@@ -54,7 +54,7 @@ Base.@kwdef mutable struct ModelSetup
     nmark_avd::Vector{Int64}        = [3, 3, 3]             
 
     "max number of same phase markers per subcell (subgrid marker control)"
-    nmark_sub::Int64                = 1        
+    nmark_sub::Int64                = 3        
 
     "Different geometric primitives that can be selected if we `msetup``=`geom`; see `geom_Sphere`"
     geom_primitives::Vector    = []
@@ -169,12 +169,11 @@ function Write_LaMEM_InputFile(io, d::ModelSetup)
             (f == :msetup)      ||
             (f == :rand_noise)  ||
             (f == :nmark_lim)   ||
+            (f == :nmark_sub)   ||
             (f == :advect)      ||
             (f == :interp)      ||
             (f == :mark_ctrl)
 
-
-
 
             if (f != :geom_primitives)
                 # only print if value differs from reference value

src/LaMEM_ModelGeneration/SolutionParams.jl

@@ -183,7 +183,9 @@ function Write_LaMEM_InputFile(io, d::SolutionParams)
         if getfield(d,f) != getfield(Reference,f) ||
             (f == :eta_ref) ||
             (f == :gravity) ||
-            (f == :act_p_shift)
+            (f == :act_p_shift) ||
+            (f == :init_guess) ||
+            (f == :p_lim_plast)
 
             # only print if value differs from reference value
             name = rpad(String(f),15)

src/LaMEM_ModelGeneration/Utils.jl

@@ -5,7 +5,9 @@ export  add_phase!, rm_phase!, rm_last_phase!, replace_phase!,
         add_vbox!, rm_last_vbox!, rm_vbox!,
         add_softening!, add_phasetransition!, add_phaseaggregate!,
         add_dike!, add_geom! , cross_section,
-        set_air, copy_phase
+        set_air, copy_phase,
+        isdefault, hasplasticity
+
 
 
 """
@@ -370,4 +372,35 @@ function add_topography!(model::Model, topography::CartData; surf_air_phase=0, s
     model.BoundaryConditions.open_top_bound=open_top_bound
 
     return nothing
+end
+
+
+"""
+    isdefault(s1::S, s_default::S) 
+
+Checks whether a struct `s1` has default parameters `s_default`
+"""
+function isdefault(s1, s_default)
+
+    default = true
+    for f in fieldnames(typeof(s1))
+        if getfield(s1,f) !=   getfield(s_default,f)
+            default = false
+        end
+    end
+
+    return default
+end
+
+"""
+    hasplasticity(p::Phase)
+
+`true` if `p` contains plastic parameters (cohesion or friction angle)
+"""
+function hasplasticity(p::Phase)
+    plastic = false
+    if !isnothing(p.ch) || !isnothing(p.fr)
+        plastic = true
+    end
+    return plastic
 end

src/PlotsExt.jl

@@ -4,16 +4,21 @@ println("Adding Plots.jl plotting extensions for LaMEM")
 using LaMEM, GeophysicalModelGenerator, Statistics, DelimitedFiles
 using .Plots
 import .Plots: heatmap
-export plot_topo, plot_cross_section, plot_phasediagram
+export plot_topo, plot_cross_section, plot_phasediagram, plot_cross_section_simulation
 
 """
-    plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol=:phase,  dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
+    plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol=:phase,  
+                        timestep::Union{Nothing, Int64}=nothing,
+                        dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
 
 This plots a cross-section through the LaMEM `model`, of the field  `field`. If that is a vector or tensor field, specify `dim` to indicate which of the fields you want to see.
 If the keyword `timestep` is specified, it will load a timestep 
 """
-function plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol=:phase, timestep=nothing, dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
-
+function plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol=:phase, 
+        timestep::Union{Nothing, Int64}=nothing, dim=1, 
+        x=nothing, y=nothing, z=nothing, 
+        aspect_ratio::Union{Real, Symbol}=:equal)
+
     if !isnothing(timestep)
         # load a particular timestep
         data_cart, time = Read_LaMEM_timestep(model,timestep)
@@ -51,6 +56,24 @@ function plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol
 end
 
 
+"""
+    plot_cross_section_simulation(model::Union{Model,CartData}, args...; field::Symbol=:phase,  
+                        dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
+
+As `plot_cross_section`, but for the entire simulation instead of a single timestep.
+"""
+function plot_cross_section_simulation(model::Union{Model,CartData}, args...; field::Symbol=:phase, 
+        dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
+
+    Timesteps,_,_ = Read_LaMEM_simulation(model);
+    for timestep_val in Timesteps
+        plot_cross_section(model, args, field=field, timestep=timestep_val, dim=dim, x=x, y=y, z=z, aspect_ratio=aspect_ratio)
+    end
+
+    return nothing
+end
+
+
 """
     plot_topo(topo::CartData, args...)
 

test/mesh_refinement_test.jl

@@ -31,7 +31,7 @@ using GeophysicalModelGenerator
     # read last timestep
     data,time = Read_LaMEM_timestep(model,last=true);
 
-    @test  sum(data.fields.velocity[3][:,:,:]) ≈ 0.36892682f0 # check Vz
+    @test  sum(data.fields.velocity[3][:,:,:]) ≈ 0.3680135f0 # check Vz
 
     # cleanup the directory
     rm(model.Output.out_dir, force=true, recursive=true)

test/test_julia_setups.jl

@@ -151,7 +151,7 @@ using GeophysicalModelGenerator
     # read last timestep
     data,time = Read_LaMEM_timestep(model,last=true);
 
-    @test  sum(data.fields.phase) ≈ 29047.736f0
+    @test  sum(data.fields.phase) ≈ 29060.664f0
 
     # cleanup the directory
     rm(model.Output.out_dir, force=true, recursive=true)