Use GeoParams rheologies in LaMEM.jl

[boriskaus]

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); ] 

2) 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

3) 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.