PTsolvers · albert-de-montserrat · Oct 12, 2024 · Oct 3, 2024 · Oct 4, 2024 · Oct 4, 2024
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -15,6 +15,7 @@ jobs:
       matrix:
         version:
           - '1.10'
+          - '1.11'
           - 'pre'
           # - 'nightly'
         os:

diff --git a/Project.toml b/Project.toml
@@ -6,9 +6,7 @@ version = "0.3.2"
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CellArrays = "d35fcfd7-7af4-4c67-b1aa-d78070614af4"
-GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
 GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38"
-GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ImplicitGlobalGrid = "4d7a3746-15be-11ea-1130-334b0c4f5fa0"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"

diff --git a/docs/src/man/Blankenbach.md b/docs/src/man/Blankenbach.md
@@ -13,7 +13,7 @@ const backend_JR = CPUBackend
 For this benchmark we will use particles to track the advection of the material phases and their information. For this, we will use [JustPIC.jl](https://github.com/JuliaGeodynamics/JustPIC.jl)
 ```julia
 using JustPIC, JustPIC._2D
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 ```
 
 We will also use [ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl) to write some device-agnostic helper functions:

diff --git a/docs/src/man/subduction2D/subduction2D.md b/docs/src/man/subduction2D/subduction2D.md
@@ -56,22 +56,19 @@ For the rheology we will use the `rheology` object we created in the previous se
 
 ## Initialize particles
 ```julia
-nxcell              = 40 # initial number of particles per cell
-max_xcell           = 60 # maximum number of particles per cell
-min_xcell           = 20 # minimum number of particles per cell
-    particles = init_particles(
-        backend, nxcell, max_xcell, min_xcell, xvi...
-    )
-)
-subgrid_arrays      = SubgridDiffusionCellArrays(particles)
+nxcell          = 40 # initial number of particles per cell
+max_xcell       = 60 # maximum number of particles per cell
+min_xcell       = 20 # minimum number of particles per cell
+particles       = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...)
+subgrid_arrays  = SubgridDiffusionCellArrays(particles)
 # velocity staggered grids
-grid_vxi            = velocity_grids(xci, xvi, di)
+grid_vxi        = velocity_grids(xci, xvi, di)
 ```
 
 We will like to advect two fields, the temperature `pT` and the material phases of each particle `pPhases`. We will initialize these fields as `CellArray` objects:
 ```julia
-pPhases, pT         = init_cell_arrays(particles, Val(2))
-particle_args       = (pT, pPhases)
+pPhases, pT      = init_cell_arrays(particles, Val(2))
+particle_args    = (pT, pPhases)
 ```
 
 # Assign particles phases anomaly

diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
@@ -9,7 +9,7 @@ using JustPIC._2D
 # Threads is the default backend,
 # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script,
 # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # Load script dependencies
 using Printf, LinearAlgebra, GeoParams, GLMakie, CellArrays

diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl
@@ -9,7 +9,7 @@ using JustPIC._2D
 # Threads is the default backend,
 # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script,
 # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # Load script dependencies
 using Printf, LinearAlgebra, GeoParams, CairoMakie, CellArrays

diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl
@@ -9,7 +9,7 @@ using JustPIC._2D
 # Threads is the default backend,
 # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script,
 # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # Load script dependencies
 using Printf, LinearAlgebra, GeoParams, CairoMakie

diff --git a/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl b/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl
@@ -12,7 +12,7 @@ using JustPIC, JustPIC._2D
 # Threads is the default backend,
 # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script,
 # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # x-length of the domain
 const λ = 0.9142

diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl
@@ -2,7 +2,7 @@ using JustRelax, JustRelax.JustRelax2D
 const backend_JR = CPUBackend
 
 using JustPIC, JustPIC._2D
-const backend = CPUBackend
+const backend = JustPIC.CPUBackend
 
 using ParallelStencil, ParallelStencil.FiniteDifferences2D
 @init_parallel_stencil(Threads, Float64, 2)
@@ -109,7 +109,7 @@ function main(igg, nx, ny)
     # Initialize particles -------------------------------
     nxcell, max_xcell, min_xcell = 30, 40, 15
     particles = init_particles(
-        backend, nxcell, max_xcell, min_xcell, xvi[1], xvi[2], di[1], di[2], nx, ny
+        backend, nxcell, max_xcell, min_xcell, xvi, di, ni
     )
     # velocity grids
     grid_vx, grid_vy = velocity_grids(xci, xvi, di)
@@ -195,7 +195,7 @@ function main(igg, nx, ny)
             velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...)
             nt = 5
             fig = Figure(size = (900, 900), title = "t = $t")
