Skip to content

Commit

Permalink
clean up subduction miniapp (#168)
Browse files Browse the repository at this point in the history
  • Loading branch information
@albert-de-montserrat
albert-de-montserrat authored May 30, 2024
1 parent 559c13b commit bb04877
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 102 deletions.
74 changes: 4 additions & 70 deletions miniapps/subduction/2D/Subd2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,16 +75,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
# ----------------------------------------------------

# Physical properties using GeoParams ----------------
ρbg = 3.2e3 + 1
rheology = init_rheology_linear(; ρbg = ρbg)
rheology_augmented = init_rheology_linear(; ρbg = 0e0)

rheology = init_rheology_nonNewtonian(; ρbg = ρbg)
rheology_augmented = init_rheology_nonNewtonian(; ρbg = 0e0)

rheology = init_rheology_nonNewtonian_plastic(; ρbg = ρbg)
rheology_augmented = init_rheology_nonNewtonian_plastic(; ρbg = 0e0)

rheology = init_rheology_linear()
dt = 50e3 * 3600 * 24 * 365 # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down Expand Up @@ -133,15 +124,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk

# Buoyancy forces
ρg = ntuple(_ -> @zeros(ni...), Val(2))
# for _ in 1:2
# compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P))
# JustRelax.Stokes2D.init_P!(stokes.P, ρg[2], xci[2])
# end
compute_ρg!(ρg[2], phase_ratios, rheology_augmented, (T=thermal.Tc, P=stokes.P))
stokes.P .= PTArray(backend)(reverse(cumsum(reverse((ρg[2]).* di[2], dims=2), dims=2), dims=2))

# ρg_bg = ρbg * 9.81
# Plitho = reverse(cumsum(reverse((ρg[2] .+ ρg_bg).* di[2], dims=2), dims=2), dims=2)

# Rheology
args0 = (T=thermal.Tc, P=stokes.P, dt = Inf)
Expand All @@ -160,24 +144,6 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
)
flow_bcs!(stokes, flow_bcs) # apply boundary conditions
update_halo!(@velocity(stokes)...)

# Plot initial T and η profiles
# let
# Yv = [y for x in xvi[1], y in xvi[2]][:]
# Y = [y for x in xci[1], y in xci[2]][:]
# fig = Figure(size = (1200, 900))
# ax1 = Axis(fig[1,1], aspect = 2/3, title = "T")
# ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
# scatter!(ax1, Array(thermal.T[2:end-1,:][:]), Yv./1e3)
# scatter!(ax2, Array(log10.(stokes.viscosity.η[:])), Y./1e3)
# # scatter!(ax2, Array(stokes.P[:]), Y./1e3)
# # scatter!(ax2, Array(Plitho[:]), Y./1e3)
# ylims!(ax1, minimum(xvi[2])./1e3, 0)
# ylims!(ax2, minimum(xvi[2])./1e3, 0)
# hideydecorations!(ax2)
# # save(joinpath(figdir, "initial_profile.png"), fig)
# fig
# end

# IO -------------------------------------------------
# if it does not exist, make folder where figures are stored
Expand Down Expand Up @@ -225,12 +191,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
# Time loop
t, it = 0.0, 0

# Vertice arrays of normal components of the stress tensor
# τxx_vertex, τyy_vertex = @zeros(ni.+1...), @zeros(ni.+1...)
# τxx_o_vertex, τyy_o_vertex = @zeros(ni.+1...), @zeros(ni.+1...)

while it < 1 #000 # run only for 5 Myrs
# while (t/(1e6 * 3600 * 24 *365.25)) < 5 # run only for 5 Myrs
while it < 1000 # run only for 5 Myrs

# interpolate fields from particle to grid vertices
particle2grid!(T_buffer, pT, xvi, particles)
Expand All @@ -240,15 +201,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
thermal_bcs!(thermal, thermal_bc)
temperature2center!(thermal)

# Update buoyancy and viscosity -
# Plitho .= reverse(cumsum(reverse((ρg[2]).* di[2], dims=2), dims=2), dims=2)
# Plitho .= stokes.P .+ Plitho .- minimum(stokes.P)

# args = (; T = thermal.Tc, P = Plitho, dt=Inf, ρbg = ρbg * 9.81)
args = (; T = thermal.Tc, P = stokes.P, dt=Inf)
# compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff)
# compute_ρg!(ρg[2], phase_ratios, rheology, args)


# Stokes solver ----------------
t_stokes = @elapsed begin
out = solve!(
Expand Down Expand Up @@ -309,24 +263,9 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
advection!(particles, RungeKutta2(), @velocity(stokes), grid_vxi, dt)
# advect particles in memory
move_particles!(particles, xvi, particle_args)
# JustRelax.Stokes2D.rotate_stress_particles!(
# stokes,
# τxx_vertex,
# τyy_vertex,
# τxx_o_vertex,
# τyy_o_vertex,
# τxx_p,
# τyy_p,
# τxy_p,
# vorticity_p,
# particles,
# grid,
# dt
# )

# check if we need to inject particles
inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer, ), xvi)
# inject_particles_phase!(particles, pPhases, (pT, τxx_p, τyy_p, τxy_p, vorticity_p), (T_buffer, τxx_vertex, τyy_vertex, stokes.τ.xy, stokes.ω.xy_v), xvi)
# update phase ratios
phase_ratios_center!(phase_ratios, particles, grid, pPhases)

