diff --git a/Project.toml b/Project.toml index bd2689be..12a7a346 100644 --- a/Project.toml +++ b/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" diff --git a/src/LaMEM.jl b/src/LaMEM.jl index ae83114d..08304f93 100644 --- a/src/LaMEM.jl +++ b/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 diff --git a/src/LaMEM_ModelGeneration/ErrorChecking.jl b/src/LaMEM_ModelGeneration/ErrorChecking.jl index 968dd7ad..57e0d3f8 100644 --- a/src/LaMEM_ModelGeneration/ErrorChecking.jl +++ b/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") diff --git a/src/LaMEM_ModelGeneration/Materials.jl b/src/LaMEM_ModelGeneration/Materials.jl index 2eb0d994..e1edea24 100644 --- a/src/LaMEM_ModelGeneration/Materials.jl +++ b/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") 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] diff --git a/src/LaMEM_ModelGeneration/Utils.jl b/src/LaMEM_ModelGeneration/Utils.jl index 46a35323..ac072fa2 100644 --- a/src/LaMEM_ModelGeneration/Utils.jl +++ b/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 \ No newline at end of file +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 diff --git a/test/runtests.jl b/test/runtests.jl index 9a521b34..74fb3df9 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_GeoParams_integration.jl b/test/test_GeoParams_integration.jl new file mode 100644 index 00000000..07b3a77f --- /dev/null +++ b/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 \ No newline at end of file