diff --git a/docs/make.jl b/docs/make.jl
index 9159dcc7..ca6f4fc1 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -11,7 +11,7 @@ makedocs(;
 
     warnonly = Documenter.except(:footnote),
     pages=[
-        "Home"      => "man/index.md",
+        "Home"      => "index.md",
         "User guide"=> Any[
             "Installation" => "man/installation.md",
             "Backend"   => "man/backend.md",
diff --git a/docs/src/assets/Blankenbach_Temp.png b/docs/src/assets/Blankenbach_Temp.png
new file mode 100644
index 00000000..0f473885
Binary files /dev/null and b/docs/src/assets/Blankenbach_Temp.png differ
diff --git a/docs/src/assets/Blankenbach_Time_Series_V_Nu.png b/docs/src/assets/Blankenbach_Time_Series_V_Nu.png
new file mode 100644
index 00000000..c219a119
Binary files /dev/null and b/docs/src/assets/Blankenbach_Time_Series_V_Nu.png differ
diff --git a/docs/src/man/index.md b/docs/src/index.md
similarity index 95%
rename from docs/src/man/index.md
rename to docs/src/index.md
index 95a5ed48..233233fc 100644
--- a/docs/src/man/index.md
+++ b/docs/src/index.md
@@ -18,7 +18,7 @@ The package serves several purposes:
   * It provides a collection of solvers to be used in quickly developing new applications
   * It provides some standardization so that application codes can
 
-     - more easily handle local material properties through the use of [GeoParams.jl]((https://github.com/JuliaGeodynamics/GeoParams.jl))
+     - more easily handle local material properties through the use of [GeoParams.jl](https://github.com/JuliaGeodynamics/GeoParams.jl)
      - more easily switch between a pseudo-transient solver and another solvers (e.g. an explicit thermal solvers)
 
   * It provides a natural repository for contributions of new solvers for use by the larger community
diff --git a/docs/src/man/Blankenbach.md b/docs/src/man/Blankenbach.md
index eb358fd3..9bb4f8c2 100644
--- a/docs/src/man/Blankenbach.md
+++ b/docs/src/man/Blankenbach.md
@@ -392,7 +392,7 @@ fig
 
 ### Final model
 Temperature field
-![Temperatuere](../assets/Blankenbach/Temp.png)
+![Temperature](../assets/Blankenbach_Temp.png)
 
 And time history of the rms-velocity and Nusselt number
-![time series](../assets/Blankenbach/Time_Series_V_Nu.png)
+![time series](../assets/Blankenbach_Time_Series_V_Nu.png)
diff --git a/docs/src/man/listfunctions.md b/docs/src/man/listfunctions.md
index 181a8a49..7e98412f 100644
--- a/docs/src/man/listfunctions.md
+++ b/docs/src/man/listfunctions.md
@@ -2,5 +2,5 @@
 
 Here an overview of all functions:
 ```@autodocs
-Modules = [JustRelax]
+Modules = [JustRelax, JustRelax.JustRelax2D, JustRelax.JustRelax3D, JustRelax.DataIO]
 ```
diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
index 3bf7c5a9..c87b1812 100755
--- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
+++ b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
@@ -85,7 +85,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     pPhases,            = init_cell_arrays(particles, Val(1))
     phase_ratios        = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl
index 7f3868d8..8b9593a7 100755
--- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl
+++ b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd.jl
@@ -86,7 +86,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     particle_args       = (pT, pT0, pPhases)
     phase_ratios        = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -257,7 +257,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer, ), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         # Nusselt number, Nu = H/ΔT/L ∫ ∂T/∂z dx ----
         Nu_it   =   (ly / (1000.0*lx)) *
diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl
index 4233b273..1c5e8e5d 100755
--- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl
+++ b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_sgd_scaled.jl
@@ -78,7 +78,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     particle_args       = (pT, pT0, pPhases)
     phase_ratios        = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -247,7 +247,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer, ), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         # Nusselt number, Nu = ∫ ∂T/∂z dx ----
         Nu_it   =   sum( ((abs.(thermal.T[2:end-1,end] - thermal.T[2:end-1,end-1])) ./ di[2]) .*di[1])
diff --git a/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl b/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl
index e7a0dc0b..8236e35d 100644
--- a/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl
+++ b/miniapps/benchmarks/stokes2D/VanKeken.jl/VanKeken.jl
@@ -86,7 +86,7 @@ function main2D(igg; ny=64, nx=64, figdir="model_figs")
     particle_args        = (pPhases, )
     phase_ratios         = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
@@ -170,7 +170,7 @@ function main2D(igg; ny=64, nx=64, figdir="model_figs")
         # inject && break
         inject_particles_phase!(particles, pPhases, (), (), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl
index bfe0d97f..393b7f59 100644
--- a/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl
+++ b/miniapps/benchmarks/stokes2D/free_surface_stabilization/PlumeFreeSurface_2D.jl
@@ -120,7 +120,7 @@ function main(igg, nx, ny)
     # Elliptical temperature anomaly
     init_phases!(pPhases, particles)
     phase_ratios  = PhaseRatio(backend_JR, ni, length(rheology))
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -186,7 +186,7 @@ function main(igg, nx, ny)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (), (), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl b/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl
index f864a07f..a96b43da 100644
--- a/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl
+++ b/miniapps/benchmarks/stokes2D/free_surface_stabilization/RayleighTaylor2D.jl
@@ -118,7 +118,7 @@ function RT_2D(igg, nx, ny)
     A             = 5e3    # Amplitude of the anomaly
     phase_ratios  = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, A)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -197,7 +197,7 @@ function RT_2D(igg, nx, ny)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (), (), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl
index fd83a6b4..87b61ec6 100644
--- a/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl
+++ b/miniapps/benchmarks/stokes2D/shear_heating/Shearheating2D.jl
@@ -81,7 +81,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     r_anomaly        = 3e3    # radius of perturbation
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, xc_anomaly, yc_anomaly, r_anomaly)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -245,7 +245,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
             @show it += 1
             t        += dt
diff --git a/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl b/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl
index a9799b9b..a2352737 100644
--- a/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl
+++ b/miniapps/benchmarks/stokes2D/sinking_block/SinkingBlock2D.jl
@@ -115,7 +115,7 @@ function sinking_block2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", thermal_per
     r_anomaly    =  50e3   # radius of perturbation
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, xc_anomaly, abs(yc_anomaly), r_anomaly)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
diff --git a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl
index e8f7a0be..89d8914f 100644
--- a/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl
+++ b/miniapps/benchmarks/stokes3D/shear_heating/Shearheating3D.jl
@@ -75,7 +75,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
     r_anomaly        = 3e3    # radius of perturbation
     init_phases!(backend_JR, pPhases, particles, xc_anomaly, yc_anomaly, zc_anomaly, r_anomaly)
     phase_ratios     = PhaseRatio(ni, length(rheology))
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -227,7 +227,7 @@ function main3D(igg; ar=8, ny=16, nx=ny*8, nz=ny*8, figdir="figs3D", do_vtk =fal
             # check if we need to inject particles
             inject_particles_phase!(particles, pPhases, (pT, ), (thermal.T,), xvi)
             # update phase ratios
-            phase_ratios_center(phase_ratios, particles, grid, pPhases)
+            phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
             @show it += 1
             t        += dt
diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl
index 3d12b48d..c4dc5efe 100644
--- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl
+++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase.jl
@@ -126,7 +126,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0)
     pPhases,     = init_cell_arrays(particles, Val(1))
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, center_perturbation..., r)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios.center, args)
diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_multiphase.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_multiphase.jl
index 4f5387e8..29901841 100644
--- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_multiphase.jl
+++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion3D_multiphase.jl
@@ -145,7 +145,7 @@ function diffusion_3D(;
     pPhases,     = init_cell_arrays(particles, Val(1))
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, center_perturbation..., r)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # PT coefficients for thermal diffusion
diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl
index 94132853..b2e2b261 100644
--- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl
+++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim.jl
@@ -241,7 +241,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
     anomaly      = nondimensionalize((750 + 273)K, CharDim) # thermal perturbation (in K)
     init_phases!(pPhases, particles, x_anomaly, y_anomaly, r_anomaly, sticky_air, nondimensionalize(0.0km,CharDim), nondimensionalize(20km,CharDim))
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # Initialisation of thermal profile
     thermal     = ThermalArrays(backend_JR, ni) # initialise thermal arrays and boundary conditions
