Skip to content

Commit

Permalink
update GMG version
Browse files Browse the repository at this point in the history
  • Loading branch information
@boriskaus
boriskaus committed Apr 6, 2024
1 parent 4e9b097 commit 5e50cc3
Show file tree
Hide file tree
Showing 13 changed files with 56 additions and 58 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
[compat]
CUDA = "4 - 5"
CairoMakie = "0.8 - 0.20"
GeophysicalModelGenerator = "0.3 - 0.6"
GeophysicalModelGenerator = "0.7"
GeoParams = "0.5"
Interpolations = "0.13 - 0.15"
JLD2 = "0.4"
Expand Down
8 changes: 6 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,14 @@ This easy to use and versatile package simulates the thermal evolution of magmat
Below we give a number of example scripts that show how it can be used to simulate a number of scenarios.

## Contents
- [MagmaThermoKinematics.jl](#magmathermokinematicsjl)
- [Contents](#contents)
- [100-lines 2D example](#100-lines-2d-example)
- [100-lines 3D example](#100-lines-3d-example)
- [Dependencies](#dependencies)
- [Installation](#installation)
- [Ongoing development](#ongoing-development)
- [Related work](#related-work)

## 100-lines 2D example
A simple example that simulates the emplacement of dikes within the crust over a period of 10'000 years is shown below.
Expand Down Expand Up @@ -41,7 +45,7 @@ Num = Numeric_params(verbose=false) # No
MatParam = (
SetMaterialParams(Name="Rock", Phase=1,
Density = ConstantDensity(ρ=2800kg/m^3),
HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
Conductivity = ConstantConductivity(k=1.5Watt/K/m),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
Melting = MeltingParam_Caricchi()),
Expand Down Expand Up @@ -166,7 +170,7 @@ using WriteVTK
MatParam = (
SetMaterialParams(Name="Rock", Phase=1,
Density = ConstantDensity(ρ=2800kg/m^3),
HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
Conductivity = ConstantConductivity(k=1.5Watt/K/m),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
Melting = MeltingParam_Caricchi()),
Expand Down
2 changes: 1 addition & 1 deletion examples/Example2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ using Plots
MatParam = (
SetMaterialParams(Name="Rock", Phase=1,
Density = ConstantDensity(ρ=2800kg/m^3),
HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
Conductivity = ConstantConductivity(k=1.5Watt/K/m),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
Melting = MeltingParam_Caricchi()),
Expand Down
10 changes: 5 additions & 5 deletions examples/Example2D_ZASSy.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@
# It includes comparisons with 2D simulations done by the Geneva (Gregor Weber, Luca Caricchi) & UCLA (Oscar Lovera) Tracers_SimParams
#
#
const USE_GPU=true;
const USE_GPU=false;
using MagmaThermoKinematics
if USE_GPU
environment!(:gpu, Float64, 2) # initialize parallel stencil in 2D
CUDA.device!(1) # select the GPU you use (starts @ zero)
# CUDA.device!(1) # select the GPU you use (starts @ zero)
else
environment!(:cpu, Float64, 2) # initialize parallel stencil in 2D
end
Expand Down Expand Up @@ -446,7 +446,7 @@ if 1==0
# Conductivity = ConstantConductivity(k=3.3Watt/K/m), # in case we use constant k
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
#Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting
# Melting = MeltingParam_Caricchi()), # Caricchi melting

Expand Down Expand Up @@ -665,7 +665,7 @@ if 1==1
Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = T_HeatCapacity_Whittington(), # T-dependent cp
# Conductivity = ConstantConductivity(k=3.3Watt/K/m), # in case we use constant k
# HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
# HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = MeltingParam_Assimilation() # Quadratic parameterization as in Tierney et al.
#Melting = MeltingParam_Caricchi()
),
Expand All @@ -676,7 +676,7 @@ if 1==1
Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = T_HeatCapacity_Whittington(), # T-dependent cp
#Conductivity = ConstantConductivity(k=3.3Watt/K/m), # in case we use constant k
#HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
#HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_Quadratic(T_s=(700+273.15)K, T_l=(1100+273.15)K))
#Melting = MeltingParam_Caricchi()
)
Expand Down
2 changes: 1 addition & 1 deletion examples/Example3D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ using WriteVTK
MatParam = (
SetMaterialParams(Name="Rock", Phase=1,
Density = ConstantDensity(ρ=2800kg/m^3),
HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
Conductivity = ConstantConductivity(k=1.5Watt/K/m),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
Melting = MeltingParam_Caricchi()),
Expand Down
2 changes: 1 addition & 1 deletion examples/Example3D_v2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ using WriteVTK
MatParam = (
SetMaterialParams(Name="Rock", Phase=1,
Density = ConstantDensity(ρ=2800kg/m^3),
HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
Conductivity = ConstantConductivity(k=1.5Watt/K/m),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
Melting = MeltingParam_Caricchi()),
Expand Down
2 changes: 1 addition & 1 deletion examples/MTK_GMG_2D_example1.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ function MTK_GMG.MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::Numeri
# open pvd file if requested
if Num.Output_VTK
name = joinpath(Num.SimName,Num.SimName*".pvd")
Num.pvd = Movie_Paraview(name=name, Initialize=true);
Num.pvd = movie_paraview(name=name, Initialize=true);
end

return nothing
Expand Down
40 changes: 17 additions & 23 deletions examples/MTK_GMG_2D_example2.jl
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,4 @@
# This is a second example that shows how to use MTK in combination with the GeophysicalModelGenerator
# Compared to example 1, we show a few more features:
# - How to define a custom structure with temporal values & use it in the code
# - How to generate a model setup using GMG
# - How to overwrite some of the default functions to customize the simulation
# - How to create paraview files

# Unzen setup
const USE_GPU=false;
using MagmaThermoKinematics
if USE_GPU
Expand All @@ -18,33 +12,34 @@ using Plots
using Random
using GeophysicalModelGenerator

# Test setup
println("Example 2 of the MTK - GMG integration")

# Model setup
println(" --- Generating Setup --- ")

# Topography and project it.

# NOTE: The first time you do this, please set this to true, which will download the topography data from the internet and save it in a file
if false
using GMT, Statistics
Topo = ImportTopo(lon = [130.0, 130.5], lat=[32.55, 32.90], file="@earth_relief_03s.grd")
proj = ProjectionPoint(; Lat=mean(Topo.lat.val), Lon=mean(Topo.lon.val))
Topo_cart = Convert2CartData(Topo, proj)
save_GMG("examples/Topo_cart", Topo_cart)
using GMT, Statistics
Topo = ImportTopo(lon = [130.0, 130.5], lat=[32.55, 32.90], file="@earth_relief_03s.grd")
proj = ProjectionPoint(; Lat=mean(Topo.lat.val), Lon=mean(Topo.lon.val))
Topo_cart = Convert2CartData(Topo, proj)
Xt,Yt,Zt = xyz_grid(-23:.1:23,-19:.1:19,0)
Topo_cart = ProjectCartData(CartData(Xt,Yt,Zt,(Zt=Zt,)), Topo, proj)

save_GMG("Topo_cart", Topo_cart)
end
Topo_cart = load_GMG("Topo_cart")

# Create 3D grid of the region
X,Y,Z = XYZGrid(-23:.1:23,-19:.1:19,-20:.1:5)
X,Y,Z = xyz_grid(-23:.1:23,-19:.1:19,-20:.1:5)
Data_set3D = CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X)))); # 3D dataset

# Create 2D cross-section
Nx = 135*6; # resolution in x
Nz = 135*4;
Data_2D = CrossSection(Data_set3D, Start=(-20,4), End=(20,4), dims=(Nx, Nz))
Data_2D = AddField(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
Data_2D = AddField(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
Data_2D = addfield(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
Data_2D = addfield(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))

# Intersect with topography
Below = BelowSurface(Data_2D, Topo_cart)
Expand All @@ -63,7 +58,6 @@ x_c, z_c, r = -10, -15, 2.5
Volume = 4/3*pi*r^3 # equivalent 3D volume of the anomaly [km^3]
@views Data_2D.fields.Temp[(Data_2D.x.val .- x_c).^2 .+ (Data_2D.z.val .- z_c).^2 .< r^2] .= 800.0


println(" --- Performing MTK models --- ")

# Overwrite some of the default functions
Expand Down Expand Up @@ -167,24 +161,24 @@ MatParam = (SetMaterialParams(Name="Air", Phase=0,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=0.0J/kg),
Conductivity = ConstantConductivity(k=3Watt/K/m), # in case we use constant k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting
SetMaterialParams(Name="Crust", Phase=1,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting
SetMaterialParams(Name="Mantle", Phase=2,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K)),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K)),
SetMaterialParams(Name="Dikes", Phase=3,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())) # Marxer & Ulmer melting
)

Expand Down
14 changes: 7 additions & 7 deletions src/MTK_GMG.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@ function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParam
# Open pvd file if requested
if Num.Output_VTK
name = joinpath(Num.SimName,Num.SimName*".pvd")
Num.pvd = Movie_Paraview(name=name, Initialize=true);
Num.pvd = movie_paraview(name=name, Initialize=true);
end

return nothing
Expand Down Expand Up @@ -201,7 +201,7 @@ function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParam
# open pvd file if requested
if Num.Output_VTK
name = joinpath(Num.SimName,Num.SimName*".pvd")
Num.pvd = Movie_Paraview(name=name, Initialize=true);
Num.pvd = movie_paraview(name=name, Initialize=true);
end

return nothing
Expand All @@ -215,7 +215,7 @@ Finalize model run
"""
function MTK_finalize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::Union{Nothing,CartData})
if Num.Output_VTK & !isnothing(Num.pvd)
Movie_Paraview(pvd=Num.pvd, Finalize=true)
movie_paraview(pvd=Num.pvd, Finalize=true)
end

return nothing
Expand Down Expand Up @@ -270,9 +270,9 @@ function MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArra
Data_set3D = CartData_input
else
if length(Grid.coord1D)==3
X,Y,Z = XYZGrid(Grid.coord1D...)
X,Y,Z = xyz_grid(Grid.coord1D...)
elseif length(Grid.coord1D)==2
X,Y,Z = XYZGrid(Grid.coord1D[1], 0, Grid.coord1D[2])
X,Y,Z = xyz_grid(Grid.coord1D[1], 0, Grid.coord1D[2])
end
Data_set3D = CartData(X/1e3,Y/1e3,Z/1e3, (Z=Z,))
end
Expand All @@ -282,7 +282,7 @@ function MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArra
Data_set3D = add_data_CartData(Data_set3D, "MeltFraction", Array(Arrays.ϕ));

# Save output to CartData
Num.pvd = Write_Paraview(Data_set3D, name, pvd=Num.pvd,time=Num.time/SecYear/1e3);
Num.pvd = write_paraview(Data_set3D, name, pvd=Num.pvd,time=Num.time/SecYear/1e3);
end
end
return nothing
Expand All @@ -304,7 +304,7 @@ function add_data_CartData(d::CartData, name::String, data::Array)
else
a = data
end
d = AddField(d, name, a)
d = addfield(d, name, a)
return d
end

Expand Down
2 changes: 1 addition & 1 deletion src/MTK_GMG_2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ There are a few functions that you can overwrite in your user code to customize
if !isnothing(CartData_input)

if !hasfield(typeof(CartData_input.fields),:FlatCrossSection)
error("You should add a Field :FlatCrossSection to your data structure with Data_Cross = AddField(Data_Cross,\"FlatCrossSection\", FlattenCrossSection(Data_Cross))")
error("You should add a Field :FlatCrossSection to your data structure with Data_Cross = addfield(Data_Cross,\"FlatCrossSection\", FlattenCrossSection(Data_Cross))")
end

Num = MTK_GMG.Setup_Model_CartData(CartData_input, Num, Mat_tup)
Expand Down
16 changes: 8 additions & 8 deletions test/test_MTK_GMG_2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ MatParam = (SetMaterialParams(Name="Rock & partial melt", Phase=1,
# Conductivity = ConstantConductivity(k=3.3Watt/K/m), # in case we use constant k
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
#Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting
# Melting = MeltingParam_Caricchi()), # Caricchi melting
# add more parameters here, in case you have >1 phase in the model
Expand All @@ -75,15 +75,15 @@ Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_
Topo_cart = load_GMG("../examples/Topo_cart") # Note: Laacher seee is around [10,20]

# Create 3D grid of the region
X,Y,Z = XYZGrid(-23:.1:23,-19:.1:19,-20:.1:5)
X,Y,Z = xyz_grid(-23:.1:23,-19:.1:19,-20:.1:5)
Data_set3D = CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X)))); # 3D dataset

# Create 2D cross-section
Nx = Num.Nx; # resolution in x
Nz = Num.Nz;
Data_2D = CrossSection(Data_set3D, Start=(-20,4), End=(20,4), dims=(Nx, Nz))
Data_2D = AddField(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
Data_2D = AddField(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
Data_2D = addfield(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
Data_2D = addfield(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))

# Intersect with topography
Below = BelowSurface(Data_2D, Topo_cart)
Expand Down Expand Up @@ -141,30 +141,30 @@ MatParam = (SetMaterialParams(Name="Air", Phase=0,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=0.0J/kg),
Conductivity = ConstantConductivity(k=3Watt/K/m), # in case we use constant k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting

SetMaterialParams(Name="Crust", Phase=1,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
#Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())), # Marxer & Ulmer melting

SetMaterialParams(Name="Mantle", Phase=2,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K)),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K)),

SetMaterialParams(Name="Dikes", Phase=3,
Density = ConstantDensity(ρ=2700kg/m^3),
LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
# Conductivity = ConstantConductivity(k=3.3Watt/K/m), # in case we use constant k
Conductivity = T_Conductivity_Whittington_parameterised(), # T-dependent k
#Conductivity = T_Conductivity_Whittington(), # T-dependent k
HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
Melting = SmoothMelting(MeltingParam_4thOrder())) # Marxer & Ulmer melting

)
Expand Down
Loading

0 comments on commit 5e50cc3

Please sign in to comment.