-            ax  = Axis(fig[1,1], aspect = 1, title = " t=$(t/(1e3 * 3600 * 24 *365.25)) Kyrs")
+            ax  = Axis(fig[1,1], aspect = 1, title = " t=$(round.(t/(1e3 * 3600 * 24 *365.25); digits=3)) Kyrs")
             heatmap!(ax, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.viscosity.η)), colormap = :grayC)
             arrows!(
                 ax,

diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl
@@ -3,13 +3,13 @@ using JustRelax, JustRelax.JustRelax2D
 const backend_JR = CPUBackend
 
 using JustPIC, JustPIC._2D
-const backend = CPUBackend
+const backend = JustPIC.CPUBackend
 
 using ParallelStencil, ParallelStencil.FiniteDifferences2D
 @init_parallel_stencil(Threads, Float64, 2)
 
 # Load script dependencies
-using LinearAlgebra, GeoParams, GLMakie
+using LinearAlgebra, GeoParams, CairoMakie#GLMakie
 
 ## START OF HELPER FUNCTION ----------------------------------------------------------
 function copyinn_x!(A, B)
@@ -213,7 +213,7 @@ function RT_2D(igg, nx, ny)
             clr      = pPhases.data[:]
 
             fig = Figure(size = (900, 900), title = "t = $t")
-            ax  = Axis(fig[1,1], aspect = 1, title = " t=$(t/(1e3 * 3600 * 24 *365.25)) Kyrs")
+            ax  = Axis(fig[1,1], aspect = 1, title = " t=$(round.(t/(1e3 * 3600 * 24 *365.25); digits=3)) Kyrs")
             scatter!(
                 ax,
                 pxv, pyv,
@@ -237,7 +237,7 @@ end
 ## END OF MAIN SCRIPT ----------------------------------------------------------------
 
 # (Path)/folder where output data and figures are stored
-n        = 50
+n        = 100
 nx       = n
 ny       = n
 igg      = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid

diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D.jl
@@ -200,4 +200,4 @@ igg  = if !(JustRelax.MPI.Initialized())
 else
     igg
 end
-# main(igg; figdir = figdir, nx = nx, ny = ny);
+main(igg; figdir = figdir, nx = nx, ny = ny);
diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_MPI.jl
@@ -3,7 +3,11 @@ using JustRelax, JustRelax.JustRelax2D, GLMakie
 using ParallelStencil
 @init_parallel_stencil(Threads, Float64, 2)
 
-const backend = CPUBackend
+const backend_JR = CPUBackend
+
+using JustPIC, JustPIC._2D
+
+const backend = JustPIC.CPUBackend
 
 # HELPER FUNCTIONS ---------------------------------------------------------------
 solution(ε, t, G, η) = 2 * ε * η * (1 - exp(-G * t / η))
@@ -85,12 +89,12 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
 
     # Initialize phase ratios -------------------------------
     radius       = 0.1
-    phase_ratios = PhaseRatio(backend, ni, length(rheology))
+    phase_ratios = PhaseRatios(backend, length(rheology), ni)
     init_phases!(phase_ratios, xci, xvi, radius)
 
    # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
-    stokes    = StokesArrays(backend, ni)
+    stokes    = StokesArrays(backend_JR, ni)
     pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-6,  CFL = 0.75 / √2.1)
 
     # Buoyancy forces
@@ -107,8 +111,8 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
         free_slip = (left = true, right = true, top = true, bot = true),
         no_slip   = (left = false, right = false, top = false, bot=false),
     )
