diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index bb45dff3..9a1a0ab9 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -6,6 +6,12 @@ on:
   pull_request:
     branches:
       - main
+
+# Cancel redundant CI tests automatically
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
diff --git a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
index ee7dade3..ae2a9ff0 100755
--- a/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
+++ b/miniapps/benchmarks/stokes2D/Blankenbach2D/Benchmark2D_WENO5.jl
@@ -73,7 +73,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 1e1, figdir="figs2D", do_vtk =fal
     # ----------------------------------------------------
 
     # Weno model -----------------------------------------
-    weno = WENO5(Val(2), (nx,ny).+1) # ni.+1 for Temp
+    weno = WENO5(backend, Val(2), (nx,ny).+1) # ni.+1 for Temp
     # ----------------------------------------------------
 
     # Initialize particles -------------------------------
diff --git a/miniapps/convection/GlobalConvection2D_WENO5.jl b/miniapps/convection/GlobalConvection2D_WENO5.jl
index 76c136c5..d8a60d65 100644
--- a/miniapps/convection/GlobalConvection2D_WENO5.jl
+++ b/miniapps/convection/GlobalConvection2D_WENO5.jl
@@ -94,7 +94,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", therma
     # ----------------------------------------------------
 
     # Weno model -----------------------------------------
-    weno = WENO5(Val(2), ni.+1) # ni.+1 for Temp
+    weno = WENO5(backend, Val(2), ni.+1) # ni.+1 for Temp
     # ----------------------------------------------------
 
     # create rheology struct
diff --git a/miniapps/convection/WENO5/WENO_convection2D.jl b/miniapps/convection/WENO5/WENO_convection2D.jl
index 283b4eb5..6e57a614 100644
--- a/miniapps/convection/WENO5/WENO_convection2D.jl
+++ b/miniapps/convection/WENO5/WENO_convection2D.jl
@@ -101,7 +101,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
     # Weno model -----------------------------------------
     weno = WENO5(backend_JR, Val(2), ni.+1)
     # ----------------------------------------------------
-    
+
     # Initialize particles -------------------------------
     nxcell, max_xcell, min_xcell = 20, 40, 1
     particles        = init_particles(backend, nxcell, max_xcell, min_xcell, xvi...);
@@ -285,8 +285,8 @@ end
 figdir   = "Weno2D"
 ar       = 1 # aspect ratio
 n        = 64
-nx       = n*ar 
-ny       = n 
+nx       = n*ar
+ny       = n
 igg      = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
     IGG(init_global_grid(nx, ny, 1; init_MPI= true)...)
 else
