Skip to content

Commit

Permalink
up
Browse files Browse the repository at this point in the history
  • Loading branch information
@albert-de-montserrat
albert-de-montserrat committed Nov 26, 2024
1 parent e197839 commit c4abfa7
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 79 deletions.
155 changes: 83 additions & 72 deletions volcano2D/Caldera2D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk

# Physical properties using GeoParams ----------------
rheology = init_rheologies()
dt = 1e2 * 3600 * 24 * 365 # diffusive CFL timestep limiter
dt = 5e2 * 3600 * 24 * 365 # diffusive CFL timestep limiter
# dt = Inf # diffusive CFL timestep limiter
# ----------------------------------------------------

Expand Down Expand Up @@ -138,25 +138,24 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
# Allocate arrays needed for every Stokes problem
stokes = StokesArrays(backend, ni)
pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-4, Re = 3e0, r=0.7, CFL = 0.98 / √2.1) # Re=3π, r=0.7
# pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-5, Re = 2π√2, r=0.7, CFL = 0.9 / √2.1) # Re=3π, r=0.7
# ----------------------------------------------------

# TEMPERATURE PROFILE --------------------------------
Ttop = T_GMG[1, end]
Tbot = T_GMG[1, 1]
thermal = ThermalArrays(backend, ni)
@views thermal.T[2:end-1, :] .= PTArray(backend)(T_GMG)
thermal_bc = TemperatureBoundaryConditions(;
no_flux = (left = true, right = true, top = false, bot = false),
)
thermal_bcs!(thermal, thermal_bc)
temperature2center!(thermal)
Ttop = thermal.T[2:end-1, end]
Tbot = thermal.T[2:end-1, 1]
# ----------------------------------------------------

# Buoyancy forces
ρg = ntuple(_ -> @zeros(ni...), Val(2))
compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P))
# stokes.P .= PTArray(backend)(reverse(cumsum(reverse((ρg[2]).* di[2], dims=2), dims=2), dims=2))
stokes.P .= PTArray(backend)(reverse(cumsum(reverse((ρg[2]).* di[2], dims=2), dims=2), dims=2))

# Rheology
args0 = (T=thermal.Tc, P=stokes.P, dt = Inf)
Expand All @@ -169,12 +168,21 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
)

# Boundary conditions
# flow_bcs = DisplacementBoundaryConditions(;
flow_bcs = VelocityBoundaryConditions(;
free_slip = (left = true , right = true , top = true , bot = true),
free_surface = false,
)
εbg = 1e-15

# U = 0.02
# stokes.U.Ux .= PTArray(backend)([(x - li[1] * 0.5) * U / dt for x in xvi[1], _ in 1:ny+2])
# stokes.U.Uy .= PTArray(backend)([-y * U / dt for _ in 1:nx+2, y in xvi[2]])
# flow_bcs!(stokes, flow_bcs) # apply boundary conditions
# displacement2velocity!(stokes, dt)

εbg = 1e-15 * 1
apply_pure_shear(@velocity(stokes)..., εbg, xvi, li...)

flow_bcs!(stokes, flow_bcs) # apply boundary conditions
update_halo!(@velocity(stokes)...)

Expand Down Expand Up @@ -203,26 +211,26 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
grid2particle!(pT, xvi, T_buffer, particles)

## Plot initial T and P profile
fig = 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="Pressure")
scatter!(
ax1,
Array(thermal.T[2:(end - 1), :][:]),
Yv,
)
lines!(
ax2,
Array(stokes.P[:]),
Y,
)
hideydecorations!(ax2)
# save(joinpath(figdir, "initial_profile.png"), fig)
fig
end
# fig = 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="Pressure")
# scatter!(
# ax1,
# Array(thermal.T[2:(end - 1), :][:]),
# Yv,
# )
# lines!(
# ax2,
# Array(stokes.P[:]),
# Y,
# )
# hideydecorations!(ax2)
# # save(joinpath(figdir, "initial_profile.png"), fig)
# fig
# end

τxx_v = @zeros(ni.+1...)
τyy_v = @zeros(ni.+1...)
Expand All @@ -231,7 +239,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
t, it = 0.0, 0
thermal.Told .= thermal.T

while it < 100 #000 # run only for 5 Myrs
while it < 500 #000 # run only for 5 Myrs

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

args = (; T=thermal.Tc, P=stokes.P, dt=Inf, ΔTc=thermal.ΔTc)
# args = (; T=thermal.Tc, P=stokes.P, dt=Inf)

