Merge branch 'main' into pa-KraflaSill

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.19.0
+        uses: crate-ci/typos@v1.20.1

.github/workflows/ci.yml

@@ -25,7 +25,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
@@ -42,7 +42,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4.0.1
+      - uses: codecov/codecov-action@v4.1.1
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           fail_ci_if_error: false  # or true if you want CI to fail when Codecov fails

.github/workflows/format_check.yml

@@ -20,7 +20,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

Project.toml

@@ -26,15 +26,15 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 AMDGPU = "0.6, 0.7, 0.8"
-Adapt = "3.7.2"
+Adapt = "3"
 CUDA = "4.4.1, 5"
-CellArrays = "0.1, 0.2"
+CellArrays = "0.2"
 GeoParams = "0.5"
 HDF5 = "0.17.1"
-ImplicitGlobalGrid = "0.14.0, 0.15.0"
+ImplicitGlobalGrid = "0.15.0"
 MPI = "0.20"
 MuladdMacro = "0.2"
-ParallelStencil = "0.9, 0.10, 0.11"
+ParallelStencil = "0.12"
 Reexport = "1.2.2"
 StaticArrays = "1"
 Statistics = "1"

_typos.toml

@@ -1,3 +1,4 @@
 [default.extend-words]
 Ths = "Ths"
 oly = "oly"
+iy = "iy"

miniapps/benchmarks/Blankenbach2D/Benchmark2D_WENO5.jl → miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl

@@ -115,7 +115,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     thermal         = ThermalArrays(ni)
     thermal_bc      = TemperatureBoundaryConditions(;
         no_flux     = (left = true, right = true, top = false, bot = false),
-        periodicity = (left = false, right = false, top = false, bot = false),
     )
     # initialize thermal profile
     @parallel (@idx size(thermal.T)) init_T!(thermal.T, xvi[2])
@@ -153,7 +152,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     # Boundary conditions -------------------------------
     flow_bcs         = FlowBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(stokes.V.Vx, stokes.V.Vy)

miniapps/benchmarks/Blankenbach2D/Benchmark2D_sgd.jl → miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl

@@ -121,7 +121,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     thermal         = ThermalArrays(ni)
     thermal_bc      = TemperatureBoundaryConditions(;
         no_flux     = (left = true, right = true, top = false, bot = false),
-        periodicity = (left = false, right = false, top = false, bot = false),
     )
     # initialize thermal profile
     @parallel (@idx size(thermal.T)) init_T!(thermal.T, xvi[2])
@@ -159,7 +158,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     # Boundary conditions -------------------------------
     flow_bcs         = FlowBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(stokes.V.Vx, stokes.V.Vy)

miniapps/benchmarks/Blankenbach2D/Benchmark2D_sgd_scaled.jl → miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl

@@ -116,7 +116,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     thermal         = ThermalArrays(ni)
     thermal_bc      = TemperatureBoundaryConditions(;
         no_flux     = (left = true, right = true, top = false, bot = false),
-        periodicity = (left = false, right = false, top = false, bot = false),
     )
     # initialize thermal profile
     @parallel (@idx size(thermal.T)) init_T!(thermal.T, xvi[2])
@@ -152,7 +151,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, nit = 1e1, figdir="figs2D", save_vtk
     # Boundary conditions -------------------------------
     flow_bcs         = FlowBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(stokes.V.Vx, stokes.V.Vy)

miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl

@@ -37,25 +37,25 @@ end
 
 function init_phases!(phases, particles)
     ni = size(phases)
-    
+
     @parallel_indices (i, j) function init_phases!(phases, px, py, index)
         r=100e3
         f(x, A, λ) = A * sin(π*x/λ)
-        
+
         @inbounds for ip in JustRelax.cellaxes(phases)
             # quick escape
             JustRelax.@cell(index[ip, i, j]) == 0 && continue
 
             x = JustRelax.@cell px[ip, i, j]
-            depth = -(JustRelax.@cell py[ip, i, j]) 
+            depth = -(JustRelax.@cell py[ip, i, j])
             @cell phases[ip, i, j] = 2.0