-    stokes.V.Vx   .= PTArray(backend)([ x*εbg for x in xvi[1], _ in 1:ny+2])
-    stokes.V.Vy   .= PTArray(backend)([-y*εbg for _ in 1:nx+2, y in xvi[2]])
+    stokes.V.Vx   .= PTArray(backend_JR)([ x*εbg for x in xvi[1], _ in 1:ny+2])
+    stokes.V.Vy   .= PTArray(backend_JR)([-y*εbg for _ in 1:nx+2, y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(@velocity(stokes)...)
 
@@ -204,7 +208,7 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
 
 end
 
-N      = 30 
+N      = 30
 n      = N
 nx     = n*2  # if only 2 CPU/GPU are used nx = 67 - 2 with N =128
 ny     = n*2

diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_displacement.jl
@@ -3,7 +3,11 @@ using JustRelax, JustRelax.JustRelax2D
 using ParallelStencil
 @init_parallel_stencil(Threads, Float64, 2)
 
-const backend = CPUBackend
+const backend_JR = CPUBackend
+
+using JustPIC, JustPIC._2D
+
+const backend = JustPIC.CPUBackend
 
 # HELPER FUNCTIONS ----------------------------------- ----------------------------
 solution(ε, t, G, η) = 2 * ε * η * (1 - exp(-G * t / η))
@@ -90,12 +94,12 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
 
     # Initialize phase ratios -------------------------------
     radius       = 0.1
-    phase_ratios = PhaseRatio(backend, ni, length(rheology))
+    phase_ratios = PhaseRatios(backend, length(rheology), ni)
     init_phases!(phase_ratios, xci, xvi, radius)
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
-    stokes    = StokesArrays(backend, ni)
+    stokes    = StokesArrays(backend_JR, ni)
     pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-6,  CFL = 0.75 / √2.1)
 
     # Buoyancy forces
@@ -111,8 +115,8 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
     free_slip = (left = true, right = true, top = true, bot = true),
     no_slip   = (left = false, right = false, top = false, bot=false),
 )
