Merge pull request #285 from PTsolvers/pa-vtk_chain

New `save_maker_chain` function

.github/workflows/Documenter.yml

@@ -33,7 +33,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.10'
+          version: '1.11'
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@latest
       - name: Install dependencies

Project.toml

@@ -34,11 +34,11 @@ AMDGPU = "0.8, 0.9, 1"
 Adapt = "4"
 CUDA = "5"
 CellArrays = "0.2, 0.3"
-GeoParams = "0.6.7"
+GeoParams = "0.6.8"
 HDF5 = "0.17.1"
-ImplicitGlobalGrid = "0.15"
+ImplicitGlobalGrid = "0.15, 0.16"
 JLD2 = "0.4, 0.5"
-JustPIC = "0.5.3"
+JustPIC = "0.5"
 MPI = "0.20"
 MuladdMacro = "0.2"
 ParallelStencil = "0.13.6, 0.14"

docs/src/man/Blankenbach.md

@@ -153,13 +153,13 @@ and we define a rectangular thermal anomaly at $x \in [0, 0.05]$, $y \in [\frac{
 ```julia
 function rectangular_perturbation!(T, xc, yc, r, xvi)
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
-            T[i, j] += .2
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
+            T[i+1, j] += .2
         end
         return nothing
     end
-    ni = size(T)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, xvi...)
     return nothing
 end
 

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl

@@ -41,13 +41,13 @@ end
 # Thermal rectangular perturbation
 function rectangular_perturbation!(T, xc, yc, r, xvi)
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
-            T[i, j] += 20.0
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
+            T[i+1, j] += 20.0
         end
         return nothing
     end
-    ni = size(T)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
+    nx,ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, xvi...)
     return nothing
 end
 ## END OF HELPER FUNCTION ------------------------------------------------------------
@@ -328,8 +328,8 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
             fig
 
             fig2 = Figure(size = (900, 1200), title = "Time Series")
-            ax21 = Axis(fig2[1,1], aspect = 3, title = "V_{RMS}")
-            ax22 = Axis(fig2[2,1], aspect = 3, title = "Nu_{top}")
+            ax21 = Axis(fig2[1,1], aspect = 3, title = L"V_{RMS}")
+            ax22 = Axis(fig2[2,1], aspect = 3, title = L"Nu_{top}")
             l1 = lines!(ax21,trms./(1e6*(365.25*24*60*60)),Urms)
             l2 = lines!(ax22,trms./(1e6*(365.25*24*60*60)),Nu_top)
             save(joinpath(figdir, "Time_Series_V_Nu.png"), fig2)

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl

@@ -42,13 +42,13 @@ end
 # Thermal rectangular perturbation
 function rectangular_perturbation!(T, xc, yc, r, xvi)
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
-            T[i, j] += 20.0
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
+            T[i+1, j] += 20.0
         end
         return nothing
     end
-    ni = size(T)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, xvi...)
     return nothing
 end
 ## END OF HELPER FUNCTION ------------------------------------------------------------
@@ -355,8 +355,8 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
             fig
 
             fig2 = Figure(size = (900, 1200), title = "Time Series")
-            ax21 = Axis(fig2[1,1], aspect = 3, title = "V_{RMS}")
-            ax22 = Axis(fig2[2,1], aspect = 3, title = "Nu_{top}")
+            ax21 = Axis(fig2[1,1], aspect = 3, title = L"V_{RMS}")
+            ax22 = Axis(fig2[2,1], aspect = 3, title = L"Nu_{top}")
             l1 = lines!(ax21,trms./(1e6*(365.25*24*60*60)),Urms)
             l2 = lines!(ax22,trms./(1e6*(365.25*24*60*60)),Nu_top)
             save(joinpath(figdir, "Time_Series_V_Nu.png"), fig2)

miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl

@@ -35,13 +35,13 @@ end
 # Thermal rectangular perturbation
 function rectangular_perturbation!(T, xc, yc, r, xvi)
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
-            T[i, j] += .2
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc)^2 ≤ r^2)
+            T[i+1, j] += .2
         end
         return nothing
     end
-    ni = size(T)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, xvi...)
     return nothing
 end
 ## END OF HELPER FUNCTION ------------------------------------------------------------
@@ -347,8 +347,8 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
             fig
 
             fig2 = Figure(size = (900, 1200), title = "Time Series")
-            ax21 = Axis(fig2[1,1], aspect = 3, title = "(V_{RMS})")
-            ax22 = Axis(fig2[2,1], aspect = 3, title = "(Nu_{top})")
+            ax21 = Axis(fig2[1,1], aspect = 3, title = L"V_{RMS}")
+            ax22 = Axis(fig2[2,1], aspect = 3, title = L"Nu_{top}")
             l1 = lines!(ax21,trms,(Urms))
             l2 = lines!(ax22,trms,(Nu_top))
             save(joinpath(figdir, "Time_Series_V_Nu.png"), fig2)

miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl

@@ -103,11 +103,11 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
 
     # Rheology
     args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
-    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
 
     # PT coefficients for thermal diffusion
     pt_thermal       = PTThermalCoeffs(
-        rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=5e-2 / √3
+        backend_JR, rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=0.95/ √3.1
     )
 
     # Boundary conditions
@@ -117,9 +117,9 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
     )
     ## Compression and not extension - fix this
     εbg              = 5e-14
-    stokes.V.Vx .= PTArray([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2, _ in 1:nz+2])
-    stokes.V.Vy .= PTArray([ -(y - ly/2) * εbg for _ in 1:nx+2, y in xvi[2], _ in 1:nz+2])
-    stokes.V.Vz .= PTArray([  (lz - abs(z)) * εbg for _ in 1:nx+2, _ in 1:ny+2, z in xvi[3]])
+    stokes.V.Vx .= PTArray(backend_JR)([ -(x - lx/2) * εbg for x in xvi[1], _ in 1:ny+2, _ in 1:nz+2])
+    stokes.V.Vy .= PTArray(backend_JR)([ -(y - ly/2) * εbg for _ in 1:nx+2, y in xvi[2], _ in 1:nz+2])
+    stokes.V.Vz .= PTArray(backend_JR)([  (lz - abs(z)) * εbg for _ in 1:nx+2, _ in 1:ny+2, z in xvi[3]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(@velocity(stokes)...)
 
@@ -160,7 +160,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
     while it < 10
             # Update buoyancy and viscosity -
             args = (; T = thermal.Tc, P = stokes.P,  dt = Inf)
-            compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+            compute_viscosity!(stokes, phase_ratios, args, rheology, (-Inf, Inf))
             compute_ρg!(ρg[3], phase_ratios, rheology, args)
             # ------------------------------
 
@@ -296,10 +296,10 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
             end
             # ------------------------------
 
-      end
+        end
 
-      return nothing
-  end
+    return nothing
+end
 
 figdir   = "3D_Benchmark_Duretz_etal_2014"
 do_vtk = false # set to true to generate VTK files for ParaView

miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D.jl

@@ -26,13 +26,13 @@ end
 function elliptical_perturbation!(T, δT, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
-        @inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j]  += δT
+        if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j]  += δT
         end
         return nothing
     end
-
-    @parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
 end
 
 function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, ρ0=3.3e3, Cp0=1.2e3, K0=3.0)

miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_MPI.jl

@@ -24,13 +24,13 @@ end
 function elliptical_perturbation!(T, δT, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
-        @inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j]  += δT
+        if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j]  += δT
         end
         return nothing
     end
-
-    @parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
 end
 
 function diffusion_2D(figdir;

miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_inner_BCs.jl

@@ -24,14 +24,14 @@ end
 function elliptical_perturbation!(T, mask, Ω_T, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _elliptical_perturbation!(T, mask, Ω_T, xc, yc, r, x, y)
-        @inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j]   += Ω_T
+        if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j]   += Ω_T
             mask[i, j] = 1
         end
         return nothing
     end
-
-    @parallel _elliptical_perturbation!(T, mask, Ω_T, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _elliptical_perturbation!(T, mask, Ω_T, xc, yc, r, xvi...)
 end
 
 function init_phases!(phases, particles, xc, yc, r)
@@ -170,4 +170,3 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0)
 end
 
 thermal =  diffusion_2D()
-

miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl

@@ -25,13 +25,13 @@ end
 function elliptical_perturbation!(T, δT, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
-        @inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j]  += δT
+        if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j]  += δT
         end
         return nothing
     end
-
-    @parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
+    nx,ny = size(T)
+    @parallel (1:nx-2, 1:ny) _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
 end
 
 function init_phases!(phases, particles, xc, yc, r)

miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase_MPI.jl

@@ -25,13 +25,13 @@ end
 function elliptical_perturbation!(T, δT, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _elliptical_perturbation!(T, δT, xc, yc, r, x, y)
-        @inbounds if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j]  += δT
+        if (((x[i]-xc ))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j]  += δT
         end
         return nothing
     end
-
-    @parallel _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _elliptical_perturbation!(T, δT, xc, yc, r, xvi...)
 end
 
 function init_phases!(phases, particles, xc, yc, r)

miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl

@@ -103,7 +103,7 @@ function circular_perturbation!(T, δT, xc_anomaly, yc_anomaly, r_anomaly, xvi,
         T, δT, xc_anomaly, yc_anomaly, r_anomaly, x, y, sticky_air
     )
         depth = -y[j] - sticky_air
-        @inbounds if ((x[i] - xc_anomaly)^2 + (depth[j] + yc_anomaly)^2 ≤ r_anomaly^2)
+        if ((x[i] - xc_anomaly)^2 + (depth[j] + yc_anomaly)^2 ≤ r_anomaly^2)
             # T[i + 1, j] *= δT / 100 + 1
             T[i + 1, j] = δT
         end
@@ -286,10 +286,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
     flow_bcs!(stokes, flow_bcs)
     update_halo!(@velocity(stokes)...)
 
-    η = @ones(ni...) # initialise viscosity
-
-    compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
-    η_vep = copy(η)
+    compute_viscosity!(stokes, phase_ratios, args, rheology, 0, cutoff_visc)
 
     # Buoyancy force
     ρg = @zeros(ni...), @zeros(ni...) # ρg[1] is the buoyancy force in the x direction, ρg[2] is the buoyancy force in the y direction
@@ -360,7 +357,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
         # Update buoyancy and viscosity -
         args = (; T=thermal.Tc, P=stokes.P, dt=Inf, ΔTc=thermal.ΔTc)
         compute_ρg!(ρg[end], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P))
-        compute_viscosity!(stokes, phase_ratios, args, rheology, cutoff_visc)
+        compute_viscosity!(stokes, phase_ratios, args, rheology, 0, cutoff_visc)
 
         # Stokes solver -----------------
         solve!(
@@ -691,8 +688,8 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
 end
 
 figdir = "Thermal_stresses_around_cooling_magma"
-do_vtk = false # set to true to generate VTK files for ParaView
-n      = 64
+do_vtk = true # set to true to generate VTK files for ParaView
+n      = 128
 ar     = 2
 nx     = n * ar - 2
 ny     = n - 2

miniapps/convection/GlobalConvection2D_WENO5.jl

@@ -55,26 +55,26 @@ end
 function circular_perturbation!(T, δT, xc, yc, r, xvi)
 
     @parallel_indices (i, j) function _circular_perturbation!(T, δT, xc, yc, r, x, y)
-        @inbounds if (((x[i] - xc))^2 + ((y[j] - yc))^2) ≤ r^2
-            T[i, j] *= δT / 100 + 1
+        if (((x[i] - xc))^2 + ((y[j] - yc))^2) ≤ r^2
+            T[i+1, j] *= δT / 100 + 1
         end
         return nothing
     end
-
-    @parallel _circular_perturbation!(T, δT, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _circular_perturbation!(T, δT, xc, yc, r, xvi...)
 end
 
 function random_perturbation!(T, δT, xbox, ybox, xvi)
 
     @parallel_indices (i, j) function _random_perturbation!(T, δT, xbox, ybox, x, y)
-        @inbounds if (xbox[1] ≤ x[i] ≤ xbox[2]) && (abs(ybox[1]) ≤ abs(y[j]) ≤ abs(ybox[2]))
+        if (xbox[1] ≤ x[i] ≤ xbox[2]) && (abs(ybox[1]) ≤ abs(y[j]) ≤ abs(ybox[2]))
             δTi = δT * (rand() - 0.5) # random perturbation within ±δT [%]
-            T[i, j] *= δTi / 100 + 1
+            T[i+1, j] *= δTi / 100 + 1
         end
         return nothing
     end
-
-    @parallel (@idx size(T)) _random_perturbation!(T, δT, xbox, ybox, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _random_perturbation!(T, δT, xbox, ybox, xvi...)
 end
 
 # --------------------------------------------------------------------------------

miniapps/convection/Particles2D/Layered_convection2D.jl

@@ -74,16 +74,16 @@ end
 function rectangular_perturbation!(T, xc, yc, r, xvi, thick_air)
 
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc - thick_air)^2 ≤ r^2)
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc - thick_air)^2 ≤ r^2)
             depth       = -y[j] - thick_air
             dTdZ        = (2047 - 2017) / 50e3
             offset      = 2017
             T[i + 1, j] = (depth - 585e3) * dTdZ + offset
         end
         return nothing
     end
-    ni = length.(xvi)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, xvi...)
 
     return nothing
 end

miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl

@@ -73,7 +73,7 @@ end
 function rectangular_perturbation!(T, xc, yc, r, xvi, thick_air, CharDim)
 
     @parallel_indices (i, j) function _rectangular_perturbation!(T, xc, yc, r, CharDim, x, y)
-        @inbounds if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc - thick_air)^2 ≤ r^2)
+        if ((x[i]-xc)^2 ≤ r^2) && ((y[j] - yc - thick_air)^2 ≤ r^2)
             depth       = -y[j] - thick_air
             dTdZ        = nondimensionalize((2047 - 2017)K / 50km, CharDim)
             offset      = nondimensionalize(2017e0K, CharDim)
@@ -82,8 +82,8 @@ function rectangular_perturbation!(T, xc, yc, r, xvi, thick_air, CharDim)
         return nothing
     end
 
-    ni = length.(xvi)
-    @parallel (@idx ni) _rectangular_perturbation!(T, xc, yc, r, CharDim, xvi...)
+    nx, ny = size(T)
+    @parallel (1:nx-2, 1:ny) _rectangular_perturbation!(T, xc, yc, r, CharDim, xvi...)
 
     return nothing
 end