diff --git a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
index 61a0afdb..f21a9dae 100644
--- a/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
+++ b/miniapps/convection/Particles2D_nonDim/Layered_convection2D.jl
@@ -204,7 +204,7 @@ function main2D(igg; ar=8, ny=16, nx=ny*8, figdir="figs2D", do_vtk =false)
         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)
-        scatter!(ax2, Array(log10.(A[:])), Y)
+        scatter!(ax2, Array(log10.(η[:])), Y)
         ylims!(ax1, minimum(xvi[2]), 0)
         ylims!(ax2, minimum(xvi[2]), 0)
         hideydecorations!(ax2)
@@ -395,4 +395,4 @@ else
 end
 
 # run main script
-main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
+# main2D(igg; figdir = figdir, ar = ar, nx = nx, ny = ny, do_vtk = do_vtk);
diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl
index 3be0b181..4e045f0c 100644
--- a/src/stokes/Stokes2D.jl
+++ b/src/stokes/Stokes2D.jl
@@ -556,6 +556,7 @@ function JustRelax.solve!(
 
         wtime0 += @elapsed begin
             compute_maxloc!(ητ, η; window=(1, 1))
+            # compute_maxloc!(ητ, η_vep; window=(1, 1))
             update_halo!(ητ)
 
             @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes)..., _di...)
@@ -628,16 +629,25 @@ function JustRelax.solve!(
             )
 
             @hide_communication b_width begin # communication/computation overlap
+                # @parallel compute_V!(
+                #     @velocity(stokes)...,
+                #     Vx_on_Vy,
+                #     θ,
+                #     @stress(stokes)...,
+                #     pt_stokes.ηdτ,
+                #     ρg...,
+                #     ητ,
+                #     _di...,
+                #     dt * free_surface,
+                # )
                 @parallel compute_V!(
                     @velocity(stokes)...,
-                    Vx_on_Vy,
-                    θ,
+                    stokes.P,
                     @stress(stokes)...,
                     pt_stokes.ηdτ,
                     ρg...,
                     ητ,
                     _di...,
-                    dt * free_surface,
                 )
                 # apply boundary conditions
                 free_surface_bcs!(stokes, flow_bcs, η, rheology, phase_ratios, dt, di)
diff --git a/subduction/LAMEM2D.jl b/subduction/LAMEM2D.jl
index afcde93f..89eec037 100644
--- a/subduction/LAMEM2D.jl
+++ b/subduction/LAMEM2D.jl
@@ -1,21 +1,21 @@
-using CUDA
+# using CUDA
 using JustRelax, JustRelax.DataIO
 import JustRelax.@cell
 using ParallelStencil
-# @init_parallel_stencil(Threads, Float64, 2)
-@init_parallel_stencil(CUDA, Float64, 2)
+@init_parallel_stencil(Threads, Float64, 2)
+# @init_parallel_stencil(CUDA, 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
-const backend = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+const backend = CPUBackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+# const backend = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
 
 # setup ParallelStencil.jl environment
