first part of using GeoParams to set diffusion/dislocation creep (not…

… tested!)

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_ModelGeneration/Materials.jl

@@ -1,8 +1,8 @@
 # Specify Material properties
+using GeoParams
 export Materials, Phase, Softening, PhaseAggregate, PhaseTransition, Dike, Write_LaMEM_InputFile
 
 
-
 """
     Defines the material properties for each of the phases
 
@@ -272,6 +272,43 @@ 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 
+end
+
+
+function process_Phase(phase::Phase)
+    if !isnothing(phase.GeoParams)
+        # NOTE: this needs checking; likely that B in LaMEM is defined differently!
+        println("GeoParamsExt: adding creeplaw params")
+        for ph in phase.GeoParams
+            if isa(ph, DiffusionCreep)
+                phase.Bd = ph.A
+                phase.Ed = ph.E
+                phase.Vd = ph.V
+            elseif isa(ph, DislocationCreep)
+                phase.Bn = ph.A
+                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 +317,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
@@ -303,7 +353,6 @@ function show_short(d::Phase)
     return str
 end
 
-
 """
     Defines strain softening parameters
 

src/LaMEM_ModelGeneration/Utils.jl

@@ -65,11 +65,13 @@ Add several phases @ once.
 """
 function add_phase!(model::Model, phases...) 
     for phase in phases
-        push!(model.Materials.Phases, phase);
+
+        push!(model.Materials.Phases, process_Phase(phase));
     end
 end
 
 
+
 """
     rm_last_phase!(model::Model, phase::Phase)
 This removes the last added `phase` from `model`