add test for GP rheology

src/LaMEM.jl

@@ -37,7 +37,8 @@ export  LaMEM_Model, Model, Write_LaMEM_InputFile, create_initialsetup,
         add_phasetransition!, add_dike!, add_geom!, set_air, copy_phase,
         add_topography!, AboveSurface!, BelowSurface!,
         prepare_lamem, isdefault, hasplasticity,
-        add_geoparams_rheologies    
+        add_geoparams_rheologies,
+        stress_strainrate_0D    
 
 
 using .Run.LaMEM_jll

src/LaMEM_ModelGeneration/Materials.jl

@@ -289,6 +289,13 @@ Base.@kwdef mutable struct Phase
     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
 
 
@@ -297,7 +304,14 @@ function add_geoparams_rheologies(phase::Phase)
         # NOTE: this needs checking; likely that B in LaMEM is defined differently!
         for ph in phase.GeoParams
             if isa(ph, DiffusionCreep)
-                phase.Bd = NumValue(ph.A)*ph.FT/ph.FE*1e-2^(NumValue(ph.p))
+                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)
@@ -847,7 +861,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)) & (p != :GeoParams)
+            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
 
 
 
@@ -405,4 +406,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 = add_geoparams_rheologies(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