diff --git a/src/ext/AMDGPU/2D.jl b/src/ext/AMDGPU/2D.jl
index 8a99e040..0b5b9804 100644
--- a/src/ext/AMDGPU/2D.jl
+++ b/src/ext/AMDGPU/2D.jl
@@ -364,4 +364,8 @@ function JR2D.compute_shear_heating!(
     return nothing
 end
 
+function JR2D.WENO_advection!(u::ROCArrayArray, Vxi::NTuple, weno, di, dt)
+    return WENO_advection!(u, Vxi, weno, di, dt)
+end
+
 end
diff --git a/src/ext/AMDGPU/3D.jl b/src/ext/AMDGPU/3D.jl
index 00b53a2a..a40c754c 100644
--- a/src/ext/AMDGPU/3D.jl
+++ b/src/ext/AMDGPU/3D.jl
@@ -236,10 +236,6 @@ function JR3D.compute_melt_fraction!(
     return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
 end
 
-function JR3D.compute_melt_fraction!(ϕ::ROCArray, phase_ratios, rheology, args)
-    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
-end
-
 # Interpolations
 function JR3D.temperature2center!(::AMDGPUBackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)
@@ -385,4 +381,9 @@ function JR3D.compute_shear_heating!(
     )
     return nothing
 end
+
+function JR3D.WENO_advection!(u::ROCArray, Vxi::NTuple, weno, di, dt)
+    return WENO_advection!(u, Vxi, weno, di, dt)
+end
+
 end
diff --git a/src/ext/CUDA/2D.jl b/src/ext/CUDA/2D.jl
index cd6a1602..5f448607 100644
--- a/src/ext/CUDA/2D.jl
+++ b/src/ext/CUDA/2D.jl
@@ -348,4 +348,9 @@ function JR2D.compute_shear_heating!(
     return nothing
 end
 
+# WENO-5 advection
+function JR2D.WENO_advection!(u::CuArray, Vxi::NTuple, weno, di, dt)
+    return WENO_advection!(u, Vxi, weno, di, dt)
+end
+
 end
diff --git a/src/ext/CUDA/3D.jl b/src/ext/CUDA/3D.jl
index aa5737db..da702989 100644
--- a/src/ext/CUDA/3D.jl
+++ b/src/ext/CUDA/3D.jl
@@ -180,20 +180,6 @@ end
 function thermal_bcs!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays, bcs)
     return thermal_bcs!(thermal.T, bcs)
 end
-
-# # Phases
-# function JR3D.phase_ratios_center!(
-#     ::CUDABackendTrait, phase_ratios::JustPIC.PhaseRatios, particles, grid::Geometry, phases
-# )
-#     return _phase_ratios_center!(phase_ratios, particles, grid, phases)
-# end
-
-# function JR3D.phase_ratios_vertex!(
-#     ::CUDABackendTrait, phase_ratios::JustPIC.PhaseRatios, particles, grid::Geometry, phases
-# )
-#     return _phase_ratios_vertex!(phase_ratios, particles, grid, phases)
-# end
-
 # Rheology
 
 ## viscosity
@@ -252,10 +238,6 @@ function JR3D.compute_melt_fraction!(
     return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
 end
 
-function JR3D.compute_melt_fraction!(ϕ::CuArray, phase_ratios, rheology, args)
-    return compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
-end
-
 # Interpolations
 function JR3D.temperature2center!(::CUDABackendTrait, thermal::JustRelax.ThermalArrays)
     return _temperature2center!(thermal)
@@ -399,4 +381,8 @@ function JR3D.compute_shear_heating!(
     return nothing
 end
 
+function JR3D.WENO_advection!(u::CuArray, Vxi::NTuple, weno, di, dt)
+    return WENO_advection!(u, Vxi, weno, di, dt)
+end
+
 end
diff --git a/src/rheology/Melting.jl b/src/rheology/Melting.jl
index 7a340641..04ed5c82 100644
--- a/src/rheology/Melting.jl
+++ b/src/rheology/Melting.jl
@@ -1,31 +1,26 @@
 function compute_melt_fraction!(ϕ, rheology, args)
     ni = size(ϕ)
     @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, rheology, args)
+    return nothing
 end
 
 @parallel_indices (I...) function compute_melt_fraction_kernel!(ϕ, rheology, args)
     args_ijk = ntuple_idx(args, I...)
-    ϕ[I...] = compute_melt_frac(rheology, args_ijk)
+    ϕ[I...] = compute_meltfraction(rheology, args_ijk)
     return nothing
 end
 
-@inline function compute_melt_frac(rheology, args)
-    return compute_meltfraction(rheology, args)
-end
-
-function compute_melt_fraction!(ϕ, phase_ratios, rheology, args)
+function compute_melt_fraction!(ϕ, phase_ratios::JustPIC.PhaseRatios, rheology, args)
     ni = size(ϕ)
-    @parallel (@idx ni) compute_melt_fraction_kernel!(ϕ, phase_ratios, rheology, args)
+    @parallel (@idx ni) compute_melt_fraction_kernel!(
+        ϕ, phase_ratios.center, rheology, args
+    )
 end
 
 @parallel_indices (I...) function compute_melt_fraction_kernel!(
     ϕ, phase_ratios, rheology, args
 )
     args_ijk = ntuple_idx(args, I...)
-    ϕ[I...] = compute_melt_frac(rheology, args_ijk, phase_ratios[I...])
+    ϕ[I...] = fn_ratio(compute_meltfraction, rheology, @cell(phase_ratios[I...]), args_ijk)
     return nothing
 end
-
-@inline function compute_melt_frac(rheology, args, phase_ratios)
-    return compute_meltfraction_ratio(phase_ratios, rheology, args)
-end
diff --git a/src/types/weno.jl b/src/types/weno.jl
index f980d47d..1c31bd89 100644
--- a/src/types/weno.jl
+++ b/src/types/weno.jl
@@ -1,3 +1,4 @@
+using Adapt
 ## Weno5 advection scheme. Implementation based on the repository from
 # https://gmd.copernicus.org/preprints/gmd-2023-189/
 
@@ -59,3 +60,5 @@ struct WENO5{T,N,A,M} <: AbstractWENO
     # method
     method::M # 1:JS, 2:Z
 end
+
+Adapt.@adapt_structure WENO5
diff --git a/test/test_Blankenbach.jl b/test/test_Blankenbach.jl
index fd5ded93..bcea11d4 100644
--- a/test/test_Blankenbach.jl
+++ b/test/test_Blankenbach.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using Test, Suppressor
@@ -137,12 +135,6 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 10)
     temperature2center!(thermal)
     # ----------------------------------------------------
 
-    # Rayleigh number
-    ΔT = thermal.T[1,1] - thermal.T[1,end]
-    Ra = (rheology[1].Density[1].ρ0 * rheology[1].Gravity[1].g * rheology[1].Density[1].α * ΔT * ly^3.0 ) /
-       (κ * rheology[1].CompositeRheology[1].elements[1].η )
-    @show Ra
-
     args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
 
     # Buoyancy forces  & viscosity ----------------------
@@ -264,7 +256,7 @@ function main2D(igg; ar=1, nx=32, ny=32, nit = 10)
         # Compute U rms -----------------------------
         # U₍ᵣₘₛ₎ = H*ρ₀*c₍ₚ₎/k * √ 1/H/L * ∫∫ (vx²+vz²) dx dz
         Urms_it = let
-            JustRelax.JustRelax2D.velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=true)
+            velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=true)
             @. Vx_v .= hypot.(Vx_v, Vy_v) # we reuse Vx_v to store the velocity magnitude
             sqrt( sum( Vx_v.^2 .* prod(di)) / lx /ly ) *
                 ((ly * rheology[1].Density[1].ρ0 * rheology[1].HeatCapacity[1].Cp) / rheology[1].Conductivity[1].k )
