diff --git a/Project.toml b/Project.toml index 71cbcc2..0fd1abb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,13 @@ name = "MagmaThermoKinematics" uuid = "e7fe3d1a-f102-40b5-a66b-64b14986c4c8" authors = ["Boris Kaus "] -version = "0.4.2" +version = "0.5.0" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38" +GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -25,13 +25,13 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] CUDA = "4 - 5" -CairoMakie = "0.8 - 0.10" -GLMakie = "0.6 - 0.8" -Interpolations = "0.13,0.14" -ParallelStencil = "0.8 - 0.9" +CairoMakie = "0.8 - 0.20" +GeophysicalModelGenerator = "0.3 - 0.5" GeoParams = "0.5" +Interpolations = "0.13 - 0.15" JLD2 = "0.4" MAT = "0.10" +ParallelStencil = "0.8 - 0.10" Parameters = "0.12" Plots = "1.0" Reexport = "1" diff --git a/examples/Example2D.jl b/examples/Example2D.jl index 0775510..85dee14 100644 --- a/examples/Example2D.jl +++ b/examples/Example2D.jl @@ -1,7 +1,9 @@ using MagmaThermoKinematics const USE_GPU=false; -if USE_GPU environment!(:gpu, Float64, 2) # initialize in 2D on GPU -else environment!(:cpu, Float64, 2) # initialize in 3D on CPU +if USE_GPU + environment!(:gpu, Float64, 2) # initialize in 2D on GPU +else + environment!(:cpu, Float64, 2) # initialize in 3D on CPU end using MagmaThermoKinematics.Diffusion2D using Plots @@ -40,8 +42,11 @@ using Plots # CPU buffers Tnew_cpu = Matrix{Float64}(undef, Grid.N...) Phi_melt_cpu = similar(Tnew_cpu) - if USE_GPU; Phases = CUDA.ones(Int64,Grid.N...) - else Phases = ones(Int64,Grid.N...) end + if USE_GPU; + Phases = CUDA.ones(Int64,Grid.N...) + else + Phases = ones(Int64,Grid.N...) + end @parallel (1:Nx, 1:Nz) GridArray!(Arrays.X, Arrays.Z, Grid.coord1D[1], Grid.coord1D[2]) Tracers = StructArray{Tracer}(undef, 1) # Initialize tracers @@ -57,8 +62,11 @@ using Plots if floor(time/InjectionInterval)> dike_inj # Add new dike every X years dike_inj = floor(time/InjectionInterval) # Keeps track on what was injected already cen = (Grid.max .+ Grid.min)./2 .+ rand(-0.5:1e-3:0.5, 2).*[W_ran;H_ran]; # Randomly vary center of dike - if cen[end]<-12e3; Angle_rand = rand( 80.0:0.1:100.0) # Orientation: near-vertical @ depth - else Angle_rand = rand(-10.0:0.1:10.0); end # Orientation: near-vertical @ shallower depth + if cen[end]<-12e3; + Angle_rand = rand( 80.0:0.1:100.0) # Orientation: near-vertical @ depth + else + Angle_rand = rand(-10.0:0.1:10.0); + end # Orientation: near-vertical @ shallower depth dike = Dike(dike, Center=cen[:],Angle=[Angle_rand]); # Specify dike with random location/angle but fixed size/T Tnew_cpu .= Array(Arrays.T) Tracers, Tnew_cpu, Vol = InjectDike(Tracers, Tnew_cpu, Grid.coord1D, dike, nTr_dike); # Add dike, move hostrocks diff --git a/examples/Example3D.jl b/examples/Example3D.jl index d0ba588..9297acc 100644 --- a/examples/Example3D.jl +++ b/examples/Example3D.jl @@ -1,7 +1,9 @@ using MagmaThermoKinematics const USE_GPU=false; -if USE_GPU environment!(:gpu, Float64, 3) # initialize parallel stencil in 2D -else environment!(:cpu, Float64, 3) # initialize parallel stencil in 2D +if USE_GPU + environment!(:gpu, Float64, 3) # initialize parallel stencil in 2D +else + environment!(:cpu, Float64, 3) # initialize parallel stencil in 2D end using MagmaThermoKinematics.Diffusion3D using Plots @@ -41,8 +43,11 @@ using WriteVTK # CPU buffers Tnew_cpu = zeros(Float64, Grid.N...) Phi_melt_cpu = similar(Tnew_cpu) - if USE_GPU; Phases = CUDA.ones(Int64,Grid.N...) - else Phases = ones(Int64,Grid.N...) end + if USE_GPU; + Phases = CUDA.ones(Int64,Grid.N...) + else + Phases = ones(Int64,Grid.N...) + end @parallel (1:Nx,1:Ny,1:Nz) GridArray!(Arrays.X,Arrays.Y,Arrays.Z, Grid.coord1D[1], Grid.coord1D[2], Grid.coord1D[3]) Tracers = StructArray{Tracer}(undef, 1) # Initialize tracers @@ -50,7 +55,11 @@ using WriteVTK Arrays.T .= -Arrays.Z.*GeoT; # Initial (linear) temperature profile # Preparation of VTK/Paraview output - if isdir("viz3D_out")==false mkdir("viz3D_out") end; loadpath = "./viz3D_out/"; pvd = paraview_collection("Example3D"); + if isdir("viz3D_out")==false + mkdir("viz3D_out") + end; + loadpath = "./viz3D_out/"; + pvd = paraview_collection("Example3D"); time, dike_inj, InjectVol, Time_vec,Melt_Time = 0.0, 0.0, 0.0,zeros(nt,1),zeros(nt,1); for it = 1:nt # Time loop @@ -58,8 +67,11 @@ using WriteVTK if floor(time/InjectionInterval)> dike_inj # Add new dike every X years dike_inj = floor(time/InjectionInterval) # Keeps track on what was injected already cen = (Grid.max .+ Grid.min)./2 .+ rand(-0.5:1e-3:0.5, 3).*[W_ran;W_ran;H_ran]; # Randomly vary center of dike - if cen[end]<-12e3; Angle_rand = [rand(80.0:0.1:100.0); rand(0:360)] # Dikes at depth - else Angle_rand = [rand(-10.0:0.1:10.0); rand(0:360)] end # Sills at shallower depth + if cen[end]<-12e3; + Angle_rand = [rand(80.0:0.1:100.0); rand(0:360)] # Dikes at depth + else + Angle_rand = [rand(-10.0:0.1:10.0); rand(0:360)] + end # Sills at shallower depth dike = Dike(dike, Center=cen[:],Angle=Angle_rand); # Specify dike with random location/angle but fixed size/T Tnew_cpu .= Array(Arrays.T) Tracers, Tnew_cpu, Vol = InjectDike(Tracers, Tnew_cpu, Grid.coord1D, dike, nTr_dike); # Add dike, move hostrocks diff --git a/examples/Example3D_v2.jl b/examples/Example3D_v2.jl new file mode 100644 index 0000000..3d6b1fa --- /dev/null +++ b/examples/Example3D_v2.jl @@ -0,0 +1,107 @@ +using MagmaThermoKinematics +const USE_GPU=false; +if USE_GPU + environment!(:gpu, Float64, 3) # initialize parallel stencil in 2D +else + environment!(:cpu, Float64, 3) # initialize parallel stencil in 2D +end +using MagmaThermoKinematics.Diffusion3D +using Plots +using WriteVTK + +#------------------------------------------------------------------------------------------ +@views function MainCode_3D(); + Nx,Ny,Nz = 250,250,250 + Grid = CreateGrid(size=(Nx,Ny,Nz), extent=(30e3, 30e3, 30e3)) # grid points & domain size + Num = Numeric_params(verbose=false) # Nonlinear solver options + + # Set material parameters + MatParam = ( + SetMaterialParams(Name="Rock", Phase=1, + Density = ConstantDensity(ρ=2800kg/m^3), + HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K), + Conductivity = ConstantConductivity(k=1.5Watt/K/m), + LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg), + Melting = MeltingParam_Caricchi()), + ) + + GeoT = 20.0/1e3; # Geothermal gradient [K/km] + W_in, H_in = 5e3, 0.5e3; # Width and thickness of dike + T_in = 900; # Intrusion temperature + InjectionInterval = 0.1kyr; # Inject a new dike every X kyrs + maxTime = 15kyr; # Maximum simulation time in kyrs + H_ran, W_ran = Grid.L[1]*0.3,Grid.L[3]*0.4;# Size of domain in which we randomly place dikes and range of angles + DikeType = "ElasticDike" # Type to be injected ("ElasticDike","SquareDike") + κ = 1.2/(2800*1050); # thermal diffusivity + dt = minimum(Grid.Δ.^2)/κ/10; # stable timestep (required for explicit FD) + nt = floor(Int64,maxTime/dt); # number of required timesteps + nTr_dike = 300; # number of tracers inserted per dike + + # Array initializations + Arrays = CreateArrays(Dict( (Nx, Ny, Nz)=>(T=0,T_K=0, T_it_old=0, K=1.5, Rho=2800, Cp=1050, Tnew=0, Hr=0, Hl=0, Kc=1, P=0, X=0, Y=0, Z=0, ϕₒ=0, ϕ=0, dϕdT=0), + (Nx-1,Ny,Nz)=>(qx=0,Kx=0), (Nx, Ny-1, Nz)=>(qy=0,Ky=0 ) , (Nx, Ny, Nz-1)=>(qz=0,Kz=0 ) )) + # CPU buffers + Tnew_cpu = zeros(Float64, Grid.N...) + Phi_melt_cpu = similar(Tnew_cpu) + if USE_GPU; + Phases = CUDA.ones(Int64,Grid.N...) + else + Phases = ones(Int64,Grid.N...) + end + + @parallel (1:Nx,1:Ny,1:Nz) GridArray!(Arrays.X,Arrays.Y,Arrays.Z, Grid.coord1D[1], Grid.coord1D[2], Grid.coord1D[3]) + Tracers = StructArray{Tracer}(undef, 1) # Initialize tracers + dike = Dike(W=W_in,H=H_in,Type=DikeType,T=T_in); # "Reference" dike with given thickness,radius and T + Arrays.T .= -Arrays.Z.*GeoT; # Initial (linear) temperature profile + + # Preparation of VTK/Paraview output + if isdir("viz3D_out")==false + mkdir("viz3D_out") + end; + loadpath = "./viz3D_out/"; + pvd = paraview_collection("Example3D"); + + time, dike_inj, InjectVol, Time_vec,Melt_Time = 0.0, 0.0, 0.0,zeros(nt,1),zeros(nt,1); + for it = 1:nt # Time loop + + if floor(time/InjectionInterval)> dike_inj # Add new dike every X years + dike_inj = floor(time/InjectionInterval) # Keeps track on what was injected already + cen = (Grid.max .+ Grid.min)./2 .+ rand(-0.5:1e-3:0.5, 3).*[W_ran;W_ran;H_ran]; # Randomly vary center of dike + if cen[end]<-12e3; + Angle_rand = [rand(80.0:0.1:100.0); rand(0:360)] # Dikes at depth + else + Angle_rand = [rand(-10.0:0.1:10.0); rand(0:360)] + end # Sills at shallower depth + dike = Dike(dike, Center=cen[:],Angle=Angle_rand); # Specify dike with random location/angle but fixed size/T + Tnew_cpu .= Array(Arrays.T) + Tracers, Tnew_cpu, Vol = InjectDike(Tracers, Tnew_cpu, Grid.coord1D, dike, nTr_dike); # Add dike, move hostrocks + Arrays.T .= Data.Array(Tnew_cpu) + InjectVol += Vol # Keep track of injected volume + println("Added new dike; total injected magma volume = $(round(InjectVol/km³,digits=2)) km³; rate Q=$(round(InjectVol/(time),digits=2)) m³/s") + end + + Nonlinear_Diffusion_step_3D!(Arrays, MatParam, Phases, Grid, dt, Num) # Perform a nonlinear diffusion step + + copy_arrays_GPU2CPU!(Tnew_cpu, Phi_melt_cpu, Arrays.Tnew, Arrays.ϕ) # Copy arrays to CPU to update properties + UpdateTracers_T_ϕ!(Tracers, Grid.coord1D, Tnew_cpu, Phi_melt_cpu); # Update info on tracers + + @parallel assign!(Arrays.T, Arrays.Tnew) + @parallel assign!(Arrays.Tnew, Arrays.T) # Update temperature + time = time + dt; # Keep track of evolved time + Melt_Time[it] = sum(Arrays.ϕ)/prod(Grid.N) # Melt fraction in crust + Time_vec[it] = time; # Vector with time + println(" Timestep $it = $(round(time/kyr*100)/100) kyrs") + + if mod(it,20)==0 # Visualisation + x,y,z = Grid.coord1D[1], Grid.coord1D[2], Grid.coord1D[3] + vtkfile = vtk_grid("./viz3D_out/ex3D_$(Int32(it+1e4))", Vector(x/1e3), Vector(y/1e3), Vector(z/1e3)) # 3-D VTK file + vtkfile["Temperature"] = Array(Arrays.T); vtkfile["MeltFraction"] = Array(Arrays.ϕ); # Store fields in file + outfiles = vtk_save(vtkfile); pvd[time/kyr] = vtkfile # Save file & update pvd file + end + end + vtk_save(pvd) + return Time_vec, Melt_Time, Tracers, Grid, Arrays; +end # end of main function + +Time_vec, Melt_Time, Tracers, Grid, Arrays = MainCode_3D(); # start the main code +plot(Time_vec/kyr, Melt_Time, xlabel="Time [kyrs]", ylabel="Fraction of crust that is molten", label=:none); png("Time_vs_Melt_Example2D") # Create plot \ No newline at end of file diff --git a/examples/MTK_GMG_2D_example1.jl b/examples/MTK_GMG_2D_example1.jl new file mode 100644 index 0000000..669c200 --- /dev/null +++ b/examples/MTK_GMG_2D_example1.jl @@ -0,0 +1,129 @@ +# This is a first example of how to use MagmaThermoKinematics with a real setup which can be customized with user-defined functions, +# for example for plotting or printing output. + +const USE_GPU=false; +using MagmaThermoKinematics +if USE_GPU + environment!(:gpu, Float64, 2) # initialize parallel stencil in 2D +else + environment!(:cpu, Float64, 2) # initialize parallel stencil in 2D +end +using MagmaThermoKinematics.Diffusion2D +using MagmaThermoKinematics.MTK_GMG # Allow overwriting user routines +using MagmaThermoKinematics.GeophysicalModelGenerator +using Plots # plots +using Random, GeoParams + +Random.seed!(1234); # such that we can reproduce results + +# Test setup +println("Example 1 of the MTK - GMG integration") + +# Overwrite some functions +function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + println("$(Num.it), Time=$(round(Num.time/Num.SecYear)) yrs; max(T) = $(round(maximum(Arrays.Tnew)))") + return nothing +end + +@static if USE_GPU + function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + println("$(Num.it), Time=$(round(Num.time/Num.SecYear)) yrs; max(T) = $(round(maximum(Arrays.Tnew)))") + return nothing + end +else + function MTK_GMG.MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + if mod(Num.it,Num.CreateFig_steps)==0 + x_1d = Grid.coord1D[1]/1e3; + z_1d = Grid.coord1D[2]/1e3; + temp_data = Array(Arrays.Tnew)' + ϕ_data = Array(Arrays.ϕ)' + phase_data = Array(Arrays.Phases)' + + t = Num.time/SecYear; + p = plot(layout=grid(1,2) ) + + Plots.heatmap!(p[1],x_1d, z_1d, temp_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Temperature, t=$(round(t)) yrs", aspect_ratio=:equal, ylimits=(-20,0)) + # Plots.heatmap!(p[2],x_1d, z_1d, ϕ_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Melt fraction", clims=(0,1), aspect_ratio=:equal, ylimits=(-20,0)) + Plots.heatmap!(p[2],x_1d, z_1d, phase_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Phases", aspect_ratio=:equal, ylimits=(-20,0)) + + # p = plot(ps, layout=(1,2)) + display(p) + end + return nothing + end +end + +""" + MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArrays, Dikes::DikeParameters) + +Initialize temperature and phases of the grid +""" +function MTK_GMG.MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters) + # Initalize T + Arrays.T_init .= @. Num.Tsurface_Celcius - Arrays.Z*Num.Geotherm; # Initial (linear) temperature profile + + # Initialize Phases + @views Arrays.Phases[Arrays.Z .> -5000] .= 0; + Arrays.Phases_init .= Arrays.Phases; # Initialize all as rock + + # open pvd file if requested + if Num.Output_VTK + name = joinpath(Num.SimName,Num.SimName*".pvd") + Num.pvd = Movie_Paraview(name=name, Initialize=true); + end + + return nothing +end + + +# Define numerical parameters +Num = NumParam( Nx = 135*2, + Nz = 135*2, + SimName = "Test1", + maxTime_Myrs = 0.005, + fac_dt = 0.2, + ω = 0.5, + CreateFig_steps = 20, + SaveOutput_steps = 100, + USE_GPU = USE_GPU, + AddRandomSills = true, + RandomSills_timestep = 5); + +Dike_params = DikeParam(Type = "ElasticDike", + InjectionInterval_year = 1000, + W_in = 5e3, + H_in = 250, + nTr_dike = 300*4, + DikePhase = 2, + T_in_Celsius = 1000, + SillsAbove = -12e3 # below this we have dikes; above sills + ) + +MatParam = (SetMaterialParams(Name="Host rock 1", Phase=0, + Density = ConstantDensity(ρ=2700kg/m^3), # used in the parameterisation of Whittington + LatentHeat = ConstantLatentHeat(Q_L=2.55e5J/kg), + RadioactiveHeat = ExpDepthDependentRadioactiveHeat(H_0=0e-7Watt/m^3), + Conductivity = T_Conductivity_Whittington(), # T-dependent k + HeatCapacity = T_HeatCapacity_Whittington(), # T-dependent cp + Melting = MeltingParam_Assimilation() # Quadratic parameterization as in Tierney et al. + ), + SetMaterialParams(Name="Host rock", Phase=1, + Density = ConstantDensity(ρ=2700kg/m^3), # used in the parameterisation of Whittington + LatentHeat = ConstantLatentHeat(Q_L=2.55e5J/kg), + RadioactiveHeat = ExpDepthDependentRadioactiveHeat(H_0=0e-7Watt/m^3), + Conductivity = T_Conductivity_Whittington(), # T-dependent k + HeatCapacity = T_HeatCapacity_Whittington(), # T-dependent cp + Melting = MeltingParam_Assimilation() # Quadratic parameterization as in Tierney et al. + ), + SetMaterialParams(Name="Intruded rocks", Phase=2, + Density = ConstantDensity(ρ=2700kg/m^3), # used in the parameterisation of Whittington + LatentHeat = ConstantLatentHeat(Q_L=2.67e5J/kg), + RadioactiveHeat = ExpDepthDependentRadioactiveHeat(H_0=0e-7Watt/m^3), + Conductivity = T_Conductivity_Whittington(), # T-dependent k + HeatCapacity = T_HeatCapacity_Whittington(), # T-dependent cp + Melting = SmoothMelting(MeltingParam_Quadratic(T_s=(700+273.15)K, T_l=(1100+273.15)K)) + ) + ) + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_params); \ No newline at end of file diff --git a/examples/MTK_GMG_2D_example2.jl b/examples/MTK_GMG_2D_example2.jl new file mode 100644 index 0000000..7d03655 --- /dev/null +++ b/examples/MTK_GMG_2D_example2.jl @@ -0,0 +1,192 @@ +# 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 + +const USE_GPU=false; +using MagmaThermoKinematics +if USE_GPU + environment!(:gpu, Float64, 2) # initialize parallel stencil in 2D +else + environment!(:cpu, Float64, 2) # initialize parallel stencil in 2D +end +using MagmaThermoKinematics.Diffusion2D +using MagmaThermoKinematics.MTK_GMG # Allow overwriting user routines +using Plots +using Random +using GeophysicalModelGenerator + +# Test setup +println("Example 2 of the MTK - GMG integration") + +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) +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) +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)) + +# Intersect with topography +Below = BelowSurface(Data_2D, Topo_cart) +Data_2D.fields.Phases[Below] .= 1 + +# Set Moho +@views Data_2D.fields.Phases[Data_2D.z.val .< -30.0] .= 2 + +# Set T: +gradient = 30 +Data_2D.fields.Temp .= -Data_2D.z.val*gradient +@views Data_2D.fields.Temp[Data_2D.fields.Temp .< 10.0] .= 10 + +# Set thermal anomaly +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 +@static if USE_GPU + function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + println("$(Num.it), Time=$(round(Num.time/Num.SecYear)) yrs; max(T) = $(round(maximum(Arrays.Tnew)))") + return nothing + end +else + function MTK_GMG.MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + if mod(Num.it,Num.CreateFig_steps)==0 + x_1d = Grid.coord1D[1]/1e3; + z_1d = Grid.coord1D[2]/1e3; + temp_data = Array(Arrays.Tnew)' + ϕ_data = Array(Arrays.ϕ)' + phase_data = Float64.(Array(Arrays.Phases))' + + # remove topo on plots + ind = findall(phase_data .== 0) + phase_data[ind] .= NaN + temp_data[ind] .= NaN + + t = Num.time/SecYear/1e3; + + p=plot(layout=grid(1,2) ) + + Plots.heatmap!(p[1],x_1d, z_1d, temp_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Temperature, t=$(round(t)) kyrs", aspect_ratio=:equal, ylimits=(minimum(z_1d),2)) + Plots.heatmap!(p[2],x_1d, z_1d, ϕ_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Melt fraction", clims=(0,1), aspect_ratio=:equal, ylimits=(minimum(z_1d),2)) + #Plots.heatmap!(p[2],x_1d, z_1d, phase_data, c=:viridis, xlabel="x [km]", ylabel="z [km]", title="Melt fraction", aspect_ratio=:equal, ylimits=(minimum(z_1d),2)) + + display(p) + end + return nothing + end +end + +function MTK_GMG.MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + return nothing +end + +function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + println("$(Num.it), Time=$(round(Num.time/Num.SecYear/1e3, digits=3)) kyrs; max(T) = $(round(maximum(Arrays.Tnew)))") + return nothing +end + +""" + MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + +Update time-dependent properties during a simulation +""" +function MTK_GMG.MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + push!(time_props.Time_vec, Num.time); # time + push!(time_props.MeltFraction, sum( Arrays.ϕ)/(Num.Nx*Num.Nz)); # melt fraction + + ind = findall(Arrays.T.>700); + if ~isempty(ind) + Tav_magma_Time = sum(Arrays.T[ind])/length(ind) # average T of part with magma + else + Tav_magma_Time = NaN; + end + push!(time_props.Tav_magma, Tav_magma_Time); # average magma T + push!(time_props.Tmax, maximum(Arrays.T)); # maximum magma T + return nothing +end + +# Define a new structure with time-dependent properties +@with_kw mutable struct TimeDepProps1 <: TimeDependentProperties + Time_vec::Vector{Float64} = []; # Center of dike + MeltFraction::Vector{Float64} = []; # Melt fraction over time + Tav_magma::Vector{Float64} = []; # Average magma + Tmax::Vector{Float64} = []; # Max magma temperature + Tmax_1::Vector{Float64} = []; # Another magma temperature vector +end + +# Define numerical parameters +Num = NumParam( SimName = "Unzen1", + dim = 2, + maxTime_Myrs = 0.005, + SaveOutput_steps = 25, + CreateFig_steps = 5, + USE_GPU = USE_GPU, + ω = 0.5, + AddRandomSills = true, + RandomSills_timestep= 5); + +# dike parameters +Dike_params = DikeParam(Type = "ElasticDike", + InjectionInterval_year = 1000, # flux= 14.9e-6 km3/km2/yr + W_in = 5e3, + H_in = 250, + H_ran = 5000, + W_ran = 5000, # width of random injection area + nTr_dike = 2000, + Dip_ran = 45, # angle aroun d which we randomly change the dip + DikePhase = 3, # phase of dike + SillsAbove = -10e3 # below this we have dikes; above sills + ) + +# Define parameters for the different phases +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), + 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), + 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)), + 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), + Melting = SmoothMelting(MeltingParam_4thOrder())) # Marxer & Ulmer melting + ) + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_params, CartData_input=Data_2D, time_props=TimeDepProps1()); # start the main code \ No newline at end of file diff --git a/examples/Plot_ZASSY_UCLA_vs_MTK_2D.jl b/examples/Plot_ZASSY_UCLA_vs_MTK_2D.jl new file mode 100644 index 0000000..88c9f6e --- /dev/null +++ b/examples/Plot_ZASSY_UCLA_vs_MTK_2D.jl @@ -0,0 +1,122 @@ +# This creates a plot of the MTK result vs the Geneva result run by Gregor Weber +# +# +using MAT +using CairoMakie +#using GLMakie +using DelimitedFiles + +#SimName_UCLA_final = "UCLA_data/isot1_01_000.dat" +#SimName_UCLA_final = "UCLA_data/isot1_03_075.dat" +SimName_UCLA_final = "UCLA_data/isot1_05_113.dat" +SimName_UCLA_inter = "UCLA_data/isot1_02_050.dat" + + +#SimName = "Zassy_UCLA_ellipticalIntrusion" # name of simulation +#SimName = "Zassy_UCLA_ellipticalIntrusion_2D" # name of simulation +#SimName = "Zassy_UCLA_ellipticalIntrusion_constant_k" # name of simulation +SimName = "Zassy_UCLA_ellipticalIntrusion_variable_k" +#SimName = "Zassy_UCLA_ellipticalIntrusion_constantkcp" # name of simulation +#SimName = "Zassy_UCLA_ellipticalIntrusion_variable_k_radioactiveheating" +SimName = "Zassy_UCLA_ellipticalIntrusion_constant_k_radioactiveheating" + +#it_final = 28800 +it_final = 43584 +it_inter = 19200 +#title_str = "isothermal lower BC, UCLA setup, k=3.35W/m/K" +#save_str = "Zassy_UCLA_ellipticalIntrusion_constant_k" + +#title_str = "flux lower BC with qm=76mW/m^2, UCLA setup, k=k(T), no radioactive heating" +#save_str = "Zassy_UCLA_ellipticalIntrusion_variable_k" + +#title_str = "UCLA setup, qm_bottom=76mW/m2, k=k(T), radioactive heating, H0=1e-6W/m3" +#save_str = "Zassy_UCLA_ellipticalIntrusion_variable_k_radioactiveheating" + +title_str = "UCLA setup, qm_bottom=134mW/m2, k=3.35W/K/m, radioactive heating: H0=1e-6W/m3, hr=10km" +save_str = "ZASSy_UCLA_ellipticalIntrusion_constant_k_radioactiveheating" + + +# Load Oscar's final results: +T_full = readdlm(SimName_UCLA_final) +T_o_final = T_full[:,301:end]'; +z_o_ref = Vector(0:-(20e3/(size(T_o_final,2)-1)):-20e3) +x_o_ref = Vector(0:(30e3/(size(T_o_final,1)-1)):30e3) + +T_full = readdlm(SimName_UCLA_inter) +T_o_inter = T_full[:,301:end]'; + +# Load MTK results +T_MTK = matopen("$(SimName)/$(SimName)_$(it_final).mat"); + Tnew_final = read(T_MTK, "Tnew") + #Tnew_final = read(T_MTK, "T") + x = read(T_MTK, "x") + z = read(T_MTK, "z") + time_final = read(T_MTK, "time") + dike_poly = read(T_MTK, "dike_poly") +close(T_MTK) + + +T_MTK = matopen("$(SimName)/$(SimName)_$it_inter.mat"); + Tnew_inter = read(T_MTK, "Tnew") + x = read(T_MTK, "x") + z = read(T_MTK, "z") + time_inter = read(T_MTK, "time") +close(T_MTK) + +T_lin = -(801.119-25)/20e3*z .+ 25 + + +fig = Figure(resolution = (2400,800)) + + +# 1D figure with cross-sections +SecYear = 3600*24*365.25 + +# 2D temperature plot UCLA +ax2=Axis(fig[1, 1],xlabel = "Width [km]", ylabel = "Depth [km]", title = "UCLA model") +Tcon = T_o_final; + +Tcon[Tcon.>= 999.9] .= 999.9 +co = contourf!(fig[1, 1], x_o_ref/1e3, z_o_ref/1e3, Tcon, levels = 0:50:1000,colormap = :jet) +if maximum(T_o_final)>691 + co1 = contour!(fig[1, 1], x_o_ref/1e3, z_o_ref/1e3, Tcon, levels = 690:691, color=:black) # solidus +end +limits!(ax2, 0, 30, -20, 0) + +# 2D temperature plot MTK +time_Myrs_rnd = round(time_final/SecYear/1e6,digits=3) +time_Myrs_int_rnd = round(time_inter/SecYear/1e6,digits=3) + +ax2=Axis(fig[1, 2],xlabel = "Width [km]", title = "MagmaThermoKinematics.jl, t=$time_Myrs_rnd Myrs") +co = contourf!(fig[1, 2], x/1e3, z/1e3, Tnew_final, levels = 0:50:1000,colormap = :jet) +if maximum(Tnew_final)>691 + co1 = contour!(fig[1, 2], x/1e3, z/1e3, Tnew_final, levels = 690:691, color=:black) # solidus +end +pl = lines!(fig[1, 2], dike_poly[1]/1e3, dike_poly[2]/1e3, color = :black, linestyle=:dot, linewidth=1.5) +limits!(ax2, 0, 30, -20, 0) +Colorbar(fig[1, 3], co, label = "Temperature [ᵒC]") + + + +ax3 = Axis(fig[1,4], xlabel = "Temperature [ᵒC]", ylabel = "Depth [km]", title = "$title_str") +#lines!(fig[1,4],T_lin,z/1e3,label="t=0", color=:black) + +lines!(fig[1,4],Tnew_inter[1,:], z/1e3,label="Center, MTK, t=$time_Myrs_int_rnd", color=:orange) +lines!(fig[1,4],T_o_inter[1,:], z_o_ref/1e3,label="Center, UCLA",linestyle=:dot, color=:orange) +lines!(fig[1,4],Tnew_inter[end,:],z/1e3,label="Side, MTK, t=$time_Myrs_int_rnd", color=:green) +lines!(fig[1,4],T_o_inter[end,:],z_o_ref/1e3,label="Side, UCLA",linestyle=:dot, color=:green) + +lines!(fig[1,4],Tnew_final[1,:], z/1e3,label="Center, MTK, t=$time_Myrs_rnd", color=:red) +lines!(fig[1,4],T_o_final[1,:], z_o_ref/1e3,label="Center, UCLA",linestyle=:dot, color=:red) +lines!(fig[1,4],Tnew_final[end,:],z/1e3,label="Side, MTK, t=$time_Myrs_rnd", color=:blue) +lines!(fig[1,4],T_o_final[end,:],z_o_ref/1e3,label="Side, UCLA",linestyle=:dot, color=:blue) + + + +axislegend(ax3, position = :lb) +limits!(ax3, 0, 1100, -20, 0) + +#save("Comparison_Geneva_MTK_zero_flux_$(time_Myrs_rnd).png", fig) +#save("Comparison_Geneva_MTK_isothermal_$(time_Myrs_rnd).png", fig) + +save("$(save_str)_$(time_Myrs_rnd)Myrs.png", fig) diff --git a/examples/Plot_ZirconAgeStatistics.jl b/examples/Plot_ZirconAgeStatistics.jl new file mode 100644 index 0000000..91b5b01 --- /dev/null +++ b/examples/Plot_ZirconAgeStatistics.jl @@ -0,0 +1,36 @@ +using GLMakie +using JLD2 + +# Specify directory names +dirnames = ["Zassy_UCLA_ellipticalIntrusion_constant_k_radioactiveheating", + "Zassy_Geneva_zeroFlux_variable_k_4thordermelt", + "Zassy_Geneva_zeroFlux_variable_k_CaricchiMelting"] + +col = [:red,:blue,:green] +leg = ["central injection","sill underaccretion","sill underaccretion Caricchi Melting"] + +# Create figure +fig = Figure(resolution = (2000,1000)) +ax1=Axis(fig[1, 1],xlabel = "Age [Ma]", ylabel = "T_{average} magma") +xlims!(ax1, 0, 1.5) +ylims!(ax1, 400, 1000) + +ax2=Axis(fig[1, 2],xlabel = "Age [Ma]", ylabel = "Zircon age cummulative probability %") +limits!(ax2, 0, 1.5, 0,1) + +for i = 1:length(dirnames) + + # Load file + filename = dirnames[i]*"/ZirconAges.jld2"; + Age_Ma,T_av_time,cum_PDF,norm_PDF = JLD2.load(filename, "Age_Ma","T_av_time","cum_PDF","norm_PDF"); + + # Plot T + lines!(ax1,Age_Ma, T_av_time, color=col[i]) + + # Plot PDF + lines!(ax2,Age_Ma, cum_PDF, color=col[i], label=leg[i]) + +end +axislegend(position=:rb) + +fig diff --git a/examples/Topo_cart.jld2 b/examples/Topo_cart.jld2 new file mode 100644 index 0000000..49af93e Binary files /dev/null and b/examples/Topo_cart.jld2 differ diff --git a/src/Diffusion.jl b/src/Diffusion.jl index 734b8a3..f57aa82 100644 --- a/src/Diffusion.jl +++ b/src/Diffusion.jl @@ -290,10 +290,10 @@ end In-place function that creates 3D arrays with `x`,`y` and `z` coordinates using the `Grid` information """ -function GridArray!(X::AbstractArray, Z::AbstractArray, Grid::GridData) +function GridArray!(X::AbstractArray, Y::AbstractArray, Z::AbstractArray, Grid::GridData) @parallel (1:Grid.N[1], 1:Grid.N[2], 1:Grid.N[3]) GridArray!(X, Y, Z, Grid.coord1D[1], Grid.coord1D[2], Grid.coord1D[3]) - return nothing + return end @parallel_indices (i,j,k) function GridArray!(X::AbstractArray, Y::AbstractArray, Z::AbstractArray, x::AbstractVector, y::AbstractVector,z::AbstractVector) diff --git a/src/Dikes.jl b/src/Dikes.jl index 014e45a..55d62d4 100644 --- a/src/Dikes.jl +++ b/src/Dikes.jl @@ -10,40 +10,36 @@ temperature field Structure that holds the geometrical parameters of the dike, which are slightly different depending on whether we consider a 2D or a 3D case - General form: - Dike(Width=.., Thickness=.., Center=[], Angle=[], Type="..", T=.., ΔP=.., E=.., ν=.., Q=..) - - with: - - [Width]: width of dike (optional, will be computed automatically if ΔP and Q are specified) - - [Thickness]: (maximum) thickness of dike (optional, will be computed automatically if ΔP and Q are specified) +General form: + Dike(Width=.., Thickness=.., Center=[], Angle=[], Type="..", T=.., ΔP=.., E=.., ν=.., Q=..) + +with +- [Width]: width of dike (optional, will be computed automatically if ΔP and Q are specified +- [Thickness]: (maximum) thickness of dike (optional, will be computed automatically if ΔP and Q are specified) +- Center: center of the dike + 2D - [x; z] + 3D - [x; y; z] - Center: center of the dike - 2D - [x; z] - 3D - [x; y; z] - - Angle: Dip (and strike) angle of dike - 2D - [Dip] - 3D - [Strike; Dip] - - Type: Type of dike - "SquareDike" - square dike area - "SquareDike_TopAccretion" - square dike area, which grows by underaccreting - "CylindricalDike_TopAccretion" - cylindrical dike area, which grows by underaccreting - "CylindricalDike_TopAccretion_FullModelAdvection" - cylindrical dike area, which grows by underaccreting; also material to the side of the dike is moved downwards - "ElasticDike" - penny-shaped elastic dike in elastic halfspace - "EllipticalIntrusion" - elliptical dike intrusion area with radius Width/2 and height Height/2 - - T: Temperature of the dike [Celcius] - - ν: Poison ratio of host rocks - - E: Youngs modulus of host rocks [Pa] - - [ΔP]: Overpressure of dike w.r.t. host rock [Pa], (optional in case we want to compute width/length directly) - - [Q]: Volume of magma within dike [m^3], + Angle: Dip (and strike) angle of dike + 2D - [Dip] + 3D - [Strike; Dip] + + Type: Type of dike + "SquareDike" - square dike area + "SquareDike_TopAccretion" - square dike area, which grows by underaccreting + "CylindricalDike_TopAccretion" - cylindrical dike area, which grows by underaccreting + "CylindricalDike_TopAccretion_FullModelAdvection" - cylindrical dike area, which grows by underaccreting; also material to the side of the dike is moved downwards + "ElasticDike" - penny-shaped elastic dike in elastic halfspace + "EllipticalIntrusion" - elliptical dike intrusion area with radius Width/2 and height Height/2 + + T: Temperature of the dike [Celcius] + + ν: Poison ratio of host rocks + + E: Youngs modulus of host rocks [Pa] + + [ΔP]: Overpressure of dike w.r.t. host rock [Pa], (optional in case we want to compute width/length directly + [Q]: Volume of magma within dike [m^3], All parameters can be specified through keywords as shown above. @@ -672,8 +668,10 @@ function AddDike(Tfield,Tr, Grid,dike, nTr_dike) new_tracer = Tracer(num=number, coord=pt_new, T=T, Phase=PhaseDike); # Create new tracer if !isassigned(Tr,1) - StructArrays.foreachfield(v -> deleteat!(v, 1), Tr) # Delete first (undefined) row of tracer StructArray. Assumes that Tr is defined as Tr = StructArray{Tracer}(undef, 1) + if length(Tr)==0 + StructArrays.foreachfield(v -> deleteat!(v, 1), Tr) # Delete first (undefined) row of tracer StructArray. Assumes that Tr is defined as Tr = StructArray{Tracer}(undef, 1) + end Tr = StructArray([new_tracer]); # Create tracer array else push!(Tr, new_tracer); # Add new point to existing array diff --git a/src/Grid.jl b/src/Grid.jl index b39d5d8..76eda30 100644 --- a/src/Grid.jl +++ b/src/Grid.jl @@ -1,6 +1,7 @@ module Grid # This creates a computational grid import Base: show +import GeophysicalModelGenerator: CartData export GridData, CreateGrid @@ -22,7 +23,7 @@ Creates a 1D, 2D or 3D grid of given size. Grid can be created by defining the s Spacing is assumed to be constant -Note: since this is mostly for Solid Earth geoscience applications, the second dimension is called z (vertical) +Note: since this is mostly for Solid Earth geoscience applications, in 2D the second dimension is called z (vertical) # Examples ==== @@ -55,11 +56,12 @@ function CreateGrid(; if !isnothing(extent) x,y,z = nothing, nothing, nothing x = (0., extent[1]) - if dim>1 + if dim==2 z = (-extent[2], 0.0) # vertical direction (negative) end - if dim>2 - y = (0., extent[3]) + if dim==3 + y = (0., extent[2]) + z = (-extent[3], 0.0) # vertical direction (negative) end end @@ -71,8 +73,8 @@ function CreateGrid(; L = (x[2] - x[1], z[2] - z[1]) X₁= (x[1], z[1]) else - L = (x[2] - x[1], z[2] - z[1], y[2] - y[1]) - X₁= (x[1], z[1], y[1]) + L = (x[2] - x[1], y[2] - y[1], z[2] - z[1]) + X₁= (x[1], y[1], z[1]) end Xₙ = X₁ .+ L Δ = L ./ (N .- 1) @@ -83,7 +85,7 @@ function CreateGrid(; coord1D = (coord1D..., range(X₁[idim], Xₙ[idim]; length = N[idim] )) end - # Generate 1D coordinate arrays centers in all directionbs + # Generate 1D coordinate arrays centers in all directions coord1D_cen=() for idim=1:dim coord1D_cen = (coord1D_cen..., range(X₁[idim]+Δ[idim]/2, Xₙ[idim]-Δ[idim]/2; length = N[idim]-1 )) @@ -94,6 +96,30 @@ function CreateGrid(; end + +""" + Grid = CreateGrid(d::CartData; m_to_km::Bool=true) +Creates a (regularly spaced) grid that can be used within MTK, from a 2D or 3D `CartData` object (generated within the GeophysicalModelGenerator) +""" +function CreateGrid(d::CartData; m_to_km::Bool=true) + Nx,Ny,Nz = size(d.x.val) + if Nz==1; dim=2; else dim=3; end + if m_to_km; scaling = 1e3; else scaling = 1; end + + x = extrema(d.x.val).*scaling; + y = extrema(d.y.val).*scaling; + z = extrema(d.z.val).*scaling; + + if dim == 3 + Grid = CreateGrid(size=(Nx,Ny,Nz), x=x, y=y, z=z) + else + Grid = CreateGrid(size=(Nx,Ny), x=x, z=z) + end + return Grid +end + + + # view grid object function show(io::IO, g::GridData{FT, DIM}) where {FT, DIM} diff --git a/src/MTK_GMG.jl b/src/MTK_GMG.jl new file mode 100644 index 0000000..fe41456 --- /dev/null +++ b/src/MTK_GMG.jl @@ -0,0 +1,424 @@ +# various routines that are shared between the 2D and 3D MTK_GMG routines +""" + MTK_GMG +This contains various user callback routines that are shared between the 2D and 3D MTK_GMG routines. +You can overwrite this in your own code to customize the simulation. + +""" +module MTK_GMG + +using Parameters +using GeoParams +using GeophysicalModelGenerator +using StructArrays +using MagmaThermoKinematics.Grid +using MagmaThermoKinematics.Data +import MagmaThermoKinematics: NumericalParameters, DikeParameters, TimeDependentProperties +import MagmaThermoKinematics: update_Tvec!, Dike, InjectDike, km³, kyr, Myr, CreateDikePolygon +import MagmaThermoKinematics: PhasesFromTracers! +SecYear = 3600*24*365.25; + +using CUDA + +""" + Analytical geotherm used for the UCLA setups, which includes radioactive heating +""" +function AnalyticalGeotherm!(T, Z, Tsurf, qm, qs, k, hr) + T .= @. Tsurf - (qm/k)*Z + (qs-qm)*hr/k*( 1.0 - exp(Z/hr)) + return nothing +end + +""" + Tracers = MTK_inject_dikes(Grid, Num, Arrays, Mat_tup, Dikes, Tracers, Tnew_cpu) + +Function that injects dikes once in a while +""" +function MTK_inject_dikes(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters, Tracers::StructVector, Tnew_cpu) + + if floor(Num.time/Dikes.InjectionInterval)> Dikes.dike_inj + Dikes.dike_inj = floor(Num.time/Dikes.InjectionInterval) # Keeps track on what was injected already + if Num.dim==2 + T_bottom = Tnew_cpu[:,1] + else + T_bottom = Tnew_cpu[:,:,1] + end + dike = Dike(W=Dikes.W_in, H=Dikes.H_in, Type=Dikes.Type, T=Dikes.T_in_Celsius, Center=Dikes.Center[:], Angle=Dikes.Angle, Phase=Dikes.DikePhase); # "Reference" dike with given thickness,radius and T + Tnew_cpu .= Array(Arrays.T) + + Tracers, Tnew_cpu,Vol,Dikes.dike_poly, VEL = InjectDike(Tracers, Tnew_cpu, Grid.coord1D, dike, Dikes.nTr_dike, dike_poly=Dikes.dike_poly); # Add dike, move hostrocks + + if Num.flux_bottom_BC==false + # Keep bottom T constant (advection modifies this) + if Num.dim==2 + Tnew_cpu[:,1] .= T_bottom + else + Tnew_cpu[:,:,1] .= T_bottom + end + end + + Arrays.T .= Data.Array(Tnew_cpu) + Dikes.InjectVol += Vol # Keep track of injected volume + Qrate = Dikes.InjectVol/Num.time + Dikes.Qrate_km3_yr = Qrate*SecYear/km³ + Qrate_km3_yr_km2 = Dikes.Qrate_km3_yr/(pi*(Dikes.W_in/2/1e3)^2) + println(" Added new dike; time=$(Num.time/kyr) kyrs, total injected magma volume = $(Dikes.InjectVol/km³) km³; rate Q= $(Dikes.Qrate_km3_yr) km³yr⁻¹") + + if Num.advect_polygon==true && isempty(Dikes.dike_poly) + Dikes.dike_poly = CreateDikePolygon(dike); # create dike for the 1th time + end + + if length(Mat_tup)>1 + PhasesFromTracers!(Arrays.Phases, Grid, Tracers, BackgroundPhase=Dikes.BackgroundPhase, InterpolationMethod="Constant"); # update phases from grid + + # Ensure that we keep the initial phase of the area (host rocks are not deformable) + if Num.keep_init_RockPhases==true + for i in eachindex(Arrays.Phases) + if Arrays.Phases[i] != Dikes.DikePhase + Arrays.Phases[i] = Arrays.Phases_init[i] + end + end + end + end + + end + + return Tracers +end + +""" + MTK_display_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + +Function that creates plots +""" +function MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + + return nothing +end + +""" + MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + +Function that prints output to the REPL +""" +function MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + + return nothing +end + +""" + MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + +Update time-dependent properties during a simulation +""" +function MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + push!(time_props.Time_vec, Num.time); # time + push!(time_props.MeltFraction, sum( Arrays.ϕ)/(Num.Nx*Num.Nz)); # melt fraction + + ind = findall(Arrays.T.>700); + if ~isempty(ind) + Tav_magma_Time = sum(Arrays.T[ind])/length(ind) # average T of part with magma + else + Tav_magma_Time = NaN; + end + push!(time_props.Tav_magma, Tav_magma_Time); # average magma T + push!(time_props.Tmax, maximum(Arrays.T)); # maximum magma T + + return nothing +end + +""" + MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters) + +Initialize temperature and phases +""" +function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters) + # Initalize T + Arrays.T_init .= @. Num.Tsurface_Celcius - Arrays.Z*Num.Geotherm; # Initial (linear) temperature profile + + # Open pvd file if requested + if Num.Output_VTK + name = joinpath(Num.SimName,Num.SimName*".pvd") + Num.pvd = Movie_Paraview(name=name, Initialize=true); + end + + return nothing +end + +""" + MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::CartData) + +Initialize temperature and phases +""" +function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::Union{Nothing,CartData}) + # Initalize T from CartData set + # NOTE: this almost certainly requires changes if we use GPUs + + if Num.USE_GPU + if Num.dim==2 + Arrays.T_init .= Data.Array(CartData_input.fields.Temp[:,:,1]) + Arrays.Phases .= Data.Array(CartData_input.fields.Phases[:,:,1]); + Arrays.Phases_init .= Data.Array(CartData_input.fields.Phases[:,:,1]); + else + Arrays.T_init .= Data.Array(CartData_input.fields.Temp) + Arrays.Phases .= Data.Array(CartData_input.fields.Phases); + Arrays.Phases_init .= Data.Array(CartData_input.fields.Phases); + end + else + if Num.dim==2 + Arrays.T_init .= CartData_input.fields.Temp[:,:,1]; + Arrays.Phases .= CartData_input.fields.Phases[:,:,1]; + Arrays.Phases_init .= CartData_input.fields.Phases[:,:,1]; + else + Arrays.T_init .= CartData_input.fields.Temp; + Arrays.Phases .= CartData_input.fields.Phases; + Arrays.Phases_init .= CartData_input.fields.Phases; + end + end + + # open pvd file if requested + if Num.Output_VTK + name = joinpath(Num.SimName,Num.SimName*".pvd") + Num.pvd = Movie_Paraview(name=name, Initialize=true); + end + + return nothing +end + + +""" + MTK_finalize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::CartData) + +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) + end + + return nothing +end + + +""" + MTK_update_Arrays!(Arrays::NamedTuple, Grid::GridData, Dikes::DikeParameters, Num::NumericalParameters) + +Update arrays and structs of the simulation (in case you want to change them during a simulation) +You can use this, for example, to change the size and location of an intruded dike +""" +function MTK_update_ArraysStructs!(Arrays::NamedTuple, Grid::GridData, Dikes::DikeParameters, Num::NumericalParameters) + + if Num.AddRandomSills && mod(Num.it,Num.RandomSills_timestep)==0 + # This randomly changes the location and orientation of the sills + if Num.dim==2 + Loc = [Dikes.W_ran; Dikes.H_ran] + else + Loc = [Dikes.W_ran; Dikes.L_ran; Dikes.H_ran] + end + + # Randomly change location of center of dike/sill + cen = (Grid.max .+ Grid.min)./2 .+ rand(-0.5:1e-3:0.5, Num.dim).*Loc; + + Dip = rand(-Dikes.Dip_ran/2.0 : 0.1: Dikes.Dip_ran/2.0) + Strike = rand(-Dikes.Strike_ran/2.0 : 0.1: Dikes.Strike_ran/2.0) + + if cen[end]κ_max + κ_max = κ + end + end + Num.κ_time = κ_max; + Num.Δ = [dx, dz] + Num.Δmin = minimum(Num.Δ[Num.Δ.>0]); # minimum grid spacing + + Num.dt = Num.fac_dt*(Num.Δmin^2)./Num.κ_time/4; # timestep + + Num.dx = dx; + Num.dz = dz; + + Num.nt = floor(Num.maxTime/Num.dt) + + return Num +end + +function Setup_Model_CartData_3D(d::CartData, Num::NumericalParameters, Mat_tup::Tuple) + x = extrema(d.x.val.*1e3) + y = extrema(d.y.val.*1e3) + z = extrema(d.z.val.*1e3) + + Num.W = (x[2]-x[1]) + Num.L = (y[2]-y[1]) + Num.H = (z[2]-z[1]) + Num.Nx = size(d.x)[1] + Num.Ny = size(d.x)[2] + Num.Nz = size(d.x)[3] + + dx = (x[2]-x[1])/(Num.Nx-1) + dy = (y[2]-y[1])/(Num.Ny-1) + dz = (z[2]-z[1])/(Num.Nz-1) + + # estimate maximum thermal diffusivity from Mat_tup + κ_max = Num.κ_time + for mm in Mat_tup + if hasfield(typeof(mm.Conductivity[1]),:k) + k = NumValue(mm.Conductivity[1].k) + else + k = 3; + end + if hasfield(typeof(mm.HeatCapacity[1]),:cp) + cp = NumValue(mm.HeatCapacity[1].cp) + else + cp = 1050; + end + if hasfield(typeof(mm.Density[1]),:ρ) + ρ = NumValue(mm.Density[1].ρ) + else + ρ = 2700; + end + κ = k/(cp*ρ) + if κ>κ_max + κ_max = κ + end + end + Num.κ_time = κ_max; + Num.Δ = [dx, dy, dz] + Num.Δmin = minimum(Num.Δ[Num.Δ.>0]); # minimum grid spacing + + Num.dt = Num.fac_dt*(Num.Δmin^2)./Num.κ_time/4; # timestep + Num.dx = dx; + Num.dy = dy; + Num.dz = dz; + + Num.nt = floor(Num.maxTime/Num.dt) + + return Num +end + + + +end \ No newline at end of file diff --git a/src/MTK_GMG_2D.jl b/src/MTK_GMG_2D.jl new file mode 100644 index 0000000..416712c --- /dev/null +++ b/src/MTK_GMG_2D.jl @@ -0,0 +1,199 @@ +# This contains the 2D routine to create MTK simulations (using GeoParams) + +# +module MTK_GMG_2D +using ParallelStencil +using ParallelStencil.FiniteDifferences2D +using Parameters +using StructArrays +using GeophysicalModelGenerator +using CUDA + +using MagmaThermoKinematics.Diffusion2D +using MagmaThermoKinematics.MTK_GMG +using MagmaThermoKinematics + +const SecYear = 3600*24*365.25; + +export MTK_GeoParams_2D + +#----------------------------------------------------------------------------------------- +""" + Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(Mat_tup::Tuple, Num::NumericalParameters, Dikes::DikeParameters; CartData_input=nothing, time_props::TimeDependentProperties = TimeDepProps()); + +Main routine that performs a 2D or 2D axisymmetric thermal diffusion simulation with injection of dikes. + +Parameters +==== +- `Mat_tup::Tuple`: Tuple of material properties. +- `Num::NumericalParameters`: Numerical parameters. +- `Dikes::DikeParameters`: Dike parameters. +- `CartData_input::CartData`: Optional input of a CartData structure generated with GeophysicalModelGenerator. +- `time_props::TimeDependentProperties`: Optional input of a `TimeDependentProperties` structure. + +Customizable functions +==== +There are a few functions that you can overwrite in your user code to customize the simulation: + +- `MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters)` +- `MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters)` +- `MTK_update_ArraysStructs!(Arrays::NamedTuple, Grid::GridData, Dikes::DikeParameters, Num::NumericalParameters)` +- `MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input)` +- `MTK_updateTracers(Grid::GridData, Arrays::NamedTuple, Tracers::StructArray, Dikes::DikeParameters, time_props::TimeDependentProperties, Num::NumericalParameters)` +- `MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArray, Dikes::DikeParameters, time_props::TimeDependentProperties, Num::NumericalParameters, CartData_input::CartData)` +- `MTK_inject_dikes(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters, Tracers::StructVector, Tnew_cpu)` +- `MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters)` +- `MTK_finalize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::CartData)` + +""" +@views function MTK_GeoParams_2D(Mat_tup::Tuple, Num::NumericalParameters, Dikes::DikeParameters; CartData_input::Union{Nothing,CartData}=nothing, time_props::TimeDependentProperties = TimeDepProps()); + + # Change parameters based on CartData input + Num.dim = 2; + 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))") + end + + Num = MTK_GMG.Setup_Model_CartData(CartData_input, Num, Mat_tup) + end + + + # Array & grid initializations --------------- + Arrays = CreateArrays(Dict( (Num.Nx, Num.Nz )=>(T=0,T_K=0, Tnew=0, T_init=0, T_it_old=0, Kc=1, Rho=1, Cp=1, Hr=0, Hl=0, ϕ=0, dϕdT=0,dϕdT_o=0, R=0, Z=0, P=0), + (Num.Nx-1,Num.Nz )=>(qx=0,Kx=0, Rc=0), + (Num.Nx ,Num.Nz-1)=>(qz=0,Kz=0 ) + )) + + # Set up model geometry & initial T structure + if isnothing(CartData_input) + Grid = CreateGrid(size=(Num.Nx,Num.Nz), extent=(Num.W, Num.H)) + else + Grid = CreateGrid(CartData_input) + end + GridArray!(Arrays.R, Arrays.Z, Grid) + Arrays.Rc .= (Arrays.R[2:end,:] + Arrays.R[1:end-1,:])/2 # center points in x + # -------------------------------------------- + + Tracers = StructArray{Tracer}(undef, 1) # Initialize tracers + + # Update buffer & phases arrays -------------- + if Num.USE_GPU + # CPU buffers for advection + Tnew_cpu = Matrix{Float64}(undef, Num.Nx, Num.Nz) + Phi_melt_cpu = similar(Tnew_cpu) + Phases = CUDA.ones(Int64,Num.Nx,Num.Nz) + Phases_init = CUDA.ones(Int64,Num.Nx,Num.Nz) + else + Tnew_cpu = similar(Arrays.T) + Phi_melt_cpu = similar(Arrays.ϕ) + Phases = ones(Int64,Num.Nx,Num.Nz) + Phases_init = ones(Int64,Num.Nx,Num.Nz) + end + Arrays = (Arrays..., Phases=Phases, Phases_init=Phases_init); + + # Initialize Geotherm and Phases ------------- + if isnothing(CartData_input) + MTK_GMG.MTK_initialize!(Arrays, Grid, Num, Tracers, Dikes); + else + MTK_GMG.MTK_initialize!(Arrays, Grid, Num, Tracers, Dikes, CartData_input); + end + # -------------------------------------------- + + # check errors + unique_Phases = unique(Array(Arrays.Phases)); + phase_specified = [] + for mm in Mat_tup + push!(phase_specified, mm.Phase) + end + for u in unique_Phases + if !(u in phase_specified) + error("Properties for Phase $u are not specified in Mat_tup. Please add that") + end + end + + if any(isnan.(Arrays.T)) + error("NaNs in T; something is wrong") + end + + # Optionally set initial sill in models ------ + if Dikes.Type == "CylindricalDike_TopAccretion" + ind = findall( (Arrays.R.<=Dikes.W_in/2) .& (abs.(Arrays.Z.-Dikes.Center[2]) .< Dikes.H_in/2) ); + Arrays.T_init[ind] .= Dikes.T_in_Celsius; + if Num.advect_polygon==true + dike = Dike(W=Dikes.W_in,H=Dikes.H_in,Type=Dikes.Type,T=Dikes.T_in_Celsius, Center=Dikes.Center[:], Angle=Dikes.Angle, Phase=Dikes.DikePhase); # "Reference" dike with given thickness,radius and T + Dikes.dike_poly = CreateDikePolygon(dike); + end + end + # -------------------------------------------- + + # Initialise arrays -------------------------- + @parallel assign!(Arrays.Tnew, Arrays.T_init) + @parallel assign!(Arrays.T, Arrays.T_init) + + if isdir(Num.SimName)==false + mkdir(Num.SimName) # create simulation directory if needed + end; + # -------------------------------------------- + + for Num.it = 1:Num.nt # Time loop + Num.time += Num.dt; # Keep track of evolved time + + # Add new dike every X years ----------------- + Tracers = MTK_GMG.MTK_inject_dikes(Grid, Num, Arrays, Mat_tup, Dikes, Tracers, Tnew_cpu) + # -------------------------------------------- + + # Do a diffusion step, while taking T-dependencies into account + Nonlinear_Diffusion_step_2D!(Arrays, Mat_tup, Phases, Grid, Num.dt, Num) + # -------------------------------------------- + + # Update variables --------------------------- + # copy to cpu + Tnew_cpu .= Array(Arrays.Tnew) + Phi_melt_cpu .= Array(Arrays.ϕ) + + UpdateTracers_T_ϕ!(Tracers, Grid.coord1D, Tnew_cpu, Phi_melt_cpu); # Update info on tracers + + # copy back to gpu + Arrays.Tnew .= Data.Array(Tnew_cpu) + Arrays.ϕ .= Data.Array(Phi_melt_cpu) + + @parallel assign!(Arrays.T, Arrays.Tnew) + @parallel assign!(Arrays.Tnew, Arrays.T) + # -------------------------------------------- + + # Update info on tracers --------------------- + Tracers = MTK_GMG.MTK_updateTracers(Grid, Arrays, Tracers, Dikes, time_props, Num); + # -------------------------------------------- + + # Update time-dependent properties ----------- + MTK_GMG.MTK_update_TimeDepProps!(time_props, Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + # Visualize results -------------------------- + MTK_GMG.MTK_visualize_output(Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + # Save output to disk once in a while -------- + MTK_GMG.MTK_save_output(Grid, Arrays, Tracers, Dikes, time_props, Num, CartData_input); + # -------------------------------------------- + + # Optionally update arrays and structs (such as T or Dike) ------- + MTK_GMG.MTK_update_ArraysStructs!(Arrays, Grid, Dikes, Num) + # -------------------------------------------- + + # Display output ----------------------------- + MTK_GMG.MTK_print_output(Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + end + + # Finalize simulation ------------------------ + MTK_GMG.MTK_finalize!(Arrays, Grid, Num, Tracers, Dikes, CartData_input); + # -------------------------------------------- + + return Grid, Arrays, Tracers, Dikes, time_props +end # end of main function + +end \ No newline at end of file diff --git a/src/MTK_GMG_3D.jl b/src/MTK_GMG_3D.jl new file mode 100644 index 0000000..ad32ebd --- /dev/null +++ b/src/MTK_GMG_3D.jl @@ -0,0 +1,199 @@ +# This contains the 3D routine to create MTK simulations (using GeoParams) + +# +module MTK_GMG_3D + +using ParallelStencil +using ParallelStencil.FiniteDifferences3D +using Parameters +using StructArrays +using GeophysicalModelGenerator +using CUDA + +using MagmaThermoKinematics.Diffusion3D +using MagmaThermoKinematics.MTK_GMG +using MagmaThermoKinematics + + +const SecYear = 3600*24*365.25; + +export MTK_GeoParams_3D + + +#----------------------------------------------------------------------------------------- +""" + Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_3D(Mat_tup::Tuple, Num::NumericalParameters, Dikes::DikeParameters; CartData_input=nothing, time_props::TimeDependentProperties = TimeDepProps()); + +Main routine that performs a 3D thermal diffusion simulation with injection of dikes. + +Parameters +==== +- `Mat_tup::Tuple`: Tuple of material properties. +- `Num::NumericalParameters`: Numerical parameters. +- `Dikes::DikeParameters`: Dike parameters. +- `CartData_input::CartData`: Optional input of a CartData structure generated with GeophysicalModelGenerator. +- `time_props::TimeDependentProperties`: Optional input of a `TimeDependentProperties` structure. + +Customizable functions +==== +There are a few functions that you can overwrite in your user code to customize the simulation: + +- `MTK_visualize_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters)` +- `MTK_update_TimeDepProps!(time_props::TimeDependentProperties, Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters)` +- `MTK_update_ArraysStructs!(Arrays::NamedTuple, Grid::GridData, Dikes::DikeParameters, Num::NumericalParameters)` +- `MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input)` +- `MTK_updateTracers(Grid::GridData, Arrays::NamedTuple, Tracers::StructArray, Dikes::DikeParameters, time_props::TimeDependentProperties, Num::NumericalParameters)` +- `MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArray, Dikes::DikeParameters, time_props::TimeDependentProperties, Num::NumericalParameters, CartData_input::CartData)` +- `MTK_inject_dikes(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters, Tracers::StructVector, Tnew_cpu)` +- `MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters)` +- `MTK_finalize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::CartData)` + +""" +@views function MTK_GeoParams_3D(Mat_tup::Tuple, Num::NumericalParameters, Dikes::DikeParameters; CartData_input::Union{Nothing,CartData}=nothing, time_props::TimeDependentProperties = TimeDepProps()); + + # Change parameters based on CartData input + if !isnothing(CartData_input) + Num = MTK_GMG.Setup_Model_CartData(CartData_input, Num, Mat_tup) + end + + # Array & grid initializations --------------- + if Num.dim==2 + error("This is the 3D routine - use MTK_GeoParams_2D instead") + else + Arrays = CreateArrays(Dict( (Num.Nx, Num.Ny , Num.Nz )=>(T=0,T_K=0, Tnew=0, T_init=0, T_it_old=0, Kc=1, Rho=1, Cp=1, Hr=0, Hl=0, ϕ=0, dϕdT=0,dϕdT_o=0, R=0, X=0, Y=0, Z=0, P=0), + (Num.Nx-1,Num.Ny , Num.Nz )=>(qx=0,Kx=0), + (Num.Nx ,Num.Ny-1, Num.Nz )=>(qy=0,Ky=0), + (Num.Nx ,Num.Ny , Num.Nz-1)=>(qz=0,Kz=0 ) + )) + end + + # Set up model geometry & initial T structure + if isnothing(CartData_input) + Grid = CreateGrid(size=(Num.Nx,Num.Ny,Num.Nz), x = (-Num.W/2, Num.W/2), y = (-Num.L/2, Num.L/2), z=(-Num.H, 0.0)) + else + Grid = CreateGrid(CartData_input) + end + GridArray!(Arrays.X, Arrays.Y, Arrays.Z, Grid) + # -------------------------------------------- + + Tracers = StructArray{Tracer}(undef, 1) # Initialize tracers + + # Update buffer & phases arrays -------------- + if Num.USE_GPU + # CPU buffers for advection + Tnew_cpu = zeros(Float64, Num.Nx, Num.Ny, Num.Nz) + Phi_melt_cpu = similar(Tnew_cpu) + Phases = CUDA.ones(Int64,Num.Nx,Num.Ny,Num.Nz) + Phases_init = CUDA.ones(Int64,Num.Nx,Num.Ny,Num.Nz) + else + Tnew_cpu = similar(Arrays.T) + Phi_melt_cpu = similar(Arrays.ϕ) + Phases = ones(Int64,Num.Nx,Num.Ny,Num.Nz) + Phases_init = ones(Int64,Num.Nx,Num.Ny,Num.Nz) + end + Arrays = (Arrays..., Phases=Phases, Phases_init=Phases_init); + + # Initialize Geotherm and Phases ------------- + if isnothing(CartData_input) + MTK_GMG.MTK_initialize!(Arrays, Grid, Num, Tracers, Dikes); + else + MTK_GMG.MTK_initialize!(Arrays, Grid, Num, Tracers, Dikes, CartData_input); + end + + # check errors + unique_Phases = unique(Array(Arrays.Phases)); + phase_specified = [] + for mm in Mat_tup + push!(phase_specified, mm.Phase) + end + for u in unique_Phases + if !(u in phase_specified) + error("Properties for Phase $u are not specified in Mat_tup. Please add that") + end + end + + if any(isnan.(Arrays.T)) + error("NaNs in T; something is wrong") + end + # -------------------------------------------- + + # Optionally set initial sill in models ------ + if Dikes.Type == "CylindricalDike_TopAccretion" + ind = findall( (Arrays.R.<=Dikes.W_in/2) .& (abs.(Arrays.Z.-Dikes.Center[2]) .< Dikes.H_in/2) ); + Arrays.T_init[ind] .= Dikes.T_in_Celsius; + if Num.advect_polygon==true + dike = Dike(W=Dikes.W_in,H=Dikes.H_in,Type=Dikes.Type,T=Dikes.T_in_Celsius, Center=Dikes.Center[:], Angle=Dikes.Angle, Phase=Dikes.DikePhase); # "Reference" dike with given thickness,radius and T + Dikes.dike_poly = CreateDikePolygon(dike); + end + end + # -------------------------------------------- + + # Initialise arrays -------------------------- + @parallel assign!(Arrays.Tnew, Arrays.T_init) + @parallel assign!(Arrays.T, Arrays.T_init) + + if isdir(Num.SimName)==false mkdir(Num.SimName) end; # create simulation directory if needed + # -------------------------------------------- + + for Num.it = 1:Num.nt # Time loop + Num.time += Num.dt; # Keep track of evolved time + + # Add new dike every X years ----------------- + Tracers = MTK_GMG.MTK_inject_dikes(Grid, Num, Arrays, Mat_tup, Dikes, Tracers, Tnew_cpu) + # -------------------------------------------- + + # Do a diffusion step, while taking T-dependencies into account + Nonlinear_Diffusion_step_3D!(Arrays, Mat_tup, Phases, Grid, Num.dt, Num) + # -------------------------------------------- + + # Update variables --------------------------- + # copy to cpu + Tnew_cpu .= Array(Arrays.Tnew) + Phi_melt_cpu .= Array(Arrays.ϕ) + + UpdateTracers_T_ϕ!(Tracers, Grid.coord1D, Tnew_cpu, Phi_melt_cpu); # Update info on tracers + + # copy back to gpu + Arrays.Tnew .= Data.Array(Tnew_cpu) + Arrays.ϕ .= Data.Array(Phi_melt_cpu) + + @parallel assign!(Arrays.T, Arrays.Tnew) + @parallel assign!(Arrays.Tnew, Arrays.T) + # -------------------------------------------- + + # Update info on tracers --------------------- + Tracers = MTK_GMG.MTK_updateTracers(Grid, Arrays, Tracers, Dikes, time_props, Num); + # -------------------------------------------- + + # Update time-dependent properties ----------- + MTK_GMG.MTK_update_TimeDepProps!(time_props, Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + # Visualize results -------------------------- + MTK_GMG.MTK_visualize_output(Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + # Save output to disk once in a while -------- + MTK_GMG.MTK_save_output(Grid, Arrays, Tracers, Dikes, time_props, Num, CartData_input); + # -------------------------------------------- + + # Optionally update arrays and structs (such as T or Dike) ------- + MTK_GMG.MTK_update_ArraysStructs!(Arrays, Grid, Dikes, Num) + # -------------------------------------------- + + # Display output ----------------------------- + MTK_GMG.MTK_print_output(Grid, Num, Arrays, Mat_tup, Dikes) + # -------------------------------------------- + + end + + # Finalize simulation ------------------------ + MTK_GMG.MTK_finalize!(Arrays, Grid, Num, Tracers, Dikes, CartData_input); + # -------------------------------------------- + + + return Grid, Arrays, Tracers, Dikes, time_props +end # end of main function + + +end \ No newline at end of file diff --git a/src/MTK_GMG_structs.jl b/src/MTK_GMG_structs.jl new file mode 100644 index 0000000..b5d1adc --- /dev/null +++ b/src/MTK_GMG_structs.jl @@ -0,0 +1,214 @@ +using GeophysicalModelGenerator + +""" +This mutable structure represents numerical parameters in the program. It is used to store and manage numerical values that are used throughout the program. + mutable struct NumParam <: NumericalParameters + +# Fields + +- `SimName::String`: Name of the simulation. +- `FigTitle`: Title of the figure. +- `Nx::Int64`: Number of grid points in the x direction. +- `Nz::Int64`: Number of grid points in the z direction. +- `W::Float64`: Width of the domain. +- `H::Float64`: Height of the domain. +- `dx::Float64`: Grid spacing in the x direction. +- `dz::Float64`: Grid spacing in the z direction. +- `Tsurface_Celcius::Float64`: Surface temperature in Celsius. +- `Geotherm::Float64`: Geothermal gradient in K/m. +- `maxTime_Myrs::Float64`: Maximum simulation time in Myrs. +- `SecYear::Float64`: Number of seconds in a year. +- `maxTime::Float64`: Maximum simulation time in seconds. +- `SaveOutput_steps::Int64`: Number of steps between output saves. +- `CreateFig_steps::Int64`: Number of steps between figure creations. +- `flux_bottom_BC::Bool`: Whether to apply a flux at the bottom boundary. +- `flux_bottom::Float64`: Flux at the bottom boundary in W/m^2. +- `plot_tracers::Bool`: Whether to plot tracers. +- `advect_polygon::Bool`: Whether to advect a polygon around the intrusion area. +- `axisymmetric::Bool`: Whether the simulation is axisymmetric. +- `κ_time::Float64`: Thermal diffusivity. +- `fac_dt::Float64`: Factor to multiply the time step by. +- `dt::Float64`: Time step. +- `time::Float64`: Current time. +- `nt::Int64`: Total number of time steps. +- `it::Int64`: Current iteration. +- `ω::Float64`: Relaxation parameter for nonlinear iterations. +- `max_iter::Int64`: Maximum number of nonlinear iterations. +- `verbose::Bool`: Whether to print verbose output. +- `convergence::Float64`: Convergence criterion for nonlinear iterations. +- `deactivate_La_at_depth::Bool`: Whether to deactivate latent heating at the bottom of the model box. +- `deactivationDepth::Float64`: Depth at which to deactivate latent heating. +- `USE_GPU`: Whether to use a GPU. +- `AnalyticalInitialGeo::Bool`: Whether to use an analytical initial geotherm. +- `qs_anal::Float64`: Analytical surface heat flux. +- `qm_anal::Float64`: Analytical mantle heat flux. +- `hr_anal::Float64`: Analytical radiogenic heat production. +- `k_anal::Float64`: Analytical thermal conductivity. +- `InitialEllipse::Bool`: Whether to initialize with an ellipse. +- `a_init::Float64`: Semi-major axis of initial ellipse. +- `b_init::Float64`: Semi-minor axis of initial ellipse. +- `TrackTracersOnGrid::Bool`: Whether to track tracers on the grid. + +# Examples + +```julia +np = NumParam(SimName="MySim", Nx=101, Nz=101, ...) +``` + +""" +@with_kw mutable struct NumParam <: NumericalParameters + SimName::String = "Zassy_UCLA_ellipticalIntrusion" # name of simulation + FigTitle = "UCLA setup" + Nx::Int64 = 201 + Ny::Int64 = 0 + Nz::Int64 = 201 + dim::Int64 = length([Nx, Ny, Nz].>0) + W::Float64 = 20e3 + L::Float64 = 0 + H::Float64 = 20e3 + dx::Float64 = W/(Nx-1) + dy::Float64 = L/(Ny-1) + dz::Float64 = H/(Nz-1) # grid spacing in z + Tsurface_Celcius::Float64 = 0 # Surface T in celcius + Geotherm::Float64 = 40/1e3 # in K/m + maxTime_Myrs::Float64 = 1.5 # maximum timestep + SecYear::Float64 = 3600*24*365.25; + maxTime::Float64 = maxTime_Myrs*SecYear*1e6 # maximum timestep in seconds + flux_bottom_BC::Bool = false # flux bottom BC? + flux_bottom::Float64 = 167e-3 # Flux in W/m2 in case flux_bottom_BC=true + plot_tracers::Bool = true # adds passive tracers to the plot + advect_polygon::Bool = false # adds a polygon around the intrusion area + axisymmetric::Bool = false # axisymmetric (if true) of 2D geometry? + κ_time::Float64 = 3.3/(1000*2700) # κ to determine the stable timestep + fac_dt::Float64 = 0.4; # prefactor with which dt is multiplied + Δ::Vector{Float64} = [dx, dy, dz]; # grid spacing + Δmin::Float64 = minimum(Δ[Δ.>0]); # minimum grid spacing + dt::Float64 = fac_dt*(Δmin^2)./κ_time/4; # timestep + time::Float64 = 0.0; # current time + nt::Int64 = floor(maxTime/dt); + it::Int64 = 0; # current iteration + ω::Float64 = 0.8; # relaxation parameter for nonlinear iterations + max_iter::Int64 = 5000; # max. number of nonlinear iterations + verbose::Bool = false; + convergence::Float64 = 1e-5; # nonlinear convergence criteria + USE_GPU::Bool = false; + keep_init_RockPhases::Bool = true; # keep initial rock phases (if false, all phases are initialized as Dikes.BackgroundPhase) + pvd::Union{Nothing,GeophysicalModelGenerator.WriteVTK.CollectionFile} = nothing; # pvd file info for paraview + Output_VTK::Bool = true; # output VTK files in case CartData is an input? + SaveOutput_steps::Int64 = 1e3; # saves output every x steps + CreateFig_steps::Int64 = 500; # Create a figure every X steps + + AddRandomSills::Bool = false; # Add random sills/dikes to the model? + RandomSills_timestep::Int64 = 10; # After how many timesteps do we add a new sill/dike? + + # parts that can be removed @ some stage + deactivate_La_at_depth::Bool= false # deactivate latent heating @ the bottom of the model box? + deactivationDepth::Float64 = -15e3 # deactivation depth + AnalyticalInitialGeo::Bool = false; + qs_anal::Float64 = 170e-3; + qm_anal::Float64 = 167e-3; + hr_anal::Float64 = 10e3; + k_anal::Float64 = 3.35; + InitialEllipse::Bool = false; + a_init::Float64 = 2.5e3; + b_init::Float64 = 1.5e3; + TrackTracersOnGrid::Bool = true; +end + +""" + mutable struct DikeParam <: DikeParameters + +This mutable structure represents parameters related to a dike in the simulation. It is used to store and manage values related to the dike's properties and behavior. + +# Fields + +- `Type::String`: Type of the dike. +- `Center::Vector{Float64}`: Center of the dike. +- `T_in_Celsius::Float64`: Temperature of the injected magma in Celsius. +- `W_in::Float64`: Diameter of the dike. +- `H_in::Float64`: Thickness of the dike. +- `AspectRatio::Float64`: Aspect ratio of the dike. +- `SillRadius::Float64`: Radius of the sill. +- `SillArea::Float64`: Horizontal area of the sill. +- `InjectionInterval_year::Float64`: Injection interval in years. +- `SecYear`: Number of seconds in a year. +- `InjectionInterval::Float64`: Injection interval in seconds. +- `nTr_dike::Int64`: Number of tracers in the dike. +- `InjectVol`: Injected volume into the dike. +- `Qrate_km3_yr`: Dike insertion rate in km^3/year. +- `dike_poly`: Polygon representing the dike. +- `dike_inj`: Injection into the dike. + +- `H_ran`: Zone in which we vary the vertical location of the dike (if we add random dikes) +- `L_ran`: Zone in which we vary the horizontal (x) location of the dike (if we add random dikes) +- `W_ran`: Zone in which we vary the horizontal (y) location of the dike (if we add random dikes) + +- `Dip_ran`: maximum variation of dip (if we add random dikes) +- `Strike_ran`: maximum variation of strike (if we add random dikes) + + +# Examples + +```julia +dp = DikeParam(Type="MyDike", Center=[0., -7.0e3], ...) +``` +""" +@with_kw mutable struct DikeParam <: DikeParameters + Type::String = "CylindricalDike_TopAccretion" + Center::Vector{Float64} = [0.; -7.0e3 - 0/2]; # Center of dike + Angle::Vector{Float64} = [0.0]; # Angle of dike + T_in_Celsius::Float64 = 1000; # Temperature of injected magma + W_in::Float64 = 20e3 # Diameter of dike + H_in::Float64 = 74.6269 # Thickness + AspectRatio::Float64 = H_in/W_in; # Aspect ratio + SillRadius::Float64 = W_in/2 # Sill radius + SillArea::Float64 = pi*SillRadius^2 # Horizontal area of sill + InjectionInterval_year::Float64 = 10e3; # Injection interval [years] + SecYear = 3600*24*365.25; # s/year + InjectionInterval::Float64 = InjectionInterval_year*SecYear; # Injection interval [s] + nTr_dike::Int64 = 300 # Number of tracers + InjectVol::Float64 = 0.0; # injected volume + Qrate_km3_yr::Float64 = 0.0; # Dikes insertion rate + BackgroundPhase::Int64 = 1; # Background phase (non-dikes) + DikePhase::Int64 = 2; # Dike phase + dike_poly::Vector = []; # polygon with dike + dike_inj::Float64 = 0.0 + + H_ran::Float64 = 5000.0 # Zone in which we vary the horizontal location of the dike + L_ran::Float64 = 2000.0 # Zone in which we vary the horizontal location of the dike + W_ran::Float64 = 2000.0 # Zone in which we vary the vertical location of the dike + Dip_ran::Float64 = 30.0; # maximum variation of dip + Strike_ran::Float64 = 90.0; # maximum variation of strike + SillsAbove::Float64 = -15e3; # Sills above this depth + +end + + +""" + mutable struct TimeDepProps <: TimeDependentProperties + +This mutable structure represents time-dependent properties in the simulation. It is used to store and manage values that change over time. + +# Fields + +- `Time_vec::Vector{Float64}`: Vector storing the time points. +- `MeltFraction::Vector{Float64}`: Vector storing the melt fraction at each time point. +- `Tav_magma::Vector{Float64}`: Vector storing the average magma temperature at each time point. +- `Tmax::Vector{Float64}`: Vector storing the maximum magma temperature at each time point. + +# Examples + +```julia +tdp = TimeDepProps(Time_vec=[0., 1., 2.], MeltFraction=[0.1, 0.2, 0.3], ...) +``` + +# Note: +You can use multiple dispatch on this struct in your user code as long as the new struct + +""" +@with_kw mutable struct TimeDepProps <: TimeDependentProperties + Time_vec::Vector{Float64} = []; # Center of dike + MeltFraction::Vector{Float64} = []; # Melt fraction over time + Tav_magma::Vector{Float64} = []; # Average magma + Tmax::Vector{Float64} = []; # Max magma temperature +end diff --git a/src/MagmaThermoKinematics.jl b/src/MagmaThermoKinematics.jl index d48698e..a4a9b0c 100644 --- a/src/MagmaThermoKinematics.jl +++ b/src/MagmaThermoKinematics.jl @@ -17,8 +17,21 @@ using JLD2 # Load/save data to disk @reexport using GeoParams # Material parameters calculations @reexport using ParallelStencil +abstract type NumericalParameters end +abstract type DikeParameters end +abstract type TimeDependentProperties end + include("Units.jl") # various useful units +# Few useful parameters +const SecYear = 3600*24*365.25 +const kyr = 1000*SecYear +const Myr = 1e6*SecYear +const km³ = 1000^3 +export SecYear, kyr, Myr, km³ + +export NumericalParameters, DikeParameters, TimeDependentProperties + function environment!(model_device, precision, dimension) gpu = model_device == :gpu ? true : false @@ -54,7 +67,7 @@ function environment!(model_device, precision, dimension) k = keys(args) v = getindex.(values(args), i, j) argsi = (; zip(k, v)...) - A[i, j] = $(_fn)(MatParam[Phases[i,j]], argsi) + A[i, j] = $(_fn)(MatParam, Phases[i,j], argsi) return end @@ -121,21 +134,38 @@ function environment!(model_device, precision, dimension) export Process_ZirconAges, copy_arrays_GPU2CPU!, copy_arrays_CPU2GPU! end + # GMG integration + Base.@eval begin + include(joinpath(@__DIR__, "MTK_GMG_structs.jl")) + export NumParam, DikeParam, TimeDepProps + + include(joinpath(@__DIR__, "MTK_GMG.jl")) + + include(joinpath(@__DIR__, "MTK_GMG_2D.jl")) + using .MTK_GMG_2D + export MTK_GeoParams_2D + + include(joinpath(@__DIR__, "MTK_GMG_3D.jl")) + using .MTK_GMG_3D + export MTK_GeoParams_3D + + end + end export environment! -# Few useful parameters -SecYear = 3600*24*365.25 -kyr = 1000*SecYear -Myr = 1e6*SecYear -km³ = 1000^3 -export SecYear, kyr, Myr, km³ - include("Grid.jl") using .Grid export GridData, CreateGrid +# Routines that deal with tracers +include("Tracers.jl") +export UpdateTracers, AdvectTracers!, InitializeTracers,PhaseRatioFromTracers, CorrectTracersForTopography! +export RockAssemblage, update_Tvec! +export PhaseRatioFromTracers!, PhasesFromTracers!, UpdateTracers_T_ϕ!, UpdateTracers_Field! # new routines + + include("MeltingRelationships.jl") export SolidFraction, ComputeLithostaticPressure, LoadPhaseDiagrams, PhaseDiagramData, ComputeDensityAndPressure export PhaseRatioAverage!, ComputeSeismicVelocities, SolidFraction_Parameterized! @@ -153,13 +183,6 @@ export Tracer, AddDike, HostRockVelocityFromDike, CreateDikePolygon, advect_dike include("Advection.jl") export AdvectTemperature, Interpolate!, CorrectBounds, evaluate_interp_2D, evaluate_interp_3D -# Routines that deal with tracers -include("Tracers.jl") -export UpdateTracers, AdvectTracers!, InitializeTracers,PhaseRatioFromTracers, CorrectTracersForTopography! -export RockAssemblage, update_Tvec! -export PhaseRatioFromTracers!, PhasesFromTracers!, UpdateTracers_T_ϕ!, UpdateTracers_Field! # new routines - - # Routines related to Parameters.jl, which come in handy in the main routine export @unpack, @with_kw diff --git a/test/test_MTK_GMG_2D.jl b/test/test_MTK_GMG_2D.jl new file mode 100644 index 0000000..3a0abf2 --- /dev/null +++ b/test/test_MTK_GMG_2D.jl @@ -0,0 +1,180 @@ +using Test +#using Plots + +const USE_GPU=false; + +using MagmaThermoKinematics +if USE_GPU + environment!(:gpu, Float64, 2) # initialize parallel stencil in 2D +else + environment!(:cpu, Float64, 2) # initialize parallel stencil in 2D +end +using MagmaThermoKinematics +using MagmaThermoKinematics.Diffusion2D +using MagmaThermoKinematics.MTK_GMG_2D + +using Random, GeoParams, GeophysicalModelGenerator + +const rng = Random.seed!(1234); # same seed such that we can reproduce results + +# Import a few routines, so we can overwrite them below +using MagmaThermoKinematics.MTK_GMG + +@testset "MTK_GMG_2D" begin +#= + function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + println("$(Num.it), Time=$(round(Num.time/Num.SecYear)) yrs; max(T) = $(round(maximum(Arrays.Tnew)))") + return nothing + end +=# +# Test setup +println("===============================================") +println("Testing the MTK - GMG integration") +println("===============================================") + +# These are the final simulations for the ZASSy paper, but done @ a lower resolution +Num = NumParam( #Nx=269*1, Nz=269*1, + Nx=135*1, Nz=135*1, + SimName="ZASSy_Geneva_9_1e_6", axisymmetric=false, + #maxTime_Myrs=1.5, + maxTime_Myrs=0.025, + fac_dt=0.2, ω=0.5, verbose=false, + flux_bottom_BC=false, flux_bottom=0, deactivate_La_at_depth=false, + Geotherm=30/1e3, TrackTracersOnGrid=true, + SaveOutput_steps=100000, CreateFig_steps=100000, plot_tracers=false, advect_polygon=true, + FigTitle="Geneva Models, Geotherm 30/km", + USE_GPU=USE_GPU); + +Dike_params = DikeParam(Type="CylindricalDike_TopAccretion", + InjectionInterval_year = 5000, # flux= 14.9e-6 km3/km2/yr + W_in=20e3, H_in=74.6269, + nTr_dike=300*1 + ) + +MatParam = (SetMaterialParams(Name="Rock & partial melt", Phase=1, + Density = ConstantDensity(ρ=2700kg/m^3), + LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg), + #LatentHeat = ConstantLatentHeat(Q_L=0.0J/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), + 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 + ) + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_params); # start the main code + +@test sum(Arrays.Tnew)/prod(size(Arrays.Tnew)) ≈ 315.46382940863816 rtol= 1e-2 +@test sum(time_props.MeltFraction) ≈ 0.3211217281417583 rtol= 1e-5 + +# ----------------------------- + +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) +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)) + +# Intersect with topography +Below = BelowSurface(Data_2D, Topo_cart) +Data_2D.fields.Phases[Below] .= 1 + +# Set Moho +ind = findall(Data_2D.z.val .< -30.0) +Data_2D.fields.Phases[ind] .= 2 + +# Set T: +gradient = 30 +Data_2D.fields.Temp .= -Data_2D.z.val*gradient +ind = findall(Data_2D.fields.Temp .< 10.0) +Data_2D.fields.Temp[ind] .= 10.0 + +# Set thermal anomaly +x_c, z_c, r = -10, -15, 2.5 +Volume = 4/3*pi*r^3 # equivalent 3D volume of the anomaly [km^3] +ind = findall((Data_2D.x.val .- x_c).^2 .+ (Data_2D.z.val .- z_c).^2 .< r^2) +Data_2D.fields.Temp[ind] .= 800.0 + +""" +Randomly change orientation and location of a dike +""" +function MTK_GMG.MTK_update_ArraysStructs!(Arrays::NamedTuple, Grid::GridData, Dikes::DikeParameters, Num::NumericalParameters) + if mod(Num.it,10)==0 + cen = (Grid.max .+ Grid.min)./2 .+ 0*rand(rng, -0.5:1e-3:0.5, 2).*[Dikes.W_ran; Dikes.H_ran]; # Randomly vary center of dike + if cen[end]<-15e3; Angle_rand = 0*rand(rng, 80.0:0.1:100.0) # Orientation: near-vertical @ depth + else Angle_rand = 0*rand(rng,-10.0:0.1:10.0); end + + Dikes.Center = cen; + Dikes.Angle = [Angle_rand]; + end + return nothing +end + +# Define numerical parameters +Num = NumParam( SimName="Unzen1", axisymmetric=false, + maxTime_Myrs=0.005, + fac_dt=0.2, ω=0.5, verbose=false, + SaveOutput_steps=10000, CreateFig_steps=1000, plot_tracers=false, advect_polygon=false, + USE_GPU=USE_GPU); + +# dike parameters +Dike_params = DikeParam(Type="ElasticDike", + InjectionInterval_year = 1000, # flux= 14.9e-6 km3/km2/yr + W_in=5e3, H_in=250, + nTr_dike=300*1, + H_ran = 5000, W_ran = 5000, + DikePhase=3, BackgroundPhase=1, + ) + +# Define parameters for the different phases +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), + 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), + 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)), + + 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), + Melting = SmoothMelting(MeltingParam_4thOrder())) # Marxer & Ulmer melting + + ) + + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_params, CartData_input=Data_2D); # start the main code + +@test sum(Arrays.Tnew)/prod(size(Arrays.Tnew)) ≈ 251.5482011114283 rtol= 1e-2 +@test sum(time_props.MeltFraction) ≈ 1.0066096298950455 rtol= 1e-5 + + +end diff --git a/test/test_MTK_GMG_3D.jl b/test/test_MTK_GMG_3D.jl new file mode 100644 index 0000000..c66d5a4 --- /dev/null +++ b/test/test_MTK_GMG_3D.jl @@ -0,0 +1,170 @@ +using Test +#using Plots + +const USE_GPU=false; + +using MagmaThermoKinematics +if USE_GPU + environment!(:gpu, Float64, 3) # initialize parallel stencil in 2D + CUDA.device!(1) # select the GPU you use (starts @ zero) +else + environment!(:cpu, Float64, 3) # initialize parallel stencil in 2D +end +using MagmaThermoKinematics.Diffusion3D + +using Random, GeoParams, GeophysicalModelGenerator + +const rng = Random.seed!(1234); # same seed such that we can reproduce results + +# Allow overwriting user routines +using MagmaThermoKinematics.MTK_GMG + +@testset "MTK_GMG_3D" begin + +function MTK_GMG.MTK_print_output(Grid::GridData, Num::NumericalParameters, Arrays::NamedTuple, Mat_tup::Tuple, Dikes::DikeParameters) + if mod(Num.it,10) == 0 + println("$(Num.it), $(Num.time/SecYear/1e3) kyrs; max(T)=$(maximum(Arrays.Tnew))") + end + return nothing +end + +# Test setup +println("===============================================") +println("Testing MTK - GMG integration in 3D") +println("===============================================") + +# These are the final simulations for the ZASSy paper, but done @ a lower resolution +Num = NumParam( #Nx=269*1, Nz=269*1, + Nx=65*1, Ny=65*1, Nz=65*1, + SimName="Test1", + W=20e3, H=20e3, L=20e3, + #maxTime_Myrs=1.5, + maxTime_Myrs=0.005, + fac_dt=0.2, ω=0.5, verbose=false, + flux_bottom_BC=false, flux_bottom=0, deactivate_La_at_depth=false, + Geotherm=30/1e3, TrackTracersOnGrid=true, + SaveOutput_steps=10, CreateFig_steps=100000, plot_tracers=false, advect_polygon=true, + FigTitle="Geneva Models, Geotherm 30/km", + USE_GPU=USE_GPU, + AddRandomSills = false, RandomSills_timestep=5 + ); + +Dike_params = DikeParam(Type="ElasticDike", + InjectionInterval_year = 1000, + W_in=5e3, H_in=200.0*4, # note: H must be numerically resolved + Dip_ran = 20.0, Strike_ran = 0.0, + W_ran = 10e3; H_ran = 10e3, L_ran=10e3, + nTr_dike=300*1, + SillsAbove = -10e3, + Center=[0.0,0.0, -7000], Angle=[0.0, 0.0], + ) + +MatParam = (SetMaterialParams(Name="Rock & partial melt", Phase=1, + Density = ConstantDensity(ρ=2700kg/m^3), + LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg), + #LatentHeat = ConstantLatentHeat(Q_L=0.0J/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), + 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 + ) + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_3D(MatParam, Num, Dike_params); # start the main code + +@test sum(Arrays.Tnew)/prod(size(Arrays.Tnew)) ≈ 303.1195760751668 rtol= 1e-2 +@test sum(time_props.MeltFraction) ≈ 3.2384408012256802 rtol= 1e-5 +# ----------------------------- + + +Topo_cart = load_GMG("../examples/Topo_cart") # Note: Laacher seee is around [10,20] + +# Create 3D grid of the region +Nx,Ny,Nz = 100,100,100 +X,Y,Z = XYZGrid(range(-23,23, length=Nx),range(-19,19, length=Ny),range(-20,5, length=Nz)) +Data_3D = CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X)))); # 3D dataset + +# Intersect with topography +Below = BelowSurface(Data_3D, Topo_cart) +Data_3D.fields.Phases[Below] .= 1 + +# Set Moho +ind = findall(Data_3D.z.val .< -30.0) +Data_3D.fields.Phases[ind] .= 2 + +# Set T: +gradient = 30 +Data_3D.fields.Temp .= -Data_3D.z.val*gradient +ind = findall(Data_3D.fields.Temp .< 10.0) +Data_3D.fields.Temp[ind] .= 10.0 + +# Set thermal anomaly +x_c, y_c, z_c, r = -10, -10, -15, 2.5 +Volume = 4/3*pi*r^3 # equivalent 3D volume of the anomaly [km^3] +ind = findall((Data_3D.x.val .- x_c).^2 .+ (Data_3D.y.val .- y_c).^2 .+ (Data_3D.z.val .- z_c).^2 .< r^2) +Data_3D.fields.Temp[ind] .= 800.0 + + +# Define numerical parameters +Num = NumParam( SimName="Unzen2", axisymmetric=false, + maxTime_Myrs=0.005, + fac_dt=0.2, + SaveOutput_steps=20, CreateFig_steps=1000, plot_tracers=false, advect_polygon=false, + USE_GPU=USE_GPU, + AddRandomSills = false, RandomSills_timestep=5); + +# dike parameters +Dike_params = DikeParam(Type="ElasticDike", + InjectionInterval_year = 1000, # flux= 14.9e-6 km3/km2/yr + W_in=5e3, H_in=250*4, + nTr_dike=300*1, + H_ran = 5000, W_ran = 5000, + DikePhase=3, BackgroundPhase=1, + Center=[0.0,0.0, -7000], Angle=[0.0, 0.0], + ) + +# Define parameters for the different phases +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), + 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), + 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)), + + 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), + Melting = SmoothMelting(MeltingParam_4thOrder())) # Marxer & Ulmer melting + + ) + + +# Call the main code with the specified material parameters +Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_3D(MatParam, Num, Dike_params, CartData_input=Data_3D); # start the main code + +@test sum(Arrays.Tnew)/prod(size(Arrays.Tnew)) ≈ 244.14916470514495 rtol= 1e-2 +@test sum(time_props.MeltFraction) ≈ 5.706560313331254 rtol= 1e-5 + + +end