Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use GeoParams rheologies in LaMEM.jl #40

Merged
merged 6 commits into from
Feb 19, 2024
Merged

Use GeoParams rheologies in LaMEM.jl #40

merged 6 commits into from
Feb 19, 2024

Conversation

boriskaus
Copy link
Member

@boriskaus boriskaus commented Feb 19, 2024

The GeoParams.jl package has many additional laboratory-derived creep laws implemented compared to LaMEM (much of which is thanks to @Dr3ven).
With this PR it becomes possible to directly use any of the GeoParams diffusion/dislocation creep laws in LaMEM:

julia> using GeoParams, LaMEM
julia> model  = Model( Grid(nel=(8, 8), x=[-1,1], z=[-1,1]),
                    Time(nstep_max=1, dt=1e-6, dt_max=1, dt_min=1e-10, time_end=100), 
                    BoundaryConditions(exx_strain_rates=[1e-15] ),   
                    Output(out_dir="0D_1"));
julia> g = SetDislocationCreep(GeoParams.Dislocation.serpentinite_Raleigh_1965);
julia> rheology = Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000)
Phase 0 (rheology): 
  rho       = 3000.0 
  GeoParams = [Tumut Pond Serpentinite | Raleigh and Paterson (1965); ] 
julia> add_phase!(model, rheology)    
julia> run_lamem(model)

A few remarks:

  1. the add_phase! function internally sets the corresponding values by calling add_geoparams_rheologies:
julia> rheology1 = add_geoparams_rheologies(Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000))
Phase 0 (rheology): 
  rho       = 3000.0 
  Bn        = 4.025695644673317e-23 
  En        = 66000.0 
  Vn        = 0.0 
  n         = 2.8 
  GeoParams = [Tumut Pond Serpentinite | Raleigh and Paterson (1965); ] 
  1. If you specify diffusion creep, you also have to indicate the grain size (in meters) as in:
julia> rheology1 = Phase(ID=0,Name="rheology",GeoParams=[g], rho=3000, grainsize=100e-6)
Phase 0 (rheology): 
  rho       = 3000.0 
  GeoParams = [Dry Anorthite | Rybacki et al. (2006); ] 
  grainsize = 0.0001
  1. You can specify both diffusion creep and dislocation creep rheologies in the GeoParams field:
julia> diff = SetDiffusionCreep(GeoParams.Diffusion.dry_olivine_Hirth_2003);
julia> disl = SetDislocationCreep(GeoParams.Dislocation.dry_olivine_Hirth_2003);
julia> rheology = Phase(ID=0,Name="rheology",GeoParams=[diff, disl], rho=3000, grainsize=100e-6)
julia> rm_phase!(model)
julia> add_phase!(model, rheology)
julia> model.Materials.Phases
1-element Vector{Phase}:
 Phase 0 (rheology): 
  rho       = 3000.0 
  Bd        = 0.0023773397886916654 
  Ed        = 375000.0 
  Vd        = 1.0e-5 
  Bn        = 6.514566364114834e-16 
  En        = 530000.0 
  Vn        = 1.4e-5 
  n         = 3.5 
  GeoParams = [Dry Olivine | Hirth & Kohlstedt (2003); Dry Olivine | Hirth & Kohlstedt (2003); ] 
  grainsize = 0.0001 

If you specify >1 of each type, the last one will be taken.

Internally, we simply translate the GeoParams values to the corresponding Vn,En,n,Bn and Vd,Ed,Bd values that LaMEM expected.

We also added a function stress_strainrate_0D that can be used to perform 0D rheology experiments and create stress-strain rate plots by running LaMEM at low resolution for a range of strain rates.

@boriskaus boriskaus changed the title Use GeoParams rheology in LaMEM.jl Use GeoParams rheologies in LaMEM.jl Feb 19, 2024
@boriskaus boriskaus merged commit f0ce329 into main Feb 19, 2024
20 checks passed
@boriskaus boriskaus deleted the bk-geoparams branch February 19, 2024 07:39
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant