Merge pull request #33 from JuliaGeodynamics/bk-phase-diagrams-melt

Test for phase diagrams

Project.toml

@@ -4,6 +4,7 @@ authors = ["Boris Kaus <kaus@uni-mainz.de>"]
 version = "0.2.7"
 
 [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"
@@ -15,6 +16,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DocStringExtensions = "0.9"
+DelimitedFiles = "1"
 GeoParams = "0.4, 0.5"
 GeophysicalModelGenerator = "0.4 - 0.6"
 Glob = "1"

src/LaMEM.jl

@@ -11,7 +11,7 @@ include("IO_functions.jl")
 using .IO_functions
 export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File
 export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames
-export clean_directory, changefolder
+export clean_directory, changefolder, read_phase_diagram
 
 # Functions 
 include("Run.jl")

src/LaMEM_ModelGeneration/DefaultParams.jl

@@ -0,0 +1,27 @@
+# This sets default parameters, depending  on various combinations of input patameters
+
+"""
+    model = UpdateDefaultParameters(model::Model)
+
+This updates the default parameters depending on some of the input parameters.
+If you activate passive tracers, for example, it will also activate output for that 
+"""
+function UpdateDefaultParameters(model::Model)
+
+    # If PhaseTransitions are defined, we generally want this to be activated in computations
+    if !isempty(model.Materials.PhaseTransitions)
+        model.SolutionParams.Phasetrans = 1
+    end
+
+    # If PassiveTracers are defined, we generally want this to be visualized as well:
+    if model.PassiveTracers.Passive_Tracer==1
+
+        model.Output.out_ptr=1
+        model.Output.out_ptr_ID=1
+        model.Output.out_ptr_phase=1
+        model.Output.out_ptr_Pressure=1
+        model.Output.out_ptr_Temperature=1
+    end
+
+    return model
+end

src/LaMEM_ModelGeneration/Model.jl

@@ -5,7 +5,7 @@ using LaMEM.Run.LaMEM_jll
 
 export Model, Write_LaMEM_InputFile, create_initialsetup, run_lamem
 
-""";
+"""
     Model
 
 Structure that holds all the information to create a LaMEM input file
@@ -87,6 +87,8 @@ function Model(;
                 SolutionParams, Solver, ModelSetup, Output, PassiveTracers, Materials)
 end
 
+include("DefaultParams.jl")             # main LaMEM_Model
+
 """
     Model(args...)
 
@@ -122,7 +124,10 @@ function Model(args...)
     end 
     args_tuple = NamedTuple{Symbol.(names_strip)}(args)
 
-    return Model(; args_tuple...)
+    model = Model(; args_tuple...)
+    model = UpdateDefaultParameters(model)
+
+    return model
 end
 
 # Show brief overview of Model
@@ -154,19 +159,7 @@ function Write_LaMEM_InputFile(d::Model, fname::String="input.dat"; dir=pwd())
         # If we want to write an input file 
         Write_Paraview(CartData(d.Grid.Grid, (Phases=d.Grid.Phases,Temp=d.Grid.Temp)),"Model3D")
     end
-    if !isempty(d.Materials.PhaseTransitions)
-        # If PhaseTransitions are defined, we generally want this to be activated in computations
-        d.SolutionParams.Phasetrans = 1
-    end
-    if d.PassiveTracers.Passive_Tracer==1
-        # If PassiveTracers are defined, we generally want this to be visualized as well:
-        d.Output.out_ptr=1
-        d.Output.out_ptr_ID=1
-        d.Output.out_ptr_phase=1
-        d.Output.out_ptr_Pressure=1
-        d.Output.out_ptr_Temperature=1
-    end
-
+
     io = open(fname,"w")
 
     Write_LaMEM_InputFile(io, d.Scaling)

src/LaMEM_ModelGeneration/SolutionParams.jl

@@ -182,8 +182,9 @@ function Write_LaMEM_InputFile(io, d::SolutionParams)
     for f in fields
         if getfield(d,f) != getfield(Reference,f) ||
             (f == :eta_ref) ||
-            (f == :gravity)
-
+            (f == :gravity) ||
+            (f == :act_p_shift)
+
             # only print if value differs from reference value
             name = rpad(String(f),15)
             data = getfield(d,f) 

src/PlotsExt.jl

@@ -1,10 +1,10 @@
 # Plotting extensions, that are only loaded when GLMakie is available
 println("Adding Plots.jl plotting extensions for LaMEM")
 
-using LaMEM, GeophysicalModelGenerator, Statistics
+using LaMEM, GeophysicalModelGenerator, Statistics, DelimitedFiles
 using .Plots
 import .Plots: heatmap
-export plot_topo, plot_cross_section
+export plot_topo, plot_cross_section, plot_phasediagram
 
 """
     plot_cross_section(model::Union{Model,CartData}, args...; field::Symbol=:phase,  dim=1, x=nothing, y=nothing, z=nothing, aspect_ratio::Union{Real, Symbol}=:equal)
@@ -67,4 +67,37 @@ function plot_topo(topo::CartData; kwargs...)
                 kwargs...)
 
     return hm