-            
+
             if 0e0 ≤ depth ≤ 100e3
                 @cell phases[ip, i, j] = 1.0
 
-            else 
+            else
                 @cell phases[ip, i, j] = 2.0
-                
+
                 if ((x - 250e3)^2 + (depth - 250e3)^2 ≤ r^2)
                     @cell phases[ip, i, j] = 3.0
                 end
@@ -121,7 +121,7 @@ function main(igg, nx, ny)
     pT, pPhases      = init_cell_arrays(particles, Val(2))
     particle_args    = (pT, pPhases)
 
-    # Elliptical temperature anomaly 
+    # Elliptical temperature anomaly
     init_phases!(pPhases, particles)
     phase_ratios  = PhaseRatio(ni, length(rheology))
     @parallel (@idx ni) phase_ratios_center(phase_ratios.center, pPhases)
@@ -136,7 +136,7 @@ function main(igg, nx, ny)
     # TEMPERATURE PROFILE --------------------------------
     thermal          = ThermalArrays(ni)
     # ----------------------------------------------------
-   
+
     # Buoyancy forces & rheology
     ρg               = @zeros(ni...), @zeros(ni...)
     η                = @ones(ni...)
@@ -149,7 +149,7 @@ function main(igg, nx, ny)
     η_vep            = copy(η)
 
     # Boundary conditions
-    flow_bcs         = FlowBoundaryConditions(; 
+    flow_bcs         = FlowBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
     )
 
@@ -181,7 +181,7 @@ function main(igg, nx, ny)
 
         # Stokes solver ----------------
         args = (; T = thermal.Tc, P = stokes.P,  dt=Inf)
-        
+
         @parallel (JustRelax.@idx ni) compute_ρg!(ρg[2], phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P))
         @parallel init_P!(stokes.P, ρg[2], xci[2])
         @parallel (@idx ni) compute_viscosity!(
@@ -212,13 +212,13 @@ function main(igg, nx, ny)
         # advect particles in space
         advection_RK!(particles, @velocity(stokes), grid_vx, grid_vy, dt, 2 / 3)
         # advect particles in memory
-        move_particles!(particles, xvi, particle_args)        
+        move_particles!(particles, xvi, particle_args)
         # check if we need to inject particles
         inject = check_injection(particles)
         inject && inject_particles_phase!(particles, pPhases, (), (), xvi)
         # update phase ratios
         @parallel (@idx ni) phase_ratios_center(phase_ratios.center, pPhases)
-        
+
         @show it += 1
         t        += dt
 
@@ -230,7 +230,7 @@ function main(igg, nx, ny)
             heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(η)), colormap = :grayC)
             arrows!(
                 ax,
-                xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))..., 
+                xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))...,
                 lengthscale = 25 / max(maximum(Vx_v),  maximum(Vy_v)),
                 color = :red,
             )
@@ -253,4 +253,4 @@ else
     igg
 end
 
-main(igg, nx, ny)
+main(igg, nx, ny)

miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl

@@ -30,7 +30,7 @@ macro all_j(A)
 end
 
 @parallel_indices (i, j) function init_P!(P, ρg, z)
-    P[i, j] = sum(abs(ρg[i, jj] * z[jj]) * z[jj] < 0.0 for jj in j:size(P, 2))
+    P[i, j] = sum(abs(ρg[i, jj] * z[jj]) for jj in j:size(P, 2))
     return nothing
 end
 
@@ -68,9 +68,9 @@ end
 function RT_2D(igg, nx, ny)
 
     # Physical domain ------------------------------------
-    thick_air    = 100e3              # thickness of sticky air layer
+    thick_air    = 100e3             # thickness of sticky air layer
     ly           = 500e3 + thick_air # domain length in y
-    lx           = 500e3           # domain length in x
+    lx           = 500e3             # domain length in x
     ni           = nx, ny            # number of cells
     li           = lx, ly            # domain length in x- and y-
     di           = @. li / ni        # grid step in x- and -y
