From cc128a0f9e721cbac8cee8bc1debd423174ed1e2 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 16 Dec 2024 15:28:30 +0100 Subject: [PATCH] up miniapp caldera --- .../stokes2D/Volcano2D/Caldera2D.jl | 35 +++++-------------- .../stokes2D/Volcano2D/Caldera_rheology.jl | 7 ++-- .../stokes2D/Volcano2D/Caldera_setup.jl | 4 +-- 3 files changed, 15 insertions(+), 31 deletions(-) diff --git a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera2D.jl b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera2D.jl index d1d7b745..51b71033 100644 --- a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera2D.jl +++ b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera2D.jl @@ -124,7 +124,7 @@ function thermal_anomaly!(Temp, Ω_T, phase_ratios, T_chamber, T_air, conduit_ph end ## BEGIN OF MAIN SCRIPT -------------------------------------------------------------- -function main(li, origin, phases_GMG, T_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk =false) +function main(li, origin, phases_GMG, T_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk =false, extension = 1e-15 * 0) # Physical domain ------------------------------------ ni = nx, ny # number of cells @@ -249,7 +249,7 @@ function main(li, origin, phases_GMG, T_GMG, igg; nx=16, ny=16, figdir="figs2D", # flow_bcs!(stokes, flow_bcs) # apply boundary conditions # displacement2velocity!(stokes, dt) - εbg = 1e-15 * 1 + εbg = extension apply_pure_shear(@velocity(stokes)..., εbg, xvi, li...) flow_bcs!(stokes, flow_bcs) # apply boundary conditions @@ -583,12 +583,12 @@ end ## END OF MAIN SCRIPT ---------------------------------------------------------------- const plotting = true - +extension = 1e-15 * 1 do_vtk = true # set to true to generate VTK files for ParaView -# figdir = "Caldera2D_noPguess" -figdir = "$(today())_Conduit" +figdir = "Caldera2D" n = 64 nx, ny = n, n >>> 1 + li, origin, phases_GMG, T_GMG = setup2D( nx+1, ny+1; sticky_air = 4e0, @@ -596,11 +596,11 @@ li, origin, phases_GMG, T_GMG = setup2D( flat = false, chimney = true, volcano_size = (3e0, 5e0), - conduit_radius = 0.5e0, + conduit_radius = 2e-1, chamber_T = 900e0, chamber_depth = 7e0, - chamber_radius = 1e0, - aspect_x = 6, + chamber_radius = 1.25, + aspect_x = 2, ) igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid @@ -609,21 +609,4 @@ else igg end -main(li, origin, phases_GMG, T_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk); - -# function plot_particles(particles, pPhases, chain) -# p = particles.coords -# # pp = [argmax(p) for p in phase_ratios.center] #if you want to plot it in a heatmap rather than scatter -# ppx, ppy = p -# # pxv = ustrip.(dimensionalize(ppx.data[:], km, CharDim)) -# # pyv = ustrip.(dimensionalize(ppy.data[:], km, CharDim)) -# pxv = ppx.data[:] -# pyv = ppy.data[:] -# clr = pPhases.data[:] -# # clr = pϕ.data[:] -# idxv = particles.index.data[:] -# f,ax,h=scatter(Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]), colormap=:roma, markersize=1) - -# Colorbar(f[1,2], h) -# f -# end +main(li, origin, phases_GMG, T_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk, extension = extension); diff --git a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_rheology.jl b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_rheology.jl index 7bb80e26..6aa020ff 100644 --- a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_rheology.jl +++ b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_rheology.jl @@ -46,16 +46,17 @@ function init_rheologies(; linear=false) η_reg = 1e18 C = linear ? Inf : 10e6 ϕ = 15 - soft_C = NonLinearSoftening(; ξ₀=C, Δ = C / 2) # nonlinear softening law + soft_C = NonLinearSoftening(; ξ₀=C, Δ = C / 1e5) # nonlinear softening law # soft_C = NonLinearSoftening() # nonlinear softening law - pl = DruckerPrager_regularised(; C=C, ϕ=ϕ, η_vp=η_reg, Ψ=0.0, softening_C=soft_C) + pl = DruckerPrager_regularised(; C=C*MPa, ϕ=ϕ, η_vp=(η_reg)*Pas, Ψ=0.0, softening_C=soft_C) el = ConstantElasticity(; G = 25e9, ν = 0.45) β = 1 / el.Kb.val Cp = 1200.0 magma_visc = ViscosityPartialMelt_Costa_etal_2009(η=LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K,η0=1e3Pa*s)) #dislocation laws - disl_top = SetDislocationCreep(Dislocation.dry_olivine_Karato_2003) + disl_top = DislocationCreep(; A=1.67e-24, n=3.5, E=1.87e5, V=6e-6, r=0.0, R=8.3145) + # disl_top = SetDislocationCreep(Dislocation.dry_olivine_Karato_2003) # diffusion laws disl_bot = SetDislocationCreep(Dislocation.wet_quartzite_Hirth_2001) diff --git a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_setup.jl b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_setup.jl index d88f1087..d4676c40 100644 --- a/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_setup.jl +++ b/miniapps/benchmarks/stokes2D/Volcano2D/Caldera_setup.jl @@ -80,8 +80,8 @@ function setup2D( ph = Phases[:,1,:] T = Temp[:,1,:] .+ 273 - V = 4/3 * π * (chamber_radius*aspect_x) * chamber_radius * 1.0 - printstyled("Magma volume of the initial chamber: $(round(V; digits=3)) km³ \n"; bold=true, color=:red) + V = 4/3 * π * (chamber_radius*aspect_x) * chamber_radius * (chamber_radius*aspect_x) + printstyled("Magma volume of the initial chamber: $(round(V; digits=3)) km³ \n"; bold=true, color=:red, blink=true) # write_paraview(Grid, "Volcano2D") return li, origin, ph, T, Grid end