-end
+end
+
+
+
+
+""" 
+    hm = plot_phasediagram(name::String="Rhyolite", field=:ρ;  colormap=:batlow, kwargs...)
+
+This creates a plot of a LaMEM phase diagram (computed by MAGEMin or Perple_X).
+Typical available `fields`:
+- `:ρ_solid` (solid density)
+- `:ρ_melt` (melt density)
+- `:ρ` (average density)
+- `:ϕ` (melt fraction)
+
+"""
+function plot_phasediagram(name::String="Rhyolite", field=:ρ;  colormap=:batlow, kwargs...)
+    # read phase diagram
+    PD = read_phase_diagram(name)
+
+    T_C     = PD.T_K[:,1] .- 273.15   
+    P_kbar  = PD.P_bar[1,:] ./ 1e3
+    val     = PD[field]     # field
+
+    # plot
+    hm = heatmap(T_C, P_kbar, val'; 
+                xlabel="T [Celcius]",
+                ylabel="Pressure [kbar]",
+                title="Phase diagram [$field]",
+                colormap=colormap,
+                kwargs...)
+
+    return hm
+end

src/utils_IO.jl

@@ -1,6 +1,6 @@
-using Glob
+using Glob, DelimitedFiles
 
-export clean_directory, changefolder
+export clean_directory, changefolder, read_phase_diagram
 
 """ 
     clean_directory(DirName)
@@ -78,3 +78,37 @@ function changefolder()
         error("This only works on windows and mac")
     end
 end
+
+"""
+    out = read_phase_diagram(name::String)
+
+Reads a phase diagram from a file `name` and returns a NamedTuple with temperature `T`, pressure `P`, melt density `ρ_melt`, solid density `ρ_solid`, density `ρ` and melt fraction `ϕ`
+"""
+function read_phase_diagram(name::String)
+
+    f = open(name)
+
+    # Read dimensions
+    for i = 1:49; readline(f); end
+    minT = parse(Float64,   readline(f))
+    ΔT   = parse(Float64,   readline(f))
+    nT   = parse(Int64,     readline(f))
+    minP = parse(Float64,   readline(f))
+    ΔP   = parse(Float64,   readline(f))
+    nP   = parse(Int64,     readline(f))
+    close(f)
+
+    data = readdlm(name, skipstart=55);     # read numerical data
+
+    # reshape
+    ρ_melt  = reshape(data[:,1],(nT,nP));
+    ϕ       = reshape(data[:,2],(nT,nP));
+    ρ_solid = reshape(data[:,3],(nT,nP));
+    T_K     = reshape(data[:,4],(nT,nP));
+    P_bar   = reshape(data[:,5],(nT,nP));
+    ρ       = ρ_melt.*ϕ .+ ρ_solid.*(1.0 .- ϕ);
+
+
+    return (;T_K,P_bar,ρ_melt,ρ_solid,ϕ, ρ)
+end
+

test/CreateMarkers_Subduction_Linear_FreeSlip_parallel.jl

@@ -44,6 +44,6 @@ Model3D     =   CartData(Grid_LaMEM, (Phases=Phases,Temp=Temp))   # Create LaMEM
 Write_Paraview(Model3D,"LaMEM_ModelSetup")                  # Save model to paraview   (load with opening LaMEM_ModelSetup.vts in paraview)  
 
 # Save LaMEM markers
-dir =   joinpath(pkg_dir,"test");
+dir =   joinpath(pkg_dir,"test","markers");
 Save_LaMEMMarkersParallel(Model3D, directory=dir)                          # Create LaMEM marker input on 1 core