Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean up subduction miniapp #168

Merged
merged 1 commit into from
May 30, 2024
Merged

Clean up subduction miniapp #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
74 changes: 4 additions & 70 deletions miniapps/subduction/2D/Subd2D.jl
View file
Open in desktop
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
Open in desktop
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
Loading