diff --git a/test/test_Interpolations.jl b/test/test_Interpolations.jl
index 08a2fca8..3bc0b317 100644
--- a/test/test_Interpolations.jl
+++ b/test/test_Interpolations.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using Test
@@ -27,7 +25,7 @@ end
 
 
 @testset "Interpolations" begin
-
+    if backend_JR == CPUBackend
         # Set up mock data
         # Physical domain ------------------------------------
         ly           = 1.0       # domain length in y
@@ -75,4 +73,7 @@ end
         velocity2vertex!(Vx_v, Vy_v, stokes.V.Vx, stokes.V.Vy; ghost_nodes=false)
         @test iszero(Vx_v[1,1])
         @test Vy_v[1,1] == 10
+    else
+        @test true == true
+    end
 end
diff --git a/test/test_WENO5.jl b/test/test_WENO5.jl
index ddd2b008..ccd7c2d4 100644
--- a/test/test_WENO5.jl
+++ b/test/test_WENO5.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using Test, Suppressor
@@ -112,7 +110,7 @@ function thermal_convection2D(igg; ar=8, ny=16, nx=ny*8, thermal_perturbation =
     # ----------------------------------------------------
 
     # Weno model -----------------------------------------
-    weno = WENO5(Val(2), ni.+1) # ni.+1 for Temp
+    weno = WENO5(backend_JR, Val(2), ni.+1) # ni.+1 for Temp
     # ----------------------------------------------------
 
     # create rheology struct
diff --git a/test/test_boundary_conditions2D.jl b/test/test_boundary_conditions2D.jl
index dbb0d6b6..285c8e60 100644
--- a/test/test_boundary_conditions2D.jl
+++ b/test/test_boundary_conditions2D.jl
@@ -1,9 +1,7 @@
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using JustRelax, JustRelax.JustRelax2D
diff --git a/test/test_boundary_conditions3D.jl b/test/test_boundary_conditions3D.jl
index 88aa96d6..2630dfc9 100644
--- a/test/test_boundary_conditions3D.jl
+++ b/test/test_boundary_conditions3D.jl
@@ -1,9 +1,7 @@
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using JustRelax, JustRelax.JustRelax3D
diff --git a/test/test_shearheating2D.jl b/test/test_shearheating2D.jl
index b3335052..e71f5578 100644
--- a/test/test_shearheating2D.jl
+++ b/test/test_shearheating2D.jl
@@ -3,10 +3,8 @@ using Test, Suppressor
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 # Benchmark of Duretz et al. 2014
diff --git a/test/test_shearheating3D.jl b/test/test_shearheating3D.jl
index e79cf381..68b29d58 100644
--- a/test/test_shearheating3D.jl
+++ b/test/test_shearheating3D.jl
@@ -3,10 +3,8 @@ using Test, Suppressor
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 # Benchmark of Duretz et al. 2014
diff --git a/test/test_sinking_block.jl b/test/test_sinking_block.jl
index 3f9f9be6..50a823fa 100644
--- a/test/test_sinking_block.jl
+++ b/test/test_sinking_block.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true) #used because of velocity2vertex!
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true) #used because of velocity2vertex!
 end
 
 