@@ -438,7 +438,7 @@ function main2D(igg; figdir=figdir, nx=nx, ny=ny, do_vtk=false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         particle2grid!(T_buffer, pT, xvi, particles)
         @views T_buffer[:, end]        .= Tsurf
diff --git a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl
index 3d38c660..d071e8d7 100644
--- a/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl
+++ b/miniapps/benchmarks/thermal_stress/Thermal_Stress_Magma_Chamber_nondim3D.jl
@@ -242,7 +242,7 @@ function main3D(igg; figdir = "output", nx = 64, ny = 64, nz = 64, do_vtk = fals
     anomaly      = nondimensionalize((750 + 273)K, CharDim)               # thermal perturbation (in K)
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, x_anomaly, y_anomaly, z_anomaly, r_anomaly, sticky_air, nondimensionalize(0.0km,CharDim), nondimensionalize(20km,CharDim))
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # Initialisation of thermal profile
     thermal     = ThermalArrays(backend_JR, ni) # initialise thermal arrays and boundary conditions
@@ -420,7 +420,7 @@ function main3D(igg; figdir = "output", nx = 64, ny = 64, nz = 64, do_vtk = fals
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (thermal.T,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         particle2grid!(thermal.T, pT, xvi, particles)
         @views thermal.T[:, :, end] .= Tsurf
diff --git a/miniapps/convection/Particles2D/Layered_convection2D.jl b/miniapps/convection/Particles2D/Layered_convection2D.jl
index 441c402e..69822f38 100644
--- a/miniapps/convection/Particles2D/Layered_convection2D.jl
+++ b/miniapps/convection/Particles2D/Layered_convection2D.jl
@@ -128,7 +128,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     r_anomaly        = 25e3    # radius of perturbation
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, lx, yc_anomaly, r_anomaly, thick_air)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -292,7 +292,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
index 560fd1ea..4228ed60 100644
--- a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
+++ b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
@@ -135,7 +135,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     r_anomaly        = nondimensionalize(25km, CharDim) # radius of perturbation
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, lx, yc_anomaly, r_anomaly, thick_air, CharDim)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -298,7 +298,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/convection/Particles3D/Layered_convection3D.jl b/miniapps/convection/Particles3D/Layered_convection3D.jl
index b3878993..7ab0475e 100644
--- a/miniapps/convection/Particles3D/Layered_convection3D.jl
+++ b/miniapps/convection/Particles3D/Layered_convection3D.jl
@@ -115,7 +115,7 @@ function main3D(igg; ar=1, nx=16, ny=16, nz=16, figdir="figs3D", do_vtk =false)
     r_anomaly        = 50e3   # radius of perturbation
     init_phases!(pPhases, particles, lx, ly; d=abs(zc_anomaly), r=r_anomaly)
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
-    @parallel (@idx ni) phase_ratios_center(phase_ratios.center, particles.coords, xci, di, pPhases)
+    @parallel (@idx ni) phase_ratios_center!(phase_ratios.center, particles.coords, xci, di, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -264,7 +264,7 @@ function main3D(igg; ar=1, nx=16, ny=16, nz=16, figdir="figs3D", do_vtk =false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/miniapps/convection/Rifting2D/Rifting2D.jl b/miniapps/convection/Rifting2D/Rifting2D.jl
new file mode 100644
index 00000000..3e9c82d1
--- /dev/null
+++ b/miniapps/convection/Rifting2D/Rifting2D.jl
@@ -0,0 +1,399 @@
+using JustRelax, JustRelax.DataIO
+import JustRelax.@cell
+using ParallelStencil
+@init_parallel_stencil(Threads, Float64, 2) #or (CUDA, Float64, 2) or (AMDGPU, Float64, 2)
+
+using JustPIC
+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
+
+# setup ParallelStencil.jl environment
+model = PS_Setup(:Threads, Float64, 2) #or (:CUDA, Float64, 2) or (:AMDGPU, Float64, 2)
+environment!(model)
+
+# Load script dependencies
+using Printf, LinearAlgebra, GeoParams, GLMakie, CellArrays
+
+# Load file with all the rheology configurations
+include("Rifting_rheology.jl")
+
+## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT --------------------------------
+
+function copyinn_x!(A, B)
+
+    @parallel function f_x(A, B)
+        @all(A) = @inn_x(B)
+        return nothing
+    end
+
+    @parallel f_x(A, B)
+end
+
+import ParallelStencil.INDICES
+const idx_j = INDICES[2]
+macro all_j(A)
+    esc(:($A[$idx_j]))
+end
+
+# Initial pressure profile - not accurate
+@parallel function init_P!(P, ρg, z)
+    @all(P) = abs(@all(ρg) * @all_j(z)) * <(@all_j(z), 0.0)
+    return nothing
+end
+
+# Initial thermal profile
+@parallel_indices (i, j) function init_T!(T, y, thick_air)
+    depth = -y[j] - thick_air
+
+    # (depth - 15e3) because we have 15km of sticky air
+    if depth < 0e0
+        T[i + 1, j] = 273e0
+
+    elseif 0e0 ≤ (depth) < 35e3
+        dTdZ        = (923-273)/35e3
+        offset      = 273e0
+        T[i + 1, j] = (depth) * dTdZ + offset
+
+    elseif 110e3 > (depth) ≥ 35e3
+        dTdZ        = (1492-923)/75e3
+        offset      = 923
+        T[i + 1, j] = (depth - 35e3) * dTdZ + offset
+
+    elseif (depth) ≥ 110e3
+        dTdZ        = (1837 - 1492)/590e3
+        offset      = 1492e0
+        T[i + 1, j] = (depth - 110e3) * dTdZ + offset
+
+    end
+
+    return nothing
+end
+
+# Thermal rectangular perturbation
+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)
+            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...)
+
+    return nothing
+end
+## END OF HELPER FUNCTION ------------------------------------------------------------
+
+function pureshear!(stokes, εbg, xvi)
+    midpoint = (xvi[1][end] - xvi[1][1]) / 2
+    stokes.V.Vx .= PTArray([ (x-midpoint)*εbg for x in xvi[1], _ in 1:ny+2])
+    stokes.V.Vy .= PTArray([-y*εbg for _ in 1:nx+2, y in xvi[2]])
+    return nothing
+end
+
+## BEGIN OF MAIN SCRIPT --------------------------------------------------------------
+function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
+
+    # Physical domain ------------------------------------
+    thick_air    = 0                 # thickness of sticky air layer
+    ly           = 300e3 + thick_air # domain length in y
+    lx           = ly * ar           # 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
+    origin       = 0.0, -ly          # origin coordinates (15km f sticky air layer)
+    grid         = Geometry(ni, li; origin = origin)
+    (; xci, xvi) = grid # nodes at the center and vertices of the cells
+    # ----------------------------------------------------
+
+    # Physical properties using GeoParams ----------------
+    εbg          = 1e-15 # background strain rate
+    rheology     = init_rheologies(; is_plastic = true)
+    κ            = (4 / (rheology[4].HeatCapacity[1].Cp * rheology[4].Density[1].ρ0))
+    dt = dt_diff = 0.5 * min(di...)^2 / κ / 2.01 # diffusive CFL timestep limiter
+    # ----------------------------------------------------
+
+    # Initialize particles -------------------------------
+    nxcell, max_xcell, min_xcell = 25, 30, 8
+    particles        = init_particles(
+        backend, nxcell, max_xcell, min_xcell, xvi..., di..., ni...
+    )
+    subgrid_arrays   = SubgridDiffusionCellArrays(particles)
+    # velocity grids
+    grid_vx, grid_vy = velocity_grids(xci, xvi, di)
+    # temperature
+    pT, pPhases      = init_cell_arrays(particles, Val(2))
+    particle_args    = (pT, pPhases)
+
+    # Elliptical temperature anomaly
+    xc_anomaly       =  lx/2   # origin of thermal anomaly
+    yc_anomaly       = -120e3  # origin of thermal anomaly
+    r_anomaly        =  5e3    # radius of perturbation
+    init_phases!(pPhases, particles, lx, yc_anomaly, r_anomaly, thick_air)
+    phase_ratios     = PhaseRatio(ni, length(rheology))
+    @parallel (@idx ni) phase_ratios_center!(phase_ratios.center, pPhases)
+    # ----------------------------------------------------
+
+    # STOKES ---------------------------------------------
+    # Allocate arrays needed for every Stokes problem
+    stokes           = StokesArrays(ni, ViscoElastic)
+    pt_stokes        = PTStokesCoeffs(li, di; ϵ=1e-4,  CFL = 0.1 / √2.1)
+    pureshear!(stokes, εbg * 0, xvi)
+    # 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)
+    # ----------------------------------------------------
+
+    # TEMPERATURE PROFILE --------------------------------
+    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 - Half space cooling
+    @parallel (@idx ni .+ 1) init_T!(thermal.T, xvi[2], thick_air)
+    thermal_bcs!(thermal.T, thermal_bc)
+    Tbot = thermal.T[1, 1]
+    rectangular_perturbation!(thermal.T, xc_anomaly, yc_anomaly, r_anomaly, xvi, thick_air)
+    @parallel (JustRelax.@idx size(thermal.Tc)...) temperature2center!(thermal.Tc, thermal.T)
+    # ----------------------------------------------------
+
+    # Buoyancy forces
+    ρg               = @zeros(ni...), @zeros(ni...)
+    for _ in 1:1
+        @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])
+    end
+    # Rheology
+    η                = @ones(ni...)
+    args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
+    @parallel (@idx ni) compute_viscosity!(
+        η, 1.0, phase_ratios.center, @strain(stokes)..., args, rheology, (1e18, 1e24)
+    )
+    η_vep            = copy(η)
+
+    # PT coefficients for thermal diffusion
+    pt_thermal       = PTThermalCoeffs(
+        rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL= 1e-2 / √2.1
+    )
+
+    # IO ------------------------------------------------
+    # if it does not exist, make folder where figures are stored
+    if do_vtk
+        vtk_dir      = joinpath(figdir, "vtk")
+        take(vtk_dir)
+    end
+    take(figdir)
+    # ----------------------------------------------------
+
+    # Plot initial T and η profiles
+    let
+        Yv  = [y for x in xvi[1], y in xvi[2]][:]
+        Y   = [y for x in xci[1], y in xci[2]][:]
+        fig = Figure(size = (1200, 900))
+        ax1 = Axis(fig[1,1], aspect = 2/3, title = "T")
+        ax2 = Axis(fig[1,2], aspect = 2/3, title = "log10(η)")
+        scatter!(ax1, Array(thermal.T[2:end-1,:][:]), Yv./1e3)
+        scatter!(ax2, Array(log10.(η[:])), Y./1e3)
+        ylims!(ax1, minimum(xvi[2])./1e3, 0)
+        ylims!(ax2, minimum(xvi[2])./1e3, 0)
+        hideydecorations!(ax2)
+        save(joinpath(figdir, "initial_profile.png"), fig)
+        fig
+    end
+
+    T_buffer    = @zeros(ni.+1)
+    Told_buffer = similar(T_buffer)
+    dt₀         = similar(stokes.P)
+    for (dst, src) in zip((T_buffer, Told_buffer), (thermal.T, thermal.Told))
+        copyinn_x!(dst, src)
+    end
+    grid2particle!(pT, xvi, T_buffer, particles)
+
+    local Vx_v, Vy_v
+    if do_vtk
+        Vx_v = @zeros(ni.+1...)
+        Vy_v = @zeros(ni.+1...)
+    end
+    # Time loop
+    t, it = 0.0, 0
+    while (t/(1e6 * 3600 * 24 *365.25)) < 5 # run only for 5 Myrs
+
+        # interpolate fields from particle to grid vertices
+        particle2grid!(T_buffer, pT, xvi, particles)
+        @views T_buffer[:, end]      .= 273.0
+        @views T_buffer[:, 1]        .= Tbot
+        @views thermal.T[2:end-1, :] .= T_buffer
+        thermal_bcs!(thermal.T, thermal_bc)
+        temperature2center!(thermal)
+ 
+        # Update buoyancy and viscosity -
+        args = (; T = thermal.Tc, P = stokes.P,  dt=Inf)
+        @parallel (@idx ni) compute_viscosity!(
+            η, 1.0, phase_ratios.center, @strain(stokes)..., args, rheology, (1e18, 1e24)
+        )
+        @parallel (JustRelax.@idx ni) compute_ρg!(ρg[2], phase_ratios.center, rheology, args)
+        # ------------------------------
+
+        # Stokes solver ----------------
+        @edit solve!(
+            stokes,
+            pt_stokes,
+            di,
+            flow_bcs,
+            ρg,
+            η,
+            η_vep,
+            phase_ratios,
+            rheology,
+            args,
+            0.1,
+            igg;
+            iterMax = 1e3,
+            nout = 1e1,
+            viscosity_cutoff=(1e18, 1e24)
+        )
+        @parallel (JustRelax.@idx ni) JustRelax.Stokes2D.tensor_invariant!(stokes.ε.II, @strain(stokes)...)
+        dt   = compute_dt(stokes, di, dt_diff)
+        # ------------------------------
+
+        # Thermal solver ---------------
+        heatdiffusion_PT!(
+            thermal,
+            pt_thermal,
+            thermal_bc,
+            rheology,
+            args,
+            dt,
+            di;
+            igg     = igg,
+            phase   = phase_ratios,
+            iterMax = 10e3,
+            nout    = 1e2,
+            verbose = true,
+        )
+        for (dst, src) in zip((T_buffer, Told_buffer), (thermal.T, thermal.Told))
+            copyinn_x!(dst, src)
+        end
+        subgrid_characteristic_time!(
+            subgrid_arrays, particles, dt₀, phase_ratios, rheology, thermal, stokes, xci, di
+        )
+        centroid2particle!(subgrid_arrays.dt₀, xci, dt₀, particles)
+        subgrid_diffusion!(
+            pT, T_buffer, thermal.ΔT[2:end-1, :], subgrid_arrays, particles, xvi,  di, dt
+        )
+        # ------------------------------
+
+        # Advection --------------------
+        # 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)
+        # check if we need to inject particles
+        inject = check_injection(particles)
+        inject && inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
+        # update phase ratios
+        @parallel (@idx ni) phase_ratios_center!(phase_ratios.center, pPhases)
+
+        @show it += 1
+        t        += dt
+
+        # Data I/O and plotting ---------------------
+        if it == 1 || rem(it, 1) == 0
+            checkpointing(figdir, stokes, thermal.T, η, t)
+
+            if do_vtk
+                JustRelax.velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...)
+                data_v = (;
+                    T   = Array(thermal.T[2:end-1, :]),
+                    τxy = Array(stokes.τ.xy),
+                    εxy = Array(stokes.ε.xy),
+                    Vx  = Array(Vx_v),
+                    Vy  = Array(Vy_v),
+                )
+                data_c = (;
+                    P   = Array(stokes.P),
+                    τxx = Array(stokes.τ.xx),
+                    τyy = Array(stokes.τ.yy),
+                    εxx = Array(stokes.ε.xx),
+                    εyy = Array(stokes.ε.yy),
+                    η   = Array(η),
+                )
+                do_vtk(
+                    joinpath(vtk_dir, "vtk_" * lpad("$it", 6, "0")),
+                    xvi,
+                    xci,
+                    data_v,
+                    data_c
+                )
+            end
+
+            # Make particles plottable
+            p        = particles.coords
+            ppx, ppy = p
+            pxv      = ppx.data[:]./1e3
+            pyv      = ppy.data[:]./1e3
+            clr      = pPhases.data[:]
+            clr      = pT.data[:]
+            idxv     = particles.index.data[:];
+
+            # Make Makie figure
+            fig = Figure(size = (900, 900), title = "t = $t")
+            ax1 = Axis(fig[1,1], aspect = ar, title = "T [K]  (t=$(t/(1e6 * 3600 * 24 *365.25)) Myrs)")
+            ax2 = Axis(fig[2,1], aspect = ar, title = "Vy [m/s]")
+            ax3 = Axis(fig[1,3], aspect = ar, title = "log10(εII)")
+            ax4 = Axis(fig[2,3], aspect = ar, title = "log10(η)")
+            # Plot temperature
+            h1  = heatmap!(ax1, xvi[1].*1e-3, xvi[2].*1e-3, Array(thermal.T[2:end-1,:]) , colormap=:batlow)
+            # Plot particles phase
+            h2  = scatter!(ax2, Array(pxv[idxv]), Array(pyv[idxv]), color=Array(clr[idxv]))
+            # Plot 2nd invariant of strain rate
+            h3  = heatmap!(ax3, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(stokes.ε.II)) , colormap=:batlow)
+            # Plot effective viscosity
+            h4  = heatmap!(ax4, xci[1].*1e-3, xci[2].*1e-3, Array(log10.(η_vep)) , colormap=:batlow)
+            hidexdecorations!(ax1)
+            hidexdecorations!(ax2)
+            hidexdecorations!(ax3)
+            Colorbar(fig[1,2], h1)
+            Colorbar(fig[2,2], h2)
+            Colorbar(fig[1,4], h3)
+            Colorbar(fig[2,4], h4)
+            linkaxes!(ax1, ax2, ax3, ax4)
+            save(joinpath(figdir, "$(it).png"), fig)
+            fig
+        end
+        # ------------------------------
+
+    end
+
+    return nothing
+end
+## END OF MAIN SCRIPT ----------------------------------------------------------------
+
+# (Path)/folder where output data and figures are stored
+figdir   = "Plume2D"
+do_vtk = false # set to true to generate VTK files for ParaView
+ar       = 1 # aspect ratio
+n        = 128
+nx       = n*ar - 2
+ny       = n - 2
+igg      = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
+    IGG(init_global_grid(nx, ny, 1; init_MPI= true)...)
+else
+    igg
+end
+
+# run main script
+
+# main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
\ No newline at end of file
diff --git a/miniapps/convection/Rifting2D/Rifting_rheology.jl b/miniapps/convection/Rifting2D/Rifting_rheology.jl
new file mode 100644
index 00000000..c5a6f6e6
--- /dev/null
+++ b/miniapps/convection/Rifting2D/Rifting_rheology.jl
@@ -0,0 +1,149 @@
+# from "Fingerprinting secondary mantle plumes", Cloetingh et al. 2022
+
+function init_rheologies(; is_plastic = true)
+
+    # Dislocation and Diffusion creep
+    disl_upper_crust            = DislocationCreep(A=5.07e-18, n=2.3, E=154e3, V=6e-6,  r=0.0, R=8.3145)
+    disl_lower_crust            = DislocationCreep(A=2.08e-23, n=3.2, E=238e3, V=6e-6,  r=0.0, R=8.3145)
+    disl_lithospheric_mantle    = DislocationCreep(A=2.51e-17, n=3.5, E=530e3, V=6e-6,  r=0.0, R=8.3145)
+    disl_sublithospheric_mantle = DislocationCreep(A=2.51e-17, n=3.5, E=530e3, V=6e-6,  r=0.0, R=8.3145)
+    diff_lithospheric_mantle    = DislocationCreep(A=2.51e-17, n=1.0, E=530e3, V=6e-6,  r=0.0, R=8.3145)
+    diff_sublithospheric_mantle = DislocationCreep(A=2.51e-17, n=1.0, E=530e3, V=6e-6,  r=0.0, R=8.3145)
+
+    # Elasticity
+    el_upper_crust              = SetConstantElasticity(; G=25e9, ν=0.5)
+    el_lower_crust              = SetConstantElasticity(; G=25e9, ν=0.5)
+    el_lithospheric_mantle      = SetConstantElasticity(; G=67e9, ν=0.5)
+    el_sublithospheric_mantle   = SetConstantElasticity(; G=67e9, ν=0.5)
+    β_upper_crust               = inv(get_Kb(el_upper_crust))
+    β_lower_crust               = inv(get_Kb(el_lower_crust))
+    β_lithospheric_mantle       = inv(get_Kb(el_lithospheric_mantle))
+    β_sublithospheric_mantle    = inv(get_Kb(el_sublithospheric_mantle))
+
+    # Physical properties using GeoParams ----------------
+    η_reg     = 1e16
+    cohesion  = 3e6
+    friction  = asind(0.2)
+    pl_crust  = if is_plastic
+        DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    else
+        DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    end
+    friction  = asind(0.3)
+    pl        = if is_plastic
+        DruckerPrager_regularised(; C = cohesion, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    else
+        DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    end
+    pl_wz     = if is_plastic
+        DruckerPrager_regularised(; C = 2e6, ϕ=2.0, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    else
+        DruckerPrager_regularised(; C = Inf, ϕ=friction, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+    end
+
+    # crust
+    K_crust = TP_Conductivity(;
+        a = 0.64,
+        b = 807e0,
+        c = 0.77,
+        d = 0.00004*1e-6,
+    )
+
+    K_mantle = TP_Conductivity(;
+        a = 0.73,
+        b = 1293e0,
+        c = 0.77,
+        d = 0.00004*1e-6,
+    )
+
+    # Define rheolgy struct
+    rheology = (
+        # Name              = "UpperCrust",
+        SetMaterialParams(;
+            Phase             = 1,
+            Density           = PT_Density(; ρ0=2.75e3, β=β_upper_crust, T0=0.0, α = 3.5e-5),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=7.5e2),
+            Conductivity      = K_crust,
+            CompositeRheology = CompositeRheology((disl_upper_crust, el_upper_crust, pl_crust)),
+            Elasticity        = el_upper_crust,
+            RadioactiveHeat   = ConstantRadioactiveHeat(0.0),
+            Gravity           = ConstantGravity(; g=9.81),
+        ),
+        # Name              = "LowerCrust",
+        SetMaterialParams(;
+            Phase             = 2,
+            Density           = PT_Density(; ρ0=3e3, β=β_lower_crust, T0=0.0, α = 3.5e-5),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=7.5e2),
+            Conductivity      = K_crust,
+            RadioactiveHeat   = ConstantRadioactiveHeat(0.0),
+            CompositeRheology = CompositeRheology((disl_lower_crust, el_lower_crust, pl_crust)),
+            Elasticity        = el_lower_crust,
+        ),
+        # Name              = "LithosphericMantle",
+        SetMaterialParams(;
+            Phase             = 3,
+            Density           = PT_Density(; ρ0=3.3e3, β=β_lithospheric_mantle, T0=0.0, α = 3e-5),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=1.25e3),
+            Conductivity      = K_mantle,
+            RadioactiveHeat   = ConstantRadioactiveHeat(0.0),
+            CompositeRheology = CompositeRheology((disl_lithospheric_mantle, diff_lithospheric_mantle, el_lithospheric_mantle, pl)),
+            Elasticity        = el_lithospheric_mantle,
+        ),
+        SetMaterialParams(;
+            Phase             = 4,
+            Density           = PT_Density(; ρ0=3.3e3-50, β=β_sublithospheric_mantle, T0=0.0, α = 3e-5),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=1.25e3),
+            Conductivity      = K_mantle,
+            RadioactiveHeat   = ConstantRadioactiveHeat(0.0),
+            CompositeRheology = CompositeRheology((disl_sublithospheric_mantle, diff_sublithospheric_mantle, el_sublithospheric_mantle)),
+            Elasticity        = el_sublithospheric_mantle,
+        ),
+        # Name              = "StickyAir",
+        SetMaterialParams(;
+            Phase             = 5,
+            Density           = ConstantDensity(; ρ=1e3), # water density
+            HeatCapacity      = ConstantHeatCapacity(; Cp=3e3),
+            RadioactiveHeat   = ConstantRadioactiveHeat(0.0),
+            Conductivity      = ConstantConductivity(; k=1.0),
+            CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
+        ),
+    )
+end
+
+function init_phases!(phases, particles, Lx, d, r, thick_air)
+    ni = size(phases)
+
+    @parallel_indices (i, j) function init_phases!(phases, px, py, index, r, Lx)
+        @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]) - thick_air
+            if 0e0 ≤ depth ≤ 21e3
+                @cell phases[ip, i, j] = 1.0
+
+            elseif 35e3 ≥ depth > 21e3
+                @cell phases[ip, i, j] = 2.0
+
+            elseif 90e3 ≥ depth > 35e3
+                @cell phases[ip, i, j] = 3.0
+
+            elseif depth > 90e3
+                @cell phases[ip, i, j] = 3.0
+
+            elseif depth < 0e0
+                @cell phases[ip, i, j] = 5.0
+
+            end
+
+            # plume - rectangular
+            if ((x - Lx * 0.5)^2 ≤ r^2) && (((JustRelax.@cell py[ip, i, j]) - d - thick_air)^2 ≤ r^2)
+                JustRelax.@cell phases[ip, i, j] = 4.0
+            end
+        end
+        return nothing
+    end
+
+    @parallel (JustRelax.@idx ni) init_phases!(phases, particles.coords..., particles.index, r, Lx)
+end
diff --git a/miniapps/convection/WENO5/WENO_convection2D.jl b/miniapps/convection/WENO5/WENO_convection2D.jl
index 74f2073e..895af432 100644
--- a/miniapps/convection/WENO5/WENO_convection2D.jl
+++ b/miniapps/convection/WENO5/WENO_convection2D.jl
@@ -120,7 +120,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     r_anomaly        = 25e3   # radius of perturbation
     init_phases!(pPhases, particles, lx; d=abs(yc_anomaly), r=r_anomaly)
     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, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -261,7 +261,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
         @show it += 1
         t        += dt
 