-    stokes.U.Ux .= PTArray(backend)([ x*εbg*lx*dt for x in xvi[1], _ in 1:ny+2])
-    stokes.U.Uy .= PTArray(backend)([-y*εbg*ly*dt for _ in 1:nx+2, y in xvi[2]])
+    stokes.U.Ux .= PTArray(backend_JR)([ x*εbg*lx*dt for x in xvi[1], _ in 1:ny+2])
+    stokes.U.Uy .= PTArray(backend_JR)([-y*εbg*ly*dt for _ in 1:nx+2, y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     displacement2velocity!(stokes, dt)
     update_halo!(@velocity(stokes)...)

diff --git a/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl b/miniapps/benchmarks/stokes2D/shear_band/ShearBand2D_softening.jl
@@ -3,7 +3,12 @@ using JustRelax, JustRelax.JustRelax2D
 using ParallelStencil
 @init_parallel_stencil(Threads, Float64, 2)
 
-const backend = CPUBackend
+const backend_JR = CPUBackend
+
+using JustPIC, JustPIC._2D
+
+const backend = JustPIC.CPUBackend
+
 # HELPER FUNCTIONS ---------------------------------------------------------------
 solution(ε, t, G, η) = 2 * ε * η * (1 - exp(-G * t / η))
 
@@ -87,12 +92,12 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
 
     # Initialize phase ratios -------------------------------
     radius       = 0.1
-    phase_ratios = PhaseRatio(backend, ni, length(rheology))
+    phase_ratios = PhaseRatios(backend, length(rheology), ni)
     init_phases!(phase_ratios, xci, xvi, radius)
 
      # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
-    stokes    = StokesArrays(backend, ni)
+    stokes    = StokesArrays(backend_JR, ni)
     pt_stokes = PTStokesCoeffs(li, di; ϵ=1e-6,  CFL = 0.9 / √2.1)
 
     # Buoyancy forces
@@ -108,8 +113,8 @@ function main(igg; nx=64, ny=64, figdir="model_figs")
         free_slip = (left = true, right = true, top = true, bot = true),
         no_slip   = (left = false, right = false, top = false, bot=false),
     )
-    stokes.V.Vx .= PTArray(backend)([ x*εbg for x in xvi[1], _ in 1:ny+2])
-    stokes.V.Vy .= PTArray(backend)([-y*εbg for _ in 1:nx+2, y in xvi[2]])
+    stokes.V.Vx .= PTArray(backend_JR)([ x*εbg for x in xvi[1], _ in 1:ny+2])
+    stokes.V.Vy .= PTArray(backend_JR)([-y*εbg for _ in 1:nx+2, y in xvi[2]])
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(stokes.V.Vx, stokes.V.Vy)
 
@@ -192,4 +197,3 @@ else
     igg
 end
 main(igg; figdir = figdir, nx = nx, ny = ny);
-
diff --git a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl
@@ -10,7 +10,7 @@ using JustPIC, JustPIC._2D
 # Threads is the default backend,
 # to run on a CUDA GPU load CUDA.jl (i.e. "using CUDA") at the beginning of the script,
 # and to run on an AMD GPU load AMDGPU.jl (i.e. "using AMDGPU") at the beginning of the script.
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # Load script dependencies
 using GeoParams, GLMakie
@@ -326,7 +326,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
   end
 
 figdir   = "Benchmark_Duretz_etal_2014"
-do_vtk = false # set to true to generate VTK files for ParaView
+do_vtk   = false # set to true to generate VTK files for ParaView
 ar       = 1 # aspect ratio
 n        = 64
 nx       = n*ar - 2

diff --git a/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl b/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl
@@ -5,9 +5,9 @@ using ParallelStencil, ParallelStencil.FiniteDifferences2D
 @init_parallel_stencil(Threads, Float64, 2)
 
 using JustPIC, JustPIC._2D
-const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = JustPIC.CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
-using GeoParams
+using GeoParams, GLMakie
 
 ## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT --------------------------------
 
@@ -171,7 +171,7 @@ function sinking_block2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_per
     # Plotting ---------------------
     f, _, h = heatmap(velocity, colormap=:vikO)
     Colorbar(f[1,2], h)
-    f
+    display(f)
     # ------------------------------
 
     return nothing

diff --git a/miniapps/benchmarks/stokes2D/solvi/vizSolVi.jl b/miniapps/benchmarks/stokes2D/solvi/vizSolVi.jl
@@ -61,7 +61,7 @@ function solvi_solution(geometry, η0, ηi, εbg, rc)
     return (p=ps, vx=-vxs, vy=-vys)
 end
 
-function Li_error(geometry, stokes::StokesArrays, Δη, εbg, rc, ; order=2)
+function Li_error(geometry, stokes::JustRelax.StokesArrays, Δη, εbg, rc, ; order=2)
 
     # analytical solution
     sol = solvi_solution(geometry, 1, Δη, εbg, rc)
@@ -76,7 +76,7 @@ function Li_error(geometry, stokes::StokesArrays, Δη, εbg, rc, ; order=2)
     return L2_vx, L2_vy, L2_p
 end
 
-function plot_solVi_error(geometry, stokes::StokesArrays, Δη, εbg, rc)
+function plot_solVi_error(geometry, stokes::JustRelax.StokesArrays, Δη, εbg, rc)
 
     # analytical solution
     sol = solvi_solution(geometry, 1, Δη, εbg, rc)

diff --git a/miniapps/benchmarks/stokes3D/RunStokesBench3D.jl b/miniapps/benchmarks/stokes3D/RunStokesBench3D.jl
@@ -1,15 +1,14 @@
 using LinearAlgebra, CairoMakie
-using JustRelax
+using JustRelax, JustRelax.JustRelax3D
 using MPI: MPI
+
 using ParallelStencil
 @init_parallel_stencil(Threads, Float64, 3)
 
-# setup ParallelStencil.jl environment
-model = PS_Setup(:cpu, Float64, 3)
-environment!(model)
+const backend_JR = CPUBackend
 
 # choose benchmark
-benchmark = :solvi
+benchmark = :Burstedde
 
 # model resolution (number of gridpoints)
 nx, ny, nz = 16, 16, 16

diff --git a/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl b/miniapps/benchmarks/stokes3D/burstedde/Burstedde.jl
@@ -195,7 +195,7 @@ function burstedde(; nx=16, ny=16, nz=16, init_MPI=true, finalize_MPI=false)
 
     ## Allocate arrays needed for every Stokes problem
     # general stokes arrays
-    stokes    = StokesArrays(ni, ViscoElastic)
+    stokes    = StokesArrays(backend, ni)
     # general numerical coeffs for PT stokes
     pt_stokes = PTStokesCoeffs(li, di; CFL=1 / √3)
 

diff --git a/miniapps/benchmarks/stokes3D/burstedde/vizBurstedde.jl b/miniapps/benchmarks/stokes3D/burstedde/vizBurstedde.jl
@@ -52,7 +52,7 @@ function analytical_solution(xci, xvi)
     return Vx, Vy, Vz, P
 end
 
-function plot(stokes::StokesArrays, geometry; cmap=:vik)
+function plot(stokes::JustRelax.StokesArrays, geometry; cmap=:vik)
     xci, xvi = geometry.xci, geometry.xvi
     vx, vy, vz, p = analytical_solution(xci, xvi)
 
@@ -92,7 +92,7 @@ function plot(stokes::StokesArrays, geometry; cmap=:vik)
     return f
 end
 
-function error(stokes::StokesArrays, geometry)
+function error(stokes::JustRelax.StokesArrays, geometry)
     gridsize = foldl(*, geometry.di)
     vx, vy, vz, p = analytical_solution(geometry.xci, geometry.xvi)
-Original file line number
+Diff line change
@@ Expand Up / @@ -15,6 +15,7 @@ jobs: @@
           matrix:
             version:
               - '1.10'
+              - '1.11'
               - 'pre'
               # - 'nightly'
             os:
@@ Expand Down @@