diff --git a/test/test_stokes_solvi3D.jl b/test/test_stokes_solvi3D.jl
index 5e24415b..ea227fa8 100644
--- a/test/test_stokes_solvi3D.jl
+++ b/test/test_stokes_solvi3D.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using Test, Suppressor
diff --git a/test/test_thermalstresses.jl b/test/test_thermalstresses.jl
index ca538c7e..9dc806e4 100644
--- a/test/test_thermalstresses.jl
+++ b/test/test_thermalstresses.jl
@@ -2,10 +2,8 @@ push!(LOAD_PATH, "..")
 
 @static if ENV["JULIA_JUSTRELAX_BACKEND"] === "AMDGPU"
     using AMDGPU
-    AMDGPU.allowscalar(true)
 elseif ENV["JULIA_JUSTRELAX_BACKEND"] === "CUDA"
     using CUDA
-    CUDA.allowscalar(true)
 end
 
 using Test, Suppressor
@@ -222,10 +220,7 @@ function init_rheology(CharDim; is_compressible = false, steady_state=true)
 end
 
 
-function main2D(; nx=32, ny=32)
-
-    init_mpi = JustRelax.MPI.Initialized() ? false : true
-    igg      = IGG(init_global_grid(nx, ny, 1; init_MPI = init_mpi)...)
+function main2D(igg; nx=32, ny=32)
 
     # Characteristic lengths
     CharDim      = GEO_units(;length=12.5km, viscosity=1e21, temperature = 1e3C)
@@ -317,7 +312,7 @@ function main2D(; nx=32, ny=32)
 
     ϕ = @zeros(ni...)
     compute_melt_fraction!(
-        ϕ, phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P)
+        ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P)
     )
     # 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
@@ -343,8 +338,8 @@ function main2D(; nx=32, ny=32)
     @copy stokes.P0 stokes.P
     thermal.Told .= thermal.T
     P_init        = deepcopy(stokes.P)
-    Tsurf         = thermal.T[1, end]
-    Tbot          = thermal.T[1, 1]
+    Tsurf         = nondimensionalize(273K,CharDim)
+    Tbot          = nondimensionalize(648K,CharDim)
     local ϕ, stokes, thermal
 
     while it < 1
@@ -416,7 +411,7 @@ function main2D(; nx=32, ny=32)
         )
         # ------------------------------
         compute_melt_fraction!(
-            ϕ, phase_ratios.center, rheology, (T=thermal.Tc, P=stokes.P)
+            ϕ, phase_ratios, rheology, (T=thermal.Tc, P=stokes.P)
         )
         # Advection --------------------
         # advect particles in space
@@ -448,9 +443,17 @@ end
 
 @testset "thermal stresses" begin
     @suppress begin
-        ϕ, stokes, thermal = main2D(; nx=32, ny=32)
-        nx_T, ny_T = size(thermal.T) .>>> 1
-        @test Array(thermal.T)[nx_T  + 1, ny_T + 1] ≈ 0.5369 rtol = 1e-2
-        @test Array(ϕ)[nx_T + 1, ny_T + 1] ≤ 1e-8 
+        nx, ny   = 32, 32           # number of cells
+        igg      = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
+            IGG(init_global_grid(nx, ny, 1; init_MPI= true)...)
+        else
+            igg
+        end
+
+        ϕ, stokes, thermal = main2D(igg; nx=32, ny=32)
+
+        nx_T, ny_T = size(thermal.T)
+        @test  Array(thermal.T)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 0.5369 rtol = 1e-2
+        @test  Array(ϕ)[nx_T >>> 1 + 1, ny_T >>> 1 + 1] ≈ 9.351e-9 rtol = 1e-1
     end
 end
\ No newline at end of file