Expand Down Expand Up @@ -408,14 +347,9 @@ end

## END OF MAIN SCRIPT ----------------------------------------------------------------
do_vtk = true # set to true to generate VTK files for ParaView
figdir = "Buiter_2D"
# nx, ny = 512, 256
# nx, ny = 512, 128
figdir = "Subduction2D"
n = 128
nx, ny = n*6, n
nx, ny = 512, 256
# nx, ny = 128, 128
# nx, ny = 32*6, 32
li, origin, phases_GMG, T_GMG = GMG_subduction_2D(nx+1, ny+1)
igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
IGG(init_global_grid(nx, ny, 1; init_MPI= true)...)
Expand Down
43 changes: 11 additions & 32 deletions miniapps/subduction/2D/Subd2D_rheology.jl
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
using GeoParams.Dislocation
using GeoParams.Diffusion

function init_rheology_nonNewtonian(; ρbg = 0e0)
function init_rheology_nonNewtonian()
#dislocation laws
disl_wet_olivine = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
# diffusion laws
diff_wet_olivine = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003)

lithosphere_rheology = CompositeRheology( (disl_wet_olivine, diff_wet_olivine))
init_rheologies(lithosphere_rheology; ρbg = ρbg)
init_rheologies(lithosphere_rheology)
end

function init_rheology_nonNewtonian_plastic(; ρbg = 0e0)
function init_rheology_nonNewtonian_plastic()
#dislocation laws
disl_wet_olivine = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
# diffusion laws
Expand All @@ -28,29 +28,15 @@ function init_rheology_nonNewtonian_plastic(; ρbg = 0e0)
DruckerPrager_regularised(; C = C_wet_olivine, ϕ = ϕ_wet_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
)
)
init_rheologies(lithosphere_rheology; ρbg = ρbg)
init_rheologies(lithosphere_rheology)
end

function init_rheology_linear(; ρbg = 0e0)
function init_rheology_linear()
lithosphere_rheology = CompositeRheology( (LinearViscous(; η=1e20),))
init_rheologies(lithosphere_rheology; ρbg = ρbg)
init_rheologies(lithosphere_rheology)
end

function init_rheologies(lithosphere_rheology; ρbg = 0e0)
# #dislocation laws
# disl_wet_olivine = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
# # diffusion laws
# diff_wet_olivine = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003)

# ϕ_wet_olivine = asind(0.1)
# C_wet_olivine = 1e6

# ϕ_oceanic_crust_upper = asind(0.1)
# C_oceanic_crust_upper = 0.3e6

# # soft_C = LinearSoftening((C_oceanic_litho*0.05, C_oceanic_litho), (0.1, 0.5))

# elasticity = ConstantElasticity(; G=5e10, ν=0.5)
function init_rheologies(lithosphere_rheology)
# common physical properties
α = 2.4e-5 # 1 / K
Cp = 750 # J / kg K
Expand All @@ -60,7 +46,7 @@ function init_rheologies(lithosphere_rheology; ρbg = 0e0)
# Name = "Asthenoshpere",
SetMaterialParams(;
Phase = 1,
Density = ConstantDensity(; ρ=3.2e3-ρbg),
Density = ConstantDensity(; ρ=3.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e20),)),
Expand All @@ -69,30 +55,23 @@ function init_rheologies(lithosphere_rheology; ρbg = 0e0)
# Name = "Oceanic lithosphere",
SetMaterialParams(;
Phase = 2,
Density = PT_Density(; ρ0=3.2e3-ρbg, α = α, β = 0e0, T0 = 273+1474),
Density = PT_Density(; ρ0=3.2e3, α = α, β = 0e0, T0 = 273+1474),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
# CompositeRheology = CompositeRheology(
# (
# disl_wet_olivine,
# diff_wet_olivine,
# DruckerPrager_regularised(; C = C_wet_olivine, ϕ = ϕ_wet_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
# )
# ),
CompositeRheology = lithosphere_rheology,
),
# Name = "oceanic crust",
SetMaterialParams(;
Phase = 3,
Density = ConstantDensity(; ρ=3.2e3-ρbg),
Density = ConstantDensity(; ρ=3.2e3),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e20),)),
),
# Name = "StickyAir",
SetMaterialParams(;
Phase = 4,
Density = ConstantDensity(; ρ=100-ρbg), # water density
Density = ConstantDensity(; ρ=100), # water density
HeatCapacity = ConstantHeatCapacity(; Cp=3e3),
Conductivity = ConstantConductivity(; k=1.0),
CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
Expand Down

0 comments on commit bb04877

Please sign in to comment.