-# model = PS_Setup(:cpu, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
-model = PS_Setup(:CUDA, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
+model = PS_Setup(:cpu, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
+# model = PS_Setup(:CUDA, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
 environment!(model)
 
 # Load script dependencies
@@ -93,8 +93,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
     # ----------------------------------------------------
 
     # TEMPERATURE PROFILE --------------------------------
-    Ttop             = 20
-    Tbot             = 1565.0
+    Ttop             = 20 + 273
+    Tbot             = 1565.0 + 273
     thermal          = ThermalArrays(ni)
     @views thermal.T[2:end-1, :] .= PTArray(T_GMG)
     thermal_bc       = TemperatureBoundaryConditions(;
@@ -134,22 +134,22 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
     flow_bcs!(stokes, flow_bcs) # apply boundary conditions
     update_halo!(@velocity(stokes)...)
 
-    # # 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)
-    #     # scatter!(ax2, Array(ρg[2][:]), 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
+    # 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)
+        # scatter!(ax2, Array(ρg[2][:]), 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
 
     # IO -------------------------------------------------
     # if it does not exist, make folder where figures are stored
@@ -213,8 +213,8 @@ function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk
                 args,
                 dt,
                 igg;
-                iterMax          = 150e3,
-                nout             = 2e3,
+                iterMax          = 100e3,
+                nout             = 1e3,
                 viscosity_cutoff = (1e18, 1e24),
                 free_surface     = false,
                 # viscosity_relaxation = 1e-5
@@ -354,5 +354,5 @@ else
     igg
 end
 
-main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
+# main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
 
diff --git a/subduction/LAMEM2D_nondim.jl b/subduction/LAMEM2D_nondim.jl
new file mode 100644
index 00000000..d43357c9
--- /dev/null
+++ b/subduction/LAMEM2D_nondim.jl
@@ -0,0 +1,368 @@
+# using CUDA
+using JustRelax, JustRelax.DataIO
+import JustRelax.@cell
+using ParallelStencil
+@init_parallel_stencil(Threads, Float64, 2)
+# @init_parallel_stencil(CUDA, 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
+# const backend = CUDABackend # Options: CPUBackend, CUDABackend, AMDGPUBackend
+
+# setup ParallelStencil.jl environment
+model = PS_Setup(:cpu, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
+# model = PS_Setup(:CUDA, Float64, 2) # or (:CUDA, Float64, 3) or (:AMDGPU, Float64, 3)
+environment!(model)
+
+# Load script dependencies
+using Printf, LinearAlgebra, GeoParams, GLMakie, CellArrays
+
+# Load file with all the rheology configurations
+include("LAMEM_rheology_nondim.jl")
+include("LAMEM_setup2D.jl")
+
+## SET OF HELPER FUNCTIONS PARTICULAR FOR THIS SCRIPT --------------------------------
+
+import ParallelStencil.INDICES
+const idx_k = INDICES[2]
+macro all_k(A)
+    esc(:($A[$idx_k]))
+end
+
+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
+
+
+# Initial pressure profile - not accurate
+@parallel function init_P!(P, ρg, z)
+    @all(P) = abs(@all(ρg) * @all_k(z)) * <(@all_k(z), 0.0)
+    return nothing
+end
+## END OF HELPER FUNCTION ------------------------------------------------------------
+
+## BEGIN OF MAIN SCRIPT --------------------------------------------------------------
+function main(li, origin, phases_GMG, igg; nx=16, ny=16, figdir="figs2D", do_vtk =false)
+
+    thickness    = 660 * km
+    η0           = 1e20Pa*s
+    CharDim      = GEO_units(; 
+        length = thickness, viscosity = η0, temperature = 1e3C
+    )
+
+    # Physical domain ------------------------------------
+    ni           = nx, ny           # number of cells
+    li_nd        = nondimensionalize(li[1]m, CharDim), nondimensionalize(li[2]m, CharDim)
+    origin_nd    = nondimensionalize(origin[1]m, CharDim), nondimensionalize(origin[2]m, CharDim)
+    di           = @. li_nd / ni       # grid steps
+    grid         = Geometry(ni, li_nd; origin = origin_nd)
+    (; xci, xvi) = grid # nodes at the center and vertices of the cells
+    # ----------------------------------------------------
+
+    # Physical properties using GeoParams ----------------
+    rheology     = init_rheologies(CharDim)
+    dt           = nondimensionalize(10e3 * 3600 * 24 * 365 * s, CharDim) # diffusive CFL timestep limiter
+    # ----------------------------------------------------
+
+    # Initialize particles -------------------------------
+    nxcell          = 40
+    max_xcell       = 60
+    min_xcell       = 20
+    particles       = init_particles(
+        backend, nxcell, max_xcell, min_xcell, xvi, di, ni
+    )
+    subgrid_arrays  = SubgridDiffusionCellArrays(particles)
+    # velocity grids
+    grid_vxi        = velocity_grids(xci, xvi, di)
+    # temperature
+    pPhases, pT     = init_cell_arrays(particles, Val(2))
+    particle_args   = (pPhases, pT)
+
+    # Assign particles phases anomaly
+    phases_device    = PTArray(phases_GMG)
+    init_phases!(pPhases, phases_device, particles, xvi)
+    phase_ratios     = PhaseRatio(ni, length(rheology))
+    phase_ratios_center!(phase_ratios, particles, xci, di, 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)
+    # ----------------------------------------------------
+
+    # TEMPERATURE PROFILE --------------------------------
+    Ttop             = nondimensionalize(20e0C, CharDim)
+    Tbot             = nondimensionalize(1565e0C, CharDim)
+    thermal          = ThermalArrays(ni)
+    @views thermal.T[2:end-1, :] .= PTArray(nondimensionalize(T_GMG .* K, CharDim))
+    thermal_bc       = TemperatureBoundaryConditions(;
+        no_flux      = (left = true, right = true, top = false, bot = false),
+    )
+    thermal_bcs!(thermal, thermal_bc)
+    # @views thermal.T[:, end] .= Ttop 
+    # @views thermal.T[:, 1]   .= Tbot 
+    @parallel (@idx ni) temperature2center!(thermal.Tc, thermal.T)
+    # ----------------------------------------------------
+
+    # Buoyancy forces
+    ρg               = ntuple(_ -> @zeros(ni...), Val(2))
+    # for _ in 1:2
+    #     compute_ρg!(ρg[2], phase_ratios, rheology, (T=thermal.Tc, P=stokes.P))
+    #     JustRelax.Stokes2D.init_P!(stokes.P, ρg[2], xci[2])
+    # end
+    
+    # Rheology
+    η                = @ones(ni...)
+    η_vep            = similar(η)
+    args             = (; T = thermal.Tc, P = stokes.P, dt = Inf)
+    viscosity_cutoff = nondimensionalize((1e16Pa*s, 1e24Pa*s), CharDim)
+    compute_viscosity!(
+        η, 1.0, phase_ratios, stokes, args, rheology, viscosity_cutoff 
+    )
+
+    # PT coefficients for thermal diffusion
+    pt_thermal       = PTThermalCoeffs(
+        rheology, phase_ratios, args, dt, ni, di, li; ϵ=1e-5, CFL=1e-2 / √3
+    )
+
+    # Boundary conditions
+    flow_bcs         = FlowBoundaryConditions(;
+        free_slip    = (left = true , right = true , top = true , bot = true),
+        free_surface = false,
+    )
+    flow_bcs!(stokes, flow_bcs) # apply boundary conditions
+    update_halo!(@velocity(stokes)...)
+
+    # 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)
+        # scatter!(ax2, Array(ρg[2][:]), 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
+
+    # 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)
+    # ----------------------------------------------------
+
+    local Vx_v, Vy_v
+    if do_vtk
+        Vx_v = @zeros(ni.+1...)
+        Vy_v = @zeros(ni.+1...)
+    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)
+
+    # Time loop
+    t, it = 0.0, 0
+
+    while it < 1000 # run only for 5 Myrs
+    # 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]      .= Ttop
+        # @views T_buffer[:, 1]        .= Tbot
+        @views thermal.T[2:end-1, :] .= T_buffer
+        thermal_bcs!(thermal, thermal_bc)
+        temperature2center!(thermal)
+        
+        # interpolate fields from particle to grid vertices
+        # particle2grid!(thermal.T, pT, xvi, particles)
+        # temperature2center!(thermal)
+        # Update buoyancy and viscosity -
+        args = (; T = thermal.Tc, P = stokes.P,  dt=Inf)
+        compute_viscosity!(
+            η, 1.0, phase_ratios, stokes, args, rheology, viscosity_cutoff
+        )
+        compute_ρg!(ρg[2], phase_ratios, rheology, args)
+        
+        # Stokes solver ----------------
+        t_stokes = @elapsed begin
+            out = solve!(
+                stokes,
+                pt_stokes,
+                di,
+                flow_bcs,
+                ρg,
+                η,
+                η_vep,
+                phase_ratios,
+                rheology,
+                args,
+                Inf,
+                igg;
+                iterMax          = 1,
+                nout             = 1,
+                viscosity_cutoff = viscosity_cutoff,
+                free_surface     = false,
+                viscosity_relaxation = 1e-3
+            );
+        end
+        println("Stokes solver time             ")
+        println("   Total time:      $t_stokes s")
+        println("   Time/iteration:  $(t_stokes / out.iter) s")
+        @parallel (JustRelax.@idx ni) JustRelax.Stokes2D.tensor_invariant!(stokes.ε.II, @strain(stokes)...)
+        @parallel (JustRelax.@idx ni) JustRelax.Stokes2D.tensor_invariant!(stokes.τ.II, @strain(stokes)...)
+        dt   = compute_dt(stokes, di)
+        # ------------------------------
+
+        # Thermal solver ---------------
+        heatdiffusion_PT!(
+            thermal,
+            pt_thermal,
+            thermal_bc,
+            rheology,
+            args,
+            dt,
+            di;
+            igg     = igg,
+            phase   = phase_ratios,
+            iterMax = 10e3,
+            nout    = 1e2,
+            verbose = true,
+        )
+        subgrid_characteristic_time!(
+            subgrid_arrays, particles, dt₀, phase_ratios, rheology, thermal, stokes, xci, di
+        )
+        centroid2particle!(subgrid_arrays.dt₀, xci, dt₀, particles)
+        subgrid_diffusion!(
+            pT, thermal.T, thermal.ΔT, subgrid_arrays, particles, xvi,  di, dt
+        )
+        # ------------------------------
+
+        # Advection --------------------
+        # advect particles in space
+        advection!(particles, RungeKutta2(), @velocity(stokes), grid_vxi, dt)
+        # advect particles in memory
+        move_particles!(particles, xvi, particle_args)
+        # check if we need to inject particles
+        inject_particles_phase!(particles, pPhases, (pT, ), (T_buffer,), xvi)
+        # update phase ratios
+        @parallel (@idx ni) phase_ratios_center(phase_ratios.center, particles.coords, xci, di, pPhases)
+
+        @show it += 1
+        t        += dt
+
+        # Data I/O and plotting ---------------------
+        if it == 1 || rem(it, 5) == 0
+            checkpointing(figdir, stokes, thermal.T, η, t)
+
+            if do_vtk
+                JustRelax.velocity2vertex!(Vx_v, Vy_v, @velocity(stokes)...)
+                data_v = (;
+                    T   = Array(thermal.T),
+                    τ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(η),
+                )
+                velocity_v = (
+                    Array(Vx_v),
+                    Array(Vy_v),
+                )
+                save_vtk(
+                    joinpath(vtk_dir, "vtk_" * lpad("$it", 6, "0")),
+                    xvi,
+                    xci,
+                    data_v,
+                    data_c,
+                    velocity_v
+                )
+            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
+            ar  = 3
+            fig = Figure(size = (1200, 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 ----------------------------------------------------------------
+do_vtk   = true # set to true to generate VTK files for ParaView
+figdir   = "Subduction_LAMEM_2D"
+# nx, ny   = 512, 256
+# nx, ny   = 512, 128
+n        = 64
+nx, ny   = n*6, n
+li, origin, phases_GMG, T_GMG = GMG_subduction_2D(nx+1, ny+1)
+igg      = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid
+    IGG(init_global_grid(nx, ny, 1; init_MPI= true)...)
+else
+    igg
+end
+
+# main(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, do_vtk = do_vtk);
diff --git a/subduction/LAMEM_rheology.jl b/subduction/LAMEM_rheology.jl
index d3d2c22f..de9df53a 100644
--- a/subduction/LAMEM_rheology.jl
+++ b/subduction/LAMEM_rheology.jl
@@ -71,15 +71,16 @@ function init_rheologies()
             Phase             = 3,
             Density           = PT_Density(; ρ0=3.3e3, α = α, β = 0e0, T0 = 273),
             HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
-            Conductivity      = ConstantConductivity(; k =3 ),
+            Conductivity      = ConstantConductivity(; k = 3),
             CompositeRheology = CompositeRheology( 
-                (
-                    disl_oceanic_crust,
-                    ConstantElasticity(; G=5e10, ν=0.5),
-                    DruckerPrager_regularised(; C = C_oceanic_litho, ϕ = ϕ_oceanic_litho, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity
-                ) 
-            ),
-            RadioactiveHeat   = ConstantRadioactiveHeat(2),
+                    (
+                        disl_dry_olivine, 
+                        diff_dry_olivine,
+                        ConstantElasticity(; G=5e10, ν=0.5),
+                        DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+                    ) 
+                ),
+            RadioactiveHeat   = ConstantRadioactiveHeat(6.6667e-12),
             Elasticity        = ConstantElasticity(; G=5e10, ν=0.5),
         ),
         # Name              = "continental crust",
diff --git a/subduction/LAMEM_rheology_nondim.jl b/subduction/LAMEM_rheology_nondim.jl
new file mode 100644
index 00000000..d01b83b6
--- /dev/null
+++ b/subduction/LAMEM_rheology_nondim.jl
@@ -0,0 +1,173 @@
+using GeoParams.Diffusion
+using GeoParams.Dislocation
+
+function init_rheologies(CharDim)
+    disl_dry_olivine   = SetDislocationCreep(Dislocation.dry_olivine_Hirth_2003; V = 14.5e-6m^3 / mol)
+    disl_oceanic_crust = SetDislocationCreep(Dislocation.plagioclase_An75_Ji_1993)
+    # disl_oceanic_litho = SetDislocationCreep(Dislocation.plagioclase_An75_Ji_1993)
+    disl_cont_crust    = SetDislocationCreep(Dislocation.wet_quartzite_Kirby_1983)
+
+    Transform_DislocationCreep(Dislocation.wet_quartzite_Kirby_1983)
+
+    diff_dry_olivine   = SetDiffusionCreep(Diffusion.dry_olivine_Hirth_2003; V = 14.5e-6m^3 / mol)
+
+    ϕ_dry_olivine   = sind(20)
+    C_dry_olivine   = 30e6Pa
+
+    ϕ_oceanic_crust = sind(0)
+    C_oceanic_crust = 5e6Pa
+
+    ϕ_oceanic_litho = sind(0)
+    C_oceanic_litho = 5e6Pa
+
+    ϕ_cont_crust    = sind(20)
+    C_cont_crust    = 30e6Pa
+    
+    soft_C  = LinearSoftening((C_oceanic_litho.val*0.95, C_oceanic_litho.val), (0.1, 0.5))
+    
+    # common physical properties
+    α     = 3e-5 / K
+    Cp    = 1000 * J / kg * K
+    η_reg = 1e18Pa * s 
+
+
+    # Define rheolgy struct
+    rheology = (
+        # Name              = "dry olivine - Hirth_Kohlstedt_2003",
+        SetMaterialParams(;
+            Phase             = 1,
+            Density           = PT_Density(; ρ0=3.3e3kg / m^3, α = α, β = 0e0 / Pa, T0 = 273K),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
+            Conductivity      = ConstantConductivity(; k = 3Watt/m/K),
+            CompositeRheology = CompositeRheology( 
+                    (
+                        LinearViscous(;η=1e19Pa*s),
+                        # disl_dry_olivine, 
+                        # diff_dry_olivine,
+                        # ConstantElasticity(; G=5e10Pa, ν=0.5),
+                        # DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+                    ) 
+                ),
+            RadioactiveHeat   = ConstantRadioactiveHeat(6.6667e-12Watt/kg),
+            # Elasticity        = ConstantElasticity(; G=5e10Pa, ν=0.5),
+            Gravity           = ConstantGravity(; g=9.81m/s^2),
+            CharDim           = CharDim
+        ),
+        # Name              = "oceanic crust",
+        SetMaterialParams(;
+            Phase             = 2,
+            Density           = PT_Density(; ρ0=3.3e3kg / m^3, α = α, β = 0e0 / Pa, T0 = 273K),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
+            Conductivity      = ConstantConductivity(; k = 3Watt/m/K),
+            CompositeRheology = CompositeRheology( 
+                (
+                    LinearViscous(;η=1e20Pa*s),
+                    # disl_oceanic_crust,
+                    # ConstantElasticity(; G=5e10Pa, ν=0.5),
+                    # DruckerPrager_regularised(; C = C_oceanic_crust, ϕ = ϕ_oceanic_crust, η_vp=η_reg, Ψ=0.0, softening_C = soft_C) # non-regularized plasticity
+                ) 
+            ),
+            RadioactiveHeat   = ConstantRadioactiveHeat(2.333e-10Watt/kg),
+            # Elasticity        = ConstantElasticity(; G=5e10Pa, ν=0.5),
+            CharDim           = CharDim
+        ),
+        # Name              = "oceanic lithosphere",
+        SetMaterialParams(;
+            Phase             = 3,
+            Density           = PT_Density(; ρ0=3.3e3kg / m^3, α = α, β = 0e0 / Pa, T0 = 273K),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
+            Conductivity      = ConstantConductivity(; k = 3Watt/m/K),
+            CompositeRheology = CompositeRheology( 
+                    (
+                        LinearViscous(;η=1e19Pa*s),
+                        # disl_dry_olivine, 
+                        # diff_dry_olivine,
+                        # ConstantElasticity(; G=5e10Pa, ν=0.5),
+                        # DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+                    ) 
+                ),
+            RadioactiveHeat   = ConstantRadioactiveHeat(6.6667e-12Watt/kg),
+            # Elasticity        = ConstantElasticity(; G=5e10Pa, ν=0.5),
+            CharDim           = CharDim
+        ),
+        # Name              = "continental crust",
+        SetMaterialParams(;
+            Phase             = 4,
+            Density           = PT_Density(; ρ0=2.7e3kg / m^3, α = α, β = 0e0 / Pa, T0 = 273K),
+            RadioactiveHeat   = ConstantRadioactiveHeat(5.3571e-10Watt/kg),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
+            Conductivity      = ConstantConductivity(; k = 3Watt/m/K),
+            CompositeRheology = CompositeRheology( 
+                (
+                    LinearViscous(;η=1e21Pa*s),
+                    # disl_cont_crust, 
+                    # ConstantElasticity(; G=5e10Pa, ν=0.5),
+                    # DruckerPrager_regularised(; C = C_cont_crust, ϕ = ϕ_cont_crust, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+                ) 
+            ),
+            # Elasticity        = ConstantElasticity(; G=5e10Pa, ν=0.5),
+            CharDim           = CharDim
+        ),
+        # Name              = "continental lithosphere",
+        SetMaterialParams(;
+            Phase             = 5,
+            Density           = PT_Density(; ρ0=3.3e3kg / m^3, α = α, β = 0e0 / Pa, T0 = 273K),
+            HeatCapacity      = ConstantHeatCapacity(; Cp=Cp),
+            Conductivity      = ConstantConductivity(; k = 3Watt/m/K),
+            CompositeRheology = CompositeRheology( 
+                    (
+                        LinearViscous(;η=1e19Pa*s),
+                        # disl_dry_olivine, 
+                        # diff_dry_olivine,
+                        # ConstantElasticity(; G=5e10Pa, ν=0.5),
+                        # DruckerPrager_regularised(; C = C_dry_olivine, ϕ=ϕ_dry_olivine, η_vp=η_reg, Ψ=0.0) # non-regularized plasticity
+                    ) 
+                ),
+            RadioactiveHeat   = ConstantRadioactiveHeat(6.6667e-12Watt/kg),
+            # Elasticity        = ConstantElasticity(; G=5e10Pa, ν=0.5),
+            CharDim           = CharDim
+        ),
+    )
+end
+
+function init_phases!(phases, phase_grid, particles, xvi)
+    ni = size(phases)
+    @parallel (@idx ni) _init_phases!(phases, phase_grid, particles.coords, particles.index, xvi)
+end
+
+@parallel_indices (I...) function _init_phases!(phases, phase_grid, pcoords::NTuple{N, T}, index, xvi) where {N,T}
+
+    ni = size(phases)
+
+    for ip in JustRelax.cellaxes(phases)
+        # quick escape
+        @cell(index[ip, I...]) == 0 && continue
+
+        pᵢ = ntuple(Val(N)) do i 
+            @cell pcoords[i][ip, I...]
+        end
+
+        d = Inf # distance to the nearest particle
+        particle_phase = -1
+        for offi in 0:1, offj in 0:1
+            ii = I[1] + offi
+            jj = I[2] + offj
+
+            !(ii ≤ ni[1]) && continue
+            !(jj ≤ ni[2]) && continue
+
+            xvᵢ = (
+                xvi[1][ii], 
+                xvi[2][jj], 
+            )
+            d_ijk = √(sum((pᵢ[i] - xvᵢ[i])^2 for i in 1:N))
+            if d_ijk < d
+                d = d_ijk
+                particle_phase = phase_grid[ii, jj]
+            end
+        end
+        @cell phases[ip, I...] = Float64(particle_phase)
+    end
+
+    return nothing
+end