Merge pull request #40 from JuliaGeodynamics/bk-geoparams

Use GeoParams rheologies in LaMEM.jl

Project.toml

@@ -6,8 +6,8 @@ version = "0.2.11"
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
 GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742"
+GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 LaMEM_jll = "15d6fa20-f789-5486-b71b-22b4ac8eb1c1"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"

src/LaMEM.jl

@@ -36,7 +36,9 @@ 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, isdefault, hasplasticity    
+        prepare_lamem, isdefault, hasplasticity,
+        add_geoparams_rheologies,
+        stress_strainrate_0D    
 
 
 using .Run.LaMEM_jll

src/LaMEM_ModelGeneration/ErrorChecking.jl

@@ -19,8 +19,7 @@ function Check_LaMEM_Model(m::Model)
     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, ))` ")
+        @warn "Your initial `Temp` grid is constant, as is your initial `Phases` grid. \n Is that intended? \n In most cases, you would want to set some variability in the initial conditions, \n for example with the `GeophysicalModelGenerator` function `AddSphere!(model,cen=(0.0,0.0,0.0), radius=(0.15, ))` "
     end
 
     if (m.Solver.SolverType!="direct") &&  (m.Solver.SolverType!="multigrid")

src/LaMEM_ModelGeneration/Materials.jl

@@ -1,6 +1,7 @@
 # Specify Material properties
+using GeoParams
 export Materials, Phase, Softening, PhaseAggregate, PhaseTransition, Dike, Write_LaMEM_InputFile
-
+export add_geoparams_rheologies
 
 
 """
@@ -272,6 +273,56 @@ Base.@kwdef mutable struct Phase
 
     "melt fraction viscosity correction factor (positive scalar)"
     mfc::Union{Nothing,Float64}       = nothing 
+
+    """
+    GeoParams creeplaws 
+    
+    Set diffusion or dislocation creeplaws as provided by the GeoParams package:
+    
+    ```julia
+    julia> using GeoParams
+    julia> a = SetDiffusionCreep(GeoParams.Diffusion.dry_anorthite_Rybacki_2006);
+    julia> p = Phase(ID=1,Name="test", GeoParams=[a]);
+    ```
+    Note that GeoParams should be a vector, as you could, for example, have diffusion and dislocation creep parameters
+    
+    Note also that this will overwrite any other creeplaws provided in the Phase struct.
+    """
+    GeoParams::Union{Nothing,Vector{AbstractCreepLaw}}       = nothing 
+
+    """
+    grainsize [m] (not used in LaMEM!)
+    This is not actually used in LaMEM, but is required when setting diffusion creep parameters by using GeoParams 
+    """
+    grainsize::Union{Nothing, Float64} = nothing
+
+end
+
+
+function add_geoparams_rheologies(phase::Phase)
+    if !isnothing(phase.GeoParams)
+        # NOTE: this needs checking; likely that B in LaMEM is defined differently!
+        for ph in phase.GeoParams
+            if isa(ph, DiffusionCreep)
+                d0       = phase.grainsize
+                if isnothing(d0)
+                    error("If you use GeoParams to set diffusion creep, you need to specify the grainsize in the Phase info")
+                end
+
+                #phase.Bd = NumValue(ph.A)*ph.FT/ph.FE*d0^(NumValue(ph.p))
+                phase.Bd = NumValue(ph.A)*d0^(NumValue(ph.p))*ph.FT/ph.FE
+
+                phase.Ed = ph.E
+                phase.Vd = ph.V
+            elseif isa(ph, DislocationCreep)
+                phase.Bn = NumValue(ph.A)*ph.FT^NumValue(ph.n)/ph.FE
+                phase.En = ph.E
+                phase.Vn = ph.V
+                phase.n  = ph.n
+            end
+        end
+    end
+    return phase
 end
 
 function show(io::IO, d::Phase)
@@ -280,9 +331,22 @@ function show(io::IO, d::Phase)
 
     # print fields
     for f in fields
-        if !isnothing(getfield(d,f)) & (f != :ID) & (f != :Name)
+        if !isnothing(getfield(d,f)) & (f != :ID) & (f != :Name) & (f != :GeoParams)
             printstyled(io,"  $(rpad(String(f),9)) = $(getfield(d,f)) \n")        
         end
+
+        # if we have GeoParams creep data, print it differently
+        if f == :GeoParams && !isnothing(getfield(d,f))
+            g = getfield(d,f);
+            names = "["
+            for i=1:length(g)
+                name = GeoParams.uint2str(g[i].Name)
+                names = names*"$(name); "
+            end
+            names = names*"]"
+            printstyled(io,"  $(rpad(String(f),9)) = $names \n")        
+        end
+
     end
 
     return nothing
@@ -292,18 +356,27 @@ function show_short(d::Phase)
     fields    = fieldnames(typeof(d))
     str = "Phase("
     for (i,f) in enumerate(fields)
-        if !isnothing(getfield(d,f))
+        if !isnothing(getfield(d,f)) & (f != :GeoParams) 
             str=str*"$(String(f))=$(getfield(d,f))"        
             if i<length(fields)
                 str=str*","
             end
+        elseif !isnothing(getfield(d,f)) && f == :GeoParams 
+            g = getfield(d,f);
+            names = "["
+            for i=1:length(g)
+                name = GeoParams.uint2str(g[i].Name)
+                names = names*"$(name);"
+            end
+            names = names*"]"
+            str *= "  $(rpad(String(f),9)) = $names"
         end
+
     end
     str=str*")"
     return str
 end
 
-
 """
     Defines strain softening parameters
 
@@ -798,7 +871,7 @@ function Write_LaMEM_InputFile(io, d::Materials)
         println(io, "   <MaterialStart>")
         phase_fields    = fieldnames(typeof(phase))
         for p in phase_fields
-            if !isnothing(getfield(phase,p))
+            if !isnothing(getfield(phase,p)) & (p != :GeoParams) & (p != :grainsize)
                 name = rpad(String(p),15)
                 comment = get_doc(Phase, p)
                 comment = split(comment,"\n")[1]

src/LaMEM_ModelGeneration/Utils.jl

@@ -6,7 +6,8 @@ export  add_phase!, rm_phase!, rm_last_phase!, replace_phase!,
         add_softening!, add_phasetransition!, add_phaseaggregate!,
         add_dike!, add_geom! , cross_section,
         set_air, copy_phase,
-        isdefault, hasplasticity
+        isdefault, hasplasticity,
+        stress_strainrate_0D
 
 
 
@@ -55,7 +56,8 @@ end
 This adds a `phase` (with material properties) to `model`
 """
 function add_phase!(model::Model, phase::Phase) 
-    push!(model.Materials.Phases, phase);
+    phase_added = add_geoparams_rheologies(phase)
+    push!(model.Materials.Phases, phase_added);
     return nothing
 end
 
@@ -65,11 +67,13 @@ Add several phases @ once.
 """
 function add_phase!(model::Model, phases...) 
     for phase in phases
-        push!(model.Materials.Phases, phase);
+        phase_added = add_geoparams_rheologies(phase)
+        push!(model.Materials.Phases, phase_added);
     end
 end
 
 
+
 """
     rm_last_phase!(model::Model, phase::Phase)
 This removes the last added `phase` from `model`
@@ -403,4 +407,46 @@ function hasplasticity(p::Phase)
         plastic = true
     end
     return plastic
-end
+end
+
+
+
+"""
+    τ = stress_strainrate_0D(rheology, ε_vec::Vector; n=8, T=700, nstep_max=2, clean=true)
+
+Computes the stress for a given strain rate and 0D rheology setup, for viscous creep rheologies. 
+    `n` is the resolution in `x,z`, `T` the temperature, `nstep_max` the number of time steps, `ε_vec`
+    the strainrate vector (in 1/s). 
+
+"""
+function stress_strainrate_0D(rheology, ε_vec::Vector; n=8, T=700, nstep_max=2, clean=true)
+
+    # Main model setup
+    model  = Model( Grid(nel=(n, n), x=[-1,1], z=[-1,1]),
+                    Time(nstep_max=nstep_max, dt=1e-6, dt_max=1, dt_min=1e-10, time_end=100), 
+                    BoundaryConditions(exx_strain_rates=[1e-15] ),   
+                    Output(out_dir="0D_1"))
+
+    rm_phase!(model)
+    add_phase!(model, rheology)
+    model.Grid.Temp.=T;
+
+    τ = zero(ε_vec)
+    for (i,ε) in enumerate(ε_vec)
+        # run the simulation on 1 core
+        model.Output.out_dir="0D_$i"
+        model.BoundaryConditions.exx_strain_rates = [ε]
+        run_lamem(model, 1); #run
+        data,_ = Read_LaMEM_timestep(model, last=true) # read
+
+        @show extrema(data.fields.j2_dev_stress)
+
+        τ[i] = Float64.(sum(extrema(data.fields.j2_dev_stress))/2) # in MPa
+
+        if clean
+            rm("0D_$i",recursive=true)
+        end
+    end
+
+    return τ
+end

test/runtests.jl

@@ -8,6 +8,7 @@ include("run_lamem_save_grid_test.jl")
 include("mesh_refinement_test.jl")
 include("read_logfile.jl")
 include("test_compression.jl")
+include("test_GeoParams_integration.jl")
 
 if !Sys.iswindows()
     # clean up

test/test_GeoParams_integration.jl

@@ -0,0 +1,70 @@
+using Test
+using GeoParams, LaMEM
+
+
+@testset "GeoParams 0D rheology" begin
+
+    T0 = 1100;
+    ε_vec  = 10.0.^Vector(-17.0:1.0:-13.0);
+
+
+    # compute stress for given strainrate 
+
+    # 1) Linear viscous rheology
+    rheology = Phase(ID=0,Name="matrix",eta=1e20,rho=3000)
+    @show rheology
+    τ_linear = stress_strainrate_0D(rheology, ε_vec, T=T0, nstep_max=2)
+    τ_anal = 2.0 .* rheology.eta .* ε_vec ./ 1e6;
+    @show τ_linear τ_anal
+    @test sum(τ_linear-τ_anal) ≈ 0.0 atol=1e-6 # check stress
+
+    # 2) Dislocation creep rheology
+    # Serpentinite --- 
+    rheology = Phase(ID=0,Name="matrix",disl_prof="Tumut_Pond_Serpentinite-Raleigh_Paterson_1965", rho=3000)
+    τ_num1 = stress_strainrate_0D(rheology, ε_vec; T=T0)
+
+    g = SetDislocationCreep(GeoParams.Dislocation.serpentinite_Raleigh_1965)
+    rheology1 = add_geoparams_rheologies(Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000))
+    τ_num2 = stress_strainrate_0D(rheology1, ε_vec; T=T0)
+    @show τ_num1 τ_num2
+    @test sum(τ_num1-τ_num2) ≈ 0.0 atol=1e-6 # check stress
+
+    # Olivine ---
+    rheology = Phase(ID=0,Name="matrix",disl_prof="Dry_Olivine-Ranalli_1995", Vn=0.0, rho=3000)
+    τ_num1 = stress_strainrate_0D(rheology, ε_vec; T=T0)
+
+    g = SetDislocationCreep(GeoParams.Dislocation.Dislocation.dry_olivine_Gerya_2019)
+    rheology1 = add_geoparams_rheologies(Phase(ID=0,Name="matrix",GeoParams=[g], rho=3000))
+    τ_num2 = stress_strainrate_0D(rheology1, ε_vec, T=T0)
+    @test sum(τ_num1-τ_num2) ≈ 0.0 atol=1e-6 # check stress
+
+    # Quartzite --- 
+    rheology = Phase(ID=0,Name="rheology",disl_prof="Wet_Quarzite-Ueda_et_al_2008", rho=3000)
+    τ_num1 = stress_strainrate_0D(rheology, ε_vec; T=T0)
+    #    (disl)   : Bn = 1.53539e-17 [1/Pa^n/s]  En = 154000. [J/mol]  n = 2.3 [ ]  
+    #    (disl)   : Bn = 1.53539e-17 [1/Pa^n/s]  En = 154000. [J/mol]  n = 2.3 [ ]  
+
+    g = SetDislocationCreep(GeoParams.Dislocation.wet_quartzite_Ueda_2008)
+    rheology1 = add_geoparams_rheologies(Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000))
+    τ_num2 = stress_strainrate_0D(rheology1, ε_vec; T=T0)
+    @test sum(τ_num1-τ_num2) ≈ 0.0 atol=1e-6 # check stress
+
+
+    # 3) Diffusion creep rheology
+    #
+    # Note: diffusion creep is grain size sensitive. Whereas the current version of LaMEM does not include grain size
+    # we take it into account through a prefactor (d0). In GeoParams, grainsize should be added as a parameter
+    # Yet, to convert one to the other, we need to specify grainsize in the Phase info 
+    rheology = Phase(ID=0,Name="rheology",diff_prof="Dry_Plagioclase_RybackiDresen_2000", rho=3000)
+    τ_num1 = stress_strainrate_0D(rheology, ε_vec; T=T0)
+
+    g = SetDiffusionCreep(GeoParams.Diffusion.dry_anorthite_Rybacki_2006)
+    rheology1 = Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000, grainsize=100e-6)
+    τ_num2 = stress_strainrate_0D(rheology1, ε_vec; T=T0)
+    @test sum(τ_num1-τ_num2) ≈ 0.0 atol=1e-6 # check stress
+
+    # Note: many of the other diffusion creep laws have a slightly different V exponent in GeoParams vs. LaMEM,
+    # and as of now we don't have a way yet to override this in GeoParams (to be added)
+    # Thats why we only have one test now
+
+end