diff --git a/src/common.jl b/src/common.jl
index 8c4d3d80..daf3244d 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -35,7 +35,7 @@ export FlowBoundaryConditions,
 include("MiniKernels.jl")
 
 include("phases/phases.jl")
-export fn_ratio, phase_ratios_center
+export fn_ratio, phase_ratios_center!
 
 include("rheology/BuoyancyForces.jl")
 export compute_ρg!
diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl
index dcab6325..92743b0f 100644
--- a/src/ext/AMDGPU/2D.jl
+++ b/src/ext/AMDGPU/2D.jl
@@ -92,14 +92,14 @@ function thermal_bcs!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays, bc
 end
 
 # Phases
-function JR2D.phase_ratios_center(
+function JR2D.phase_ratios_center!(
     ::AMDGPUBackendTrait,
     phase_ratios::JustRelax.PhaseRatio,
     particles,
     grid::Geometry,
     phases,
 )
-    return _phase_ratios_center(phase_ratios, particles, grid, phases)
+    return _phase_ratios_center!(phase_ratios, particles, grid, phases)
 end
 
 # Rheology
@@ -183,4 +183,42 @@ function JR2D.compute_dt(::AMDGPUBackendTrait, S::JustRelax.StokesArrays, args..
     return _compute_dt(S, args...)
 end
 
+# Subgrid diffusion
+
+function JR2D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::JustRelax.PhaseRatio,
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+)
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases.center, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
+function JR2D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::AbstractArray{Int,N},
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+) where {N}
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
 end
diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl
index 472905e1..acf39f5a 100644
--- a/src/ext/AMDGPU/3D.jl
+++ b/src/ext/AMDGPU/3D.jl
@@ -92,14 +92,14 @@ function thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs)
 end
 
 # Phases
-function JR3D.phase_ratios_center(
+function JR3D.phase_ratios_center!(
     ::CUDABackendTrait,
     phase_ratios::JustRelax.PhaseRatio,
     particles,
     grid::Geometry,
     phases,
 )
-    return _phase_ratios_center(phase_ratios, particles, grid, phases)
+    return _phase_ratios_center!(phase_ratios, particles, grid, phases)
 end
 
 # Rheology
@@ -183,4 +183,40 @@ function JR3D.compute_dt(::CUDABackendTrait, S::JustRelax.StokesArrays, args...)
     return _compute_dt(S, args...)
 end
 
+function JR3D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::JustRelax.PhaseRatio,
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+)
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases.center, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
+function JR3D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::AbstractArray{Int,N},
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+) where {N}
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
 end
diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl
index e4701169..74deef2d 100644
--- a/src/ext/CUDA/2D.jl
+++ b/src/ext/CUDA/2D.jl
@@ -92,14 +92,14 @@ function thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs)
 end
 
 # Phases
-function JR2D.phase_ratios_center(
+function JR2D.phase_ratios_center!(
     ::CUDABackendTrait,
     phase_ratios::JustRelax.PhaseRatio,
     particles,
     grid::Geometry,
     phases,
 )
-    return _phase_ratios_center(phase_ratios, particles, grid, phases)
+    return _phase_ratios_center!(phase_ratios, particles, grid, phases)
 end
 
 # Rheology
@@ -182,4 +182,42 @@ function JR2D.compute_dt(::CUDABackendTrait, S::JustRelax.StokesArrays, args...)
     return _compute_dt(S, args...)
 end
 
+# Subgrid diffusion
+
+function JR2D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::CuArray,
+    phases::JustRelax.PhaseRatio,
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+)
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases.center, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
+function JR2D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::CuArray,
+    phases::AbstractArray{Int,N},
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+) where {N}
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
 end
diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl
index 73158c3a..acf39f5a 100644
--- a/src/ext/CUDA/3D.jl
+++ b/src/ext/CUDA/3D.jl
@@ -92,43 +92,44 @@ function thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs)
 end
 
 # Phases
-function JR3D.phase_ratios_center(
+function JR3D.phase_ratios_center!(
     ::CUDABackendTrait,
     phase_ratios::JustRelax.PhaseRatio,
     particles,
     grid::Geometry,
     phases,
 )
-    return _phase_ratios_center(phase_ratios, particles, grid, phases)
+    return _phase_ratios_center!(phase_ratios, particles, grid, phases)
 end
 
 # Rheology
+
 ## viscosity
-function JR3D.compute_viscosity!(::CUDABackendTrait, stokes, ν, args, rheology, cutoff)
+function JR3D.compute_viscosity!(::AMDGPUBackendTrait, stokes, ν, args, rheology, cutoff)
     return _compute_viscosity!(stokes, ν, args, rheology, cutoff)
 end
 
 function JR3D.compute_viscosity!(
-    ::CUDABackendTrait, stokes, ν, phase_ratios, args, rheology, cutoff
+    ::AMDGPUBackendTrait, stokes, ν, phase_ratios, args, rheology, cutoff
 )
     return _compute_viscosity!(stokes, ν, phase_ratios, args, rheology, cutoff)
 end
 
-function JR3D.compute_viscosity!(η, ν, εII::CuArray, args, rheology, cutoff)
+function JR2D.compute_viscosity!(η, ν, εII::RocArray, args, rheology, cutoff)
     return compute_viscosity!(η, ν, εII, args, rheology, cutoff)
 end
 
-function compute_viscosity!(::CUDABackendTrait, stokes, ν, args, rheology, cutoff)
+function compute_viscosity!(::AMDGPUBackendTrait, stokes, ν, args, rheology, cutoff)
     return _compute_viscosity!(stokes, ν, args, rheology, cutoff)
 end
 
 function compute_viscosity!(
-    ::CUDABackendTrait, stokes, ν, phase_ratios, args, rheology, cutoff
+    ::AMDGPUBackendTrait, stokes, ν, phase_ratios, args, rheology, cutoff
 )
     return _compute_viscosity!(stokes, ν, phase_ratios, args, rheology, cutoff)
 end
 
-function compute_viscosity!(η, ν, εII::CuArray, args, rheology, cutoff)
+function compute_viscosity!(η, ν, εII::RocArray, args, rheology, cutoff)
     return compute_viscosity!(η, ν, εII, args, rheology, cutoff)
 end
 
@@ -138,10 +139,10 @@ function JR3D.tensor_invariant!(::CUDABackendTrait, A::JustRelax.SymmetricTensor
 end
 
 ## Buoyancy forces
-function JR3D.compute_ρg!(ρg::CuArray, rheology, args)
+function JR3D.compute_ρg!(ρg::RocArray, rheology, args)
     return compute_ρg!(ρg, rheology, args)
 end
-function JR3D.compute_ρg!(ρg::CuArray, phase_ratios::JustRelax.PhaseRatio, rheology, args)
+function JR3D.compute_ρg!(ρg::RocArray, phase_ratios::JustRelax.PhaseRatio, rheology, args)
     return compute_ρg!(ρg, phase_ratios, rheology, args)
 end
 
@@ -154,17 +155,17 @@ function temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArray
     return _temperature2center!(thermal)
 end
 
-function JR3D.vertex2center!(center::T, vertex::T) where {T<:CuArray}
+function JR3D.vertex2center!(center::T, vertex::T) where {T<:RocArray}
     return vertex2center!(center, vertex)
 end
 
-function JR3D.center2vertex!(vertex::T, center::T) where {T<:CuArray}
+function JR3D.center2vertex!(vertex::T, center::T) where {T<:RocArray}
     return center2vertex!(vertex, center)
 end
 
 function JR3D.center2vertex!(
     vertex_yz::T, vertex_xz::T, vertex_xy::T, center_yz::T, center_xz::T, center_xy::T
-) where {T<:CuArray}
+) where {T<:RocArray}
     return center2vertex!(vertex_yz, vertex_xz, vertex_xy, center_yz, center_xz, center_xy)
 end
 
@@ -182,4 +183,40 @@ function JR3D.compute_dt(::CUDABackendTrait, S::JustRelax.StokesArrays, args...)
     return _compute_dt(S, args...)
 end
 
+function JR3D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::JustRelax.PhaseRatio,
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+)
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases.center, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
+function JR3D.subgrid_characteristic_time!(
+    subgrid_arrays,
+    particles,
+    dt₀::RocArray,
+    phases::AbstractArray{Int,N},
+    rheology,
+    thermal::JustRelax.ThermalArrays,
+    stokes::JustRelax.StokesArrays,
+    xci,
+    di,
+) where {N}
+    ni = size(stokes.P)
+    @parallel (@idx ni) subgrid_characteristic_time!(
+        dt₀, phases, rheology, thermal.Tc, stokes.P, di
+    )
+    return nothing
+end
+
 end