particle2centroid!(stokes.τ.xx, pτxx, xci, particles)
particle2centroid!(stokes.τ.yy, pτyy, xci, particles)
particle2grid!(stokes.τ.xy, pτxy, xvi, particles)
# args = (; T=thermal.Tc, P=stokes.P, dt=Inf, ΔTc=thermal.ΔTc)
args = (; T=thermal.Tc, P=stokes.P, dt=Inf)

t_stokes = @elapsed solve_VariationalStokes!(
stokes,
Expand Down Expand Up @@ -278,13 +282,13 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
println(" Total time: $t_stokes s")
tensor_invariant!(stokes.ε)
tensor_invariant!(stokes.ε_pl)
dtmax = 1e3 * 3600 * 24 * 365.25
dt = compute_dt(stokes, di) * 0.8
dtmax = 2e3 * 3600 * 24 * 365.25
dt = compute_dt(stokes, di, dtmax)

println("dt = $(dt/(3600 * 24 *365.25)) years")
# ------------------------------

# # Thermal solver ---------------
# Thermal solver ---------------
heatdiffusion_PT!(
thermal,
pt_thermal,
Expand Down Expand Up @@ -332,11 +336,15 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases)
update_rock_ratio!(ϕ, phase_ratios, (phase_ratios.Vx, phase_ratios.Vy), air_phase)

particle2centroid!(stokes.τ_o.xx, pτxx, xci, particles)
particle2centroid!(stokes.τ_o.yy, pτyy, xci, particles)
particle2grid!(stokes.τ_o.xy, pτxy, xvi, particles)

@show it += 1
t += dt

# Data I/O and plotting ---------------------
if it == 1 || rem(it, 2) == 0
if it == 1 || rem(it, 1) == 0
# checkpointing(figdir, stokes, thermal.T, η, t)
(; η_vep, η) = stokes.viscosity
if do_vtk
Expand All @@ -345,6 +353,7 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
T = Array(T_buffer),
τII = Array(stokes.τ.II),
εII = Array(stokes.ε.II),
εII_pl = Array(stokes.ε_pl.II),
)
data_c = (;
P = Array(stokes.P),
Expand All @@ -360,7 +369,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
xci,
data_v,
data_c,
velocity_v
velocity_v;
t = 0
)
end

Expand All @@ -379,19 +389,19 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
ax2 = Axis(fig[2,1], aspect = ar, title = "Vy")
# ax2 = Axis(fig[2,1], aspect = ar, title = "Phase")
ax3 = Axis(fig[1,3], aspect = ar, title = "τII [MPa]")
ax4 = Axis(fig[2,3], aspect = ar, title = "log10(εII_plastic)")
# ax4 = Axis(fig[2,3], aspect = ar, title = "log10(η)")
# ax4 = Axis(fig[2,3], aspect = ar, title = "log10(εII)")
ax4 = Axis(fig[2,3], aspect = ar, title = "log10(η)")
# Plot temperature
h1 = heatmap!(ax1, xvi[1].*1e-3, xvi[2].*1e-3, Array(thermal.T[2:end-1,:]) , colormap=:batlow)
# Plot particles phase
# h2 = scatter!(ax2, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), markersize = 1)
h2 = heatmap!(ax2, xvi[1].*1e-3, xvi[2].*1e-3, Array(stokes.V.Vy) , colormap=:batlow)
h2 = scatter!(ax2, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), markersize = 1)
# h2 = heatmap!(ax2, xvi[1].*1e-3, xvi[2].*1e-3, Array(stokes.V.Vy) , colormap=:batlow)
# Plot 2nd invariant of strain rate
# h3 = heatmap!(ax3, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.ε_pl.II)) , colormap=:batlow)
h3 = heatmap!(ax3, xci[1].*1e-3, xci[2].*1e-3, Array(stokes.τ.II)./1e6 , colormap=:batlow)
# Plot effective viscosity
h4 = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.ε_pl.II)) , colormap=:lipari)
# h4 = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.viscosity.η_vep)), colorrange=log10.(viscosity_cutoff), colormap=:batlow)
# h4 = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.ε.II)) , colormap=:lipari)
h4 = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.viscosity.η_vep)), colorrange=log10.(viscosity_cutoff), colormap=:batlow)
hidexdecorations!(ax1)
hidexdecorations!(ax2)
hidexdecorations!(ax3)
Expand All @@ -403,27 +413,27 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
fig
save(joinpath(figdir, "$(it).png"), fig)

