updates

Volcano3D.jl

@@ -235,9 +235,9 @@ function main3D(
             dt,
             igg;
             kwargs = (;
-                iterMax          = 50e3,#50e3,
+                iterMax          = 150e3,#50e3,
                 free_surface     = false,
-                nout             = 1e3,#5e3,
+                nout             = 2e3,#5e3,
                 viscosity_cutoff = cutoff_visc,
             )
         )
@@ -294,7 +294,7 @@ function main3D(
         t += dt
 
         #  # # Plotting -------------------------------------------------------
-        if it == 1 || rem(it, 5) == 0
+        if it == 1 || rem(it, 1) == 0
             (; η) = stokes.viscosity
             # checkpointing(figdir, stokes, thermal.T, η, t)
 
@@ -316,7 +316,7 @@ function main3D(
                         εyy = Array(ustrip.(dimensionalize(stokes.ε.yy, s^-1,CharDim))),
                         εzz = Array(ustrip.(dimensionalize(stokes.ε.zz, s^-1,CharDim))),
                         εII = Array(ustrip.(dimensionalize(stokes.ε.II, s^-1,CharDim))),
-                        η   = Array(ustrip.(dimensionalize(stokes.viscosity.η,Pa*s,CharDim))),
+                        η   = Array(ustrip.(dimensionalize(stokes.viscosity.η_vep,Pa*s,CharDim))),
                     )
                     velocity_v = (
                         Array(ustrip.(dimensionalize(Vx_v,cm/yr,CharDim))),

VolcanoRheology.jl

@@ -5,8 +5,8 @@ function init_rheology(CharDim; is_compressible = false)
     # plasticity setup
     do_DP   = true          # do_DP=false: Von Mises, do_DP=true: Drucker-Prager (friction angle)
     η_reg   = 1.0e16Pa*s    # regularisation "viscosity" for Drucker-Prager
-    Coh     = 10.0MPa       # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ)
-    ϕ       = 10.0 * do_DP  # friction angle
+    Coh     = 5.0MPa       # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ)
+    ϕ       = 5.0 * do_DP  # friction angle
     G0      = 6.0e11Pa      # elastic shear modulus
     G_magma = 6.0e11Pa      # elastic shear modulus perturbation
 
@@ -24,9 +24,9 @@ function init_rheology(CharDim; is_compressible = false)
         β_rock   = inv(get_Kb(el))
         β_magma  = inv(get_Kb(el_magma))
     end
-    creep_rock  = LinearViscous(; η=1e21 * Pa * s)
+    creep_rock  = LinearViscous(; η=1e24 * Pa * s)
     creep_magma = LinearViscous(; η=1e18 * Pa * s)
-    creep_air   = LinearViscous(; η=1e20 * Pa * s)
+    creep_air   = LinearViscous(; η=1e21 * Pa * s)
     g           = 9.81m/s^2
 
     rheology = (
@@ -76,7 +76,7 @@ function init_rheology(CharDim; is_compressible = false)
         #Name="Sticky Air"
         SetMaterialParams(;
             Phase             = 4,
-            Density           = ConstantDensity(ρ = 1000kg/m^3,),
+            Density           = ConstantDensity(ρ = 100kg/m^3,),
             # Density           = ConstantDensity(ρ=1kg/m^3,),
             HeatCapacity      = ConstantHeatCapacity(; Cp=1000J / kg / K),
             Conductivity      = ConstantConductivity(; k=15Watt / K / m),

VolcanoSetup.jl

@@ -1,9 +1,43 @@
-using GeophysicalModelGenerator
+using GLMakie
+using GeophysicalModelGenerator, GMT
+using Interpolations
+
+function Etna_topo(N)
+    Topo = import_topo([14.5, 15.5, 37.2, 38.2], file="@earth_relief_01m")
+    surf = Topo.depth.val[:,:,1]
+
+    l = maximum(abs.(Topo.depth.val))
+    # heatmap(Topo.depth.val[:,:,1], colormap=:oleron, colorrange= (-l,l))
+    # GLMakie.surface(surf, colormap=:oleron, colorrange= (-l,l))
+
+    ##
+    x = LinRange(14.5, 15.5, N) |> collect
+    y = LinRange(37.2, 38.2, N) |> collect
+
+    X = (
+        Topo.lon.val[:, 1, 1],
+        Topo.lat.val[1, :, 1],
+    )
+
+    itp = interpolate(X, surf, Gridded(Linear()))
+    surf_interp = [itp(x, y) for x in x, y in y]
+
+    # GLMakie.surface(surf_interp, colormap=:oleron, colorrange= (-l,l))
+
+    lat_dist  = extrema(X[1]) |> collect |> diff |> first |> abs
+    long_dist = extrema(X[2]) |> collect |> diff |> first |> abs
+    Ly        = lat_dist * 110.574
+    Lx        = long_dist * 111.320*cos(lat_dist)
+
+    return surf_interp, Lx, Ly
+end
 
 function volcano_setup(N)
+    topo_volcano, Lx, Ly = Etna_topo(N)
+
     nx = ny = nz = N
-    x = range(-50, 50, nx);
-    y = range(-50, 50, ny);
+    x = range(-Lx / 2, Lx / 2, nx);
+    y = range(-Ly / 2, Ly / 2, ny);
     z = range(-50, 10, nz);
     Grid = CartData(xyz_grid(x,y,z));
 
@@ -23,20 +57,25 @@ function volcano_setup(N)
         T = HalfspaceCoolingTemp(Age=20)
     )
 
-    add_volcano!(Phases, Temp, Grid;
-        volcanic_phase  = 1,
-        center     = (0, 0, 0),
-        height     = 4,
-        radius     = 10,
-        crater     = 0.0,
-        base       = 0.0,
-        background = nothing,
-        T = HalfspaceCoolingTemp(Age=20)
-    )
+    id = falses(nx, ny, nz)
+    for k in axes(topo_volcano,3), j in axes(topo_volcano, 2), i in axes(topo_volcano,1)
+        id[i, j, k] = Grid.z.val[i, j, k] > topo_volcano[i, j]
+    end
+
+    # add_volcano!(Phases, Temp, Grid;
+    #     volcanic_phase  = 1,
+    #     center     = (0, 0, 0),
+    #     height     = 4,
+    #     radius     = 10,
+    #     crater     = 0.0,
+    #     base       = 0.0,
+    #     background = nothing,
+    #     T = HalfspaceCoolingTemp(Age=20)
+    # )
 
     add_ellipsoid!(Phases, Temp, Grid; 
         cen    = (0, 0, -12.5), 
-        axes   = (10, 10, 5),
+        axes   = (5, 5, 2.5),
         phase  = ConstantPhase(3),
     )