diff --git a/src/phases/phases.jl b/src/phases/phases.jl
index 22452832..a0d96ac6 100644
--- a/src/phases/phases.jl
+++ b/src/phases/phases.jl
@@ -48,11 +48,13 @@ end
     end
 end
 
-function phase_ratios_center(phase_ratios, particles, grid, phases)
-    return phase_ratios_center(backend(phase_ratios), phase_ratios, particles, grid, phases)
+function phase_ratios_center!(phase_ratios, particles, grid, phases)
+    return phase_ratios_center!(
+        backend(phase_ratios), phase_ratios, particles, grid, phases
+    )
 end
 
-function phase_ratios_center(
+function phase_ratios_center!(
     ::CPUBackendTrait, phase_ratios::JustRelax.PhaseRatio, particles, grid::Geometry, phases
 )
     return _phase_ratios_center(phase_ratios, particles, grid, phases)
@@ -89,50 +91,24 @@ end
 function phase_ratio_weights(
     pxi::NTuple{NP,C}, ph::SVector{N1,T}, cell_center, di, ::Val{NC}
 ) where {N1,NC,NP,T,C}
-    if @generated
-        quote
-            Base.@_inline_meta
-            # Initiaze phase ratio weights (note: can't use ntuple() here because of the @generated function)
-            Base.@nexprs $NC i -> w_i = zero($T)
-            w = Base.@ncall $NC tuple w
-
-            # initialie sum of weights
-            sumw = zero($T)
-            Base.@nexprs $N1 i -> begin
-                # bilinear weight (1-(xᵢ-xc)/dx)*(1-(yᵢ-yc)/dy)
-                x = bilinear_weight(cell_center, getindex.(pxi, i), di)
-                sumw += x # reduce
-                ph_local = ph[i]
-                # this is doing sum(w * δij(i, phase)), where δij is the Kronecker delta
-                # Base.@nexprs $NC j -> tmp_j = w[j] + x * δ(Int(ph_local), j)
-                Base.@nexprs $NC j -> tmp_j = w[j] + x * (ph_local == j)
-                w = Base.@ncall $NC tuple tmp
-            end
-
-            # return phase ratios weights w = sum(w * δij(i, phase)) / sum(w)
-            _sumw = inv(sum(w))
-            Base.@nexprs $NC i -> w_i = w[i] * _sumw
-            w = Base.@ncall $NC tuple w
-            return w
-        end
-    else
-        # Initiaze phase ratio weights (note: can't use ntuple() here because of the @generated function)
-        w = ntuple(_ -> zero(T), Val(NC))
-        # initialie sum of weights
-        sumw = zero(T)
-
-        for i in eachindex(pxi)
-            # bilinear weight (1-(xᵢ-xc)/dx)*(1-(yᵢ-yc)/dy)
-            x = @inline bilinear_weight(cell_center, getindex.(pxi, i), di)
-            sumw += x # reduce
-            ph_local = ph[i]
-            # this is doing sum(w * δij(i, phase)), where δij is the Kronecker delta
-            # w = w .+ x .* ntuple(j -> δ(Int(ph_local), j), Val(NC))
-            w = w .+ x .* ntuple(j -> (ph_local == j), Val(NC))
-        end
-        w = w .* inv(sum(w))
-        return w
+
+    # Initiaze phase ratio weights (note: can't use ntuple() here because of the @generated function)
+    w = ntuple(_ -> zero(T), Val(NC))
+    sumw = zero(T)
+
+    for i in eachindex(ph)
+        # bilinear weight (1-(xᵢ-xc)/dx)*(1-(yᵢ-yc)/dy)
+        p = getindex.(pxi, i)
+        isnan(first(p)) && continue
+        x = @inline bilinear_weight(cell_center, p, di)
+        sumw += x # reduce
+        ph_local = ph[i]
+        # this is doing sum(w * δij(i, phase)), where δij is the Kronecker delta
+        # w = w .+ x .* ntuple(j -> δ(Int(ph_local), j), Val(NC))
+        w = w .+ x .* ntuple(j -> (ph_local == j), Val(NC))
     end
+    w = w .* inv(sumw)
+    return w
 end
 
 @generated function bilinear_weight(
diff --git a/test/test_VanKeken.jl b/test/test_VanKeken.jl
index 48d1daaa..bf977d2d 100644
--- a/test/test_VanKeken.jl
+++ b/test/test_VanKeken.jl
@@ -94,7 +94,7 @@ function VanKeken2D(ny=32, nx=32)
     particle_args        = (pPhases, )
     phase_ratios         = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem
@@ -171,7 +171,7 @@ function VanKeken2D(ny=32, nx=32)
         # inject && break
         inject_particles_phase!(particles, pPhases, (), (), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/test/test_diffusion2D_multiphase.jl b/test/test_diffusion2D_multiphase.jl
index e02ca525..2ac57db9 100644
--- a/test/test_diffusion2D_multiphase.jl
+++ b/test/test_diffusion2D_multiphase.jl
@@ -138,7 +138,7 @@ function diffusion_2D(; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3.0)
     pPhases,     = init_cell_arrays(particles, Val(1))
     init_phases!(pPhases, particles, center_perturbation..., r)
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     @parallel (@idx ni) compute_temperature_source_terms!(thermal.H, rheology, phase_ratios.center, args)
diff --git a/test/test_diffusion3D_multiphase.jl b/test/test_diffusion3D_multiphase.jl
index 6d68a9e2..6cd4e038 100644
--- a/test/test_diffusion3D_multiphase.jl
+++ b/test/test_diffusion3D_multiphase.jl
@@ -150,7 +150,7 @@ function diffusion_3D(;
     pPhases,     = init_cell_arrays(particles, Val(1))
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, center_perturbation..., r)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # PT coefficients for thermal diffusion
diff --git a/test/test_shearheating2D.jl b/test/test_shearheating2D.jl
index 23516ba2..06c757f0 100644
--- a/test/test_shearheating2D.jl
+++ b/test/test_shearheating2D.jl
@@ -88,7 +88,7 @@ function Shearheating2D()
     r_anomaly        = 3e3    # radius of perturbation
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, xc_anomaly, yc_anomaly, r_anomaly)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -214,7 +214,7 @@ function Shearheating2D()
         # check if we need to inject particles
         inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
         # update phase ratios
-        phase_ratios_center(phase_ratios, particles, grid, pPhases)
+        phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
         @show it += 1
         t        += dt
diff --git a/test/test_shearheating3D.jl b/test/test_shearheating3D.jl
index c6ed8764..ef88d44f 100644
--- a/test/test_shearheating3D.jl
+++ b/test/test_shearheating3D.jl
@@ -82,7 +82,7 @@ function Shearheating3D(nx=16, ny=16, nz=16)
     r_anomaly        = 3e3    # radius of perturbation
     phase_ratios     = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, xc_anomaly, yc_anomaly, zc_anomaly, r_anomaly)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
     # ----------------------------------------------------
 
     # STOKES ---------------------------------------------
@@ -201,7 +201,7 @@ function Shearheating3D(nx=16, ny=16, nz=16)
             # check if we need to inject particles
             inject_particles_phase!(particles, pPhases, (pT, ), (thermal.T,), xvi)
             # update phase ratios
-            phase_ratios_center(phase_ratios, particles, grid, pPhases)
+            phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
             @show it += 1
             t        += dt
diff --git a/test/test_sinking_block.jl b/test/test_sinking_block.jl
index b960436b..927060f1 100644
--- a/test/test_sinking_block.jl
+++ b/test/test_sinking_block.jl
@@ -124,7 +124,7 @@ function Sinking_Block2D()
     r_anomaly    =  50e3   # radius of perturbation
     phase_ratios = PhaseRatio(backend_JR, ni, length(rheology))
     init_phases!(pPhases, particles, xc_anomaly, abs(yc_anomaly), r_anomaly)
-    phase_ratios_center(phase_ratios, particles, grid, pPhases)
+    phase_ratios_center!(phase_ratios, particles, grid, pPhases)
 
     # STOKES ---------------------------------------------
     # Allocate arrays needed for every Stokes problem