## Plot initial T and P profile
fig = 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="Pressure")
scatter!(
ax1,
Array(thermal.T[2:(end - 1), :][:]),
Yv,
)
lines!(
ax2,
Array(stokes.P[:]),
Y,
)
hideydecorations!(ax2)
save(joinpath(figdir, "thermal_profile_$it.png"), fig)
fig
end
# ## Plot initial T and P profile
# fig = 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="Pressure")
# scatter!(
# ax1,
# Array(thermal.T[2:(end - 1), :][:]),
# Yv,
# )
# lines!(
# ax2,
# Array(stokes.P[:]),
# Y,
# )
# hideydecorations!(ax2)
# save(joinpath(figdir, "thermal_profile_$it.png"), fig)
# fig
# end
end
# ------------------------------

Expand All @@ -434,21 +444,22 @@ end

## END OF MAIN SCRIPT ----------------------------------------------------------------
do_vtk = true # set to true to generate VTK files for ParaView
# figdir = "Caldera2D_noPguess"
figdir = "Caldera2D"
n = 256
n = 128 * 2
nx, ny = n, n >>> 1
li, origin, phases_GMG, T_GMG = setup2D(
nx+1, ny+1;
sticky_air = 2.5,
flat = true,
sticky_air = 4,
flat = false,
chimney = false,
chamber_T = 1e3,
chamber_depth = 7e0,
chamber_radius = 1e0/2,
aspect_x = 10,
chamber_radius = 0.5,
aspect_x = 6,
)

heatmap(phases_GMG)
# heatmap(phases_GMG)
# heatmap(T_GMG)

igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
Expand All @@ -457,4 +468,4 @@ else
igg
end

main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
17 changes: 12 additions & 5 deletions volcano2D/Caldera_rheology.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,11 @@ function init_rheologies(; linear=false)
el = ConstantElasticity(; G = 25e9, ν = 0.45)
β = 1 / el.Kb.val
Cp = 1200.0

#dislocation laws
disl_top = SetDislocationCreep(Dislocation.dry_olivine_Karato_2003)
# diffusion laws
disl_bot = SetDislocationCreep(Dislocation.wet_quartzite_Hirth_2001)

# Define rheolgy struct
rheology = (
Expand All @@ -61,7 +66,8 @@ function init_rheologies(; linear=false)
Density = PT_Density(; ρ0=2.7e3, T0=273.15, β=β),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e23), el, pl)),
CompositeRheology = CompositeRheology( (disl_top, el, pl)),
# CompositeRheology = CompositeRheology( (LinearViscous(; η=1e23), el, pl)),
Gravity = ConstantGravity(; g=9.81),
),
# Name = "Lower crust",
Expand All @@ -70,7 +76,8 @@ function init_rheologies(; linear=false)
Density = PT_Density(; ρ0=2.7e3, T0=273.15, β=β),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e22), el, pl)),
CompositeRheology = CompositeRheology( (disl_bot, el, pl)),
# CompositeRheology = CompositeRheology( (LinearViscous(; η=1e22), el, pl)),
Gravity = ConstantGravity(; g=9.81),
),

Expand All @@ -82,7 +89,7 @@ function init_rheologies(; linear=false)
# HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity()),
HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=350e3J/kg),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e18), )),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e18), el)),
),
# Name = "magma chamber - hot anomaly",
SetMaterialParams(;
Expand All @@ -92,12 +99,12 @@ function init_rheologies(; linear=false)
# HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity()),
HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=350e3J/kg),
LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e18), )),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e18), el, )),
),
# Name = "StickyAir",
SetMaterialParams(;
Phase = 5,
Density = ConstantDensity(; ρ=0e0),
Density = ConstantDensity(; ρ=1e0),
HeatCapacity = ConstantHeatCapacity(; Cp = Cp),
Conductivity = ConstantConductivity(; k = 2.5),
CompositeRheology = CompositeRheology( (LinearViscous(; η=1e22), el, pl)),
Expand Down
4 changes: 2 additions & 2 deletions volcano2D/Caldera_setup.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,9 @@ function setup2D(
if chimney
add_cylinder!(Phases, Temp, Grid;
base = (mean(Grid.x.val), 0, -chamber_depth),
cap = (mean(Grid.x.val), 0, 0e0),
cap = (mean(Grid.x.val), 0, 3e0),
radius = 0.05,
phase = ConstantPhase(2),
phase = ConstantPhase(3),
# T = LinearTemp(Ttop=20, Tbot=1000),
# T = ConstantTemp(T=800),
T = ConstantTemp(T=chamber_T),
Expand Down

0 comments on commit c4abfa7

Please sign in to comment.