@@ -84,8 +84,8 @@ function RT_2D(igg, nx, ny)
         # Name              = "Air",
         SetMaterialParams(;
             Phase             = 1,
-            Density           = ConstantDensity(; ρ=1e1),
-            CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
+            Density           = ConstantDensity(; ρ=1e0),
+            CompositeRheology = CompositeRheology((LinearViscous(; η=1e16),)),
             Gravity           = ConstantGravity(; g=9.81),
         ),
         # Name              = "Crust",
@@ -120,14 +120,13 @@ function RT_2D(igg, nx, ny)
     A             = 5e3    # Amplitude of the anomaly
     init_phases!(pPhases, particles, A)
     phase_ratios  = PhaseRatio(ni, length(rheology))
-    # @parallel (@idx ni) phase_ratios_center(phase_ratios.center, pPhases)
     @parallel (@idx ni) phase_ratios_center(phase_ratios.center, particles.coords, xci, di, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
     stokes           = StokesArrays(ni, ViscoElastic)
-    pt_stokes        = PTStokesCoeffs(li, di; ϵ=1e-5,  CFL = 0.95 / √2.1)
+    pt_stokes        = PTStokesCoeffs(li, di; ϵ=1e-8,  CFL = 0.95 / √2.1)
     # ----------------------------------------------------
 
     # TEMPERATURE PROFILE --------------------------------
@@ -136,19 +135,20 @@ function RT_2D(igg, nx, ny)
 
     # Buoyancy forces & rheology
     ρg               = @zeros(ni...), @zeros(ni...)
-    η                = @zeros(ni...)
+    η                = @ones(ni...)
     args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
     @parallel (JustRelax.@idx ni) compute_ρg!(ρg[2], phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P))
     @parallel init_P!(stokes.P, ρg[2], xci[2])
     @parallel (@idx ni) compute_viscosity!(
-        η, 1.0, phase_ratios.center, @strain(stokes)..., args, rheology, (1e19, Inf)
+        η, 0.0, phase_ratios.center, @strain(stokes)..., args, rheology, (1e19, 1e24)
     )
     η_vep            = copy(η)
 
     # Boundary conditions
     flow_bcs         = FlowBoundaryConditions(; 
-        free_slip    = (left = true, right=true, top=true, bot=false),
-        no_slip      = (left = false, right=false, top=false, bot=true),
+        free_slip    = (left =  true, right =  true, top =  true, bot = false),
+        no_slip      = (left = false, right = false, top = false, bot =  true),
+        free_surface = true,
     )
 
     # Plot initial T and η profiles
@@ -159,7 +159,6 @@ function RT_2D(igg, nx, ny)
         ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
         scatter!(ax1, Array(ρg[2][:]./9.81), Y./1e3)
         scatter!(ax2, Array(log10.(η[:])), Y./1e3)
-        # scatter!(ax2, Array(stokes.P[:]), Y./1e3)
         ylims!(ax1, minimum(xvi[2])./1e3, 0)
         ylims!(ax2, minimum(xvi[2])./1e3, 0)
         hideydecorations!(ax2)
@@ -174,7 +173,8 @@ function RT_2D(igg, nx, ny)
 
     # Time loop
     t, it = 0.0, 0
-    dt = dt_max = 50e3 * (3600 * 24 * 365.25)
+    dt = 1e3 * (3600 * 24 * 365.25)
+    dt_max = 50e3 * (3600 * 24 * 365.25)
     while it < 500 # run only for 5 Myrs
 
         # Stokes solver ----------------
@@ -184,7 +184,6 @@ function RT_2D(igg, nx, ny)
         @parallel (@idx ni) compute_viscosity!(
             η, 1.0, phase_ratios.center, @strain(stokes)..., args, rheology, (-Inf, Inf)
         )
-        # η .= mean(η)
         solve!(
             stokes,
             pt_stokes,
@@ -198,22 +197,22 @@ function RT_2D(igg, nx, ny)
             args,
             dt,
             igg;
-            iterMax = 75e3,
-            iterMin = 15e3,
-            viscosity_relaxation = 1e-1,
-            nout    = 5e3,
-            viscosity_cutoff=(-Inf, Inf)
+            iterMax              = 150e3,
+            iterMin              =   5e3,
+            viscosity_relaxation =  1e-2,
+            nout                 =   5e3,
+            free_surface         =  true,
+            viscosity_cutoff     = (-Inf, Inf)
         )
-        # dt = if it ≤ 10
-        #     min(compute_dt(stokes, di) / 2, 1e3 * (3600 * 24 * 365.25))
-        # elseif 10 < it ≤ 20
-        #     min(compute_dt(stokes, di) / 2, 10e3 * (3600 * 24 * 365.25))
-        # elseif 20 < it ≤ 30
-        #     min(compute_dt(stokes, di) / 2, 25e3 * (3600 * 24 * 365.25))
-        # else
-        #     min(compute_dt(stokes, di) / 2, dt_max)
-        # end
-        dt = min(compute_dt(stokes, di) / 2, dt_max)
+        dt = if it ≤ 10
+            min(compute_dt(stokes, di),  1e3 * (3600 * 24 * 365.25))
+        elseif 10 < it ≤ 20
+            min(compute_dt(stokes, di), 10e3 * (3600 * 24 * 365.25))
+        elseif 20 < it ≤ 30
+            min(compute_dt(stokes, di), 25e3 * (3600 * 24 * 365.25))
+        else
+            min(compute_dt(stokes, di), dt_max)
+        end
         # ------------------------------
 
         # Advection --------------------
@@ -241,17 +240,15 @@ function RT_2D(igg, nx, ny)
             pyv      = ppy.data[:]./1e3
             clr      = pPhases.data[:]
 
-            fig = Figure(resolution = (900, 900), title = "t = $t")
+            fig = Figure(size = (900, 900), title = "t = $t")
             ax  = Axis(fig[1,1], aspect = 1, title = " t=$(t/(1e3 * 3600 * 24 *365.25)) Kyrs")
-            # heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(η)), colormap = :grayC)
             scatter!(
                 ax, 
                 pxv, pyv, 
                 color=clr, 
                 colormap = :lajolla,
                 markersize = 3
             )
-            # scatter!(ax, pxv, pyv, color=clr, colormap = :grayC)
             arrows!(
                 ax,
                 xvi[1][1:nt:end-1]./1e3, xvi[2][1:nt:end-1]./1e3, Array.((Vx_v[1:nt:end-1, 1:nt:end-1], Vy_v[1:nt:end-1, 1:nt:end-1]))..., 

miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl

@@ -102,7 +102,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     thermal          = ThermalArrays(ni)
     thermal_bc       = TemperatureBoundaryConditions(;
         no_flux      = (left = true, right = true, top = false, bot = false),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
 
     # Initialize constant temperature
@@ -134,7 +133,6 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     # Boundary conditions
     flow_bcs         = FlowBoundaryConditions(;
         free_slip    = (left = true, right=true, top=true, bot=true),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
     ## Compression and not extension - fix this
     εbg              = 5e-14

miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl

@@ -146,7 +146,6 @@ function sinking_block2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_per
     # Boundary conditions
     flow_bcs = FlowBoundaryConditions(;
         free_slip    = (left =  true, right =  true, top =  true, bot =  true),
-        periodicity  = (left = false, right = false, top = false, bot = false),
     )
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(stokes.V.Vx, stokes.V.Vy)

miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl

@@ -211,9 +211,6 @@ function burstedde(; nx=16, ny=16, nz=16, init_MPI=true, finalize_MPI=false)
     flow_bcs = FlowBoundaryConditions(;
         free_slip   = (left=false, right=false, top=false, bot=false, back=false, front=false),
         no_slip     = (left=false, right=false, top=false, bot=false, back=false, front=false),
-        periodicity = (
-            left=false, right=false, top=false, bot=false, back=false, front=false
-        ),
     )
     # impose analytical velociity at the boundaries of the domain
     velocity!(stokes, xci, xvi, di)