up miniapp caldera

PTsolvers · Dec 16, 2024 · cc128a0 · cc128a0
1 parent 3f27143
commit cc128a0
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 31 deletions.
diff --git 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,24 +583,24 @@ 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,
     dimensions     = (30e0, 20e0), # extent in x and y in km
     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
@@ -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
@@ -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