diff --git a/scripts_future_API/tm_stokes_wip.jl b/scripts_future_API/tm_stokes_wip.jl
index 64caa219..f9280834 100644
--- a/scripts_future_API/tm_stokes_wip.jl
+++ b/scripts_future_API/tm_stokes_wip.jl
@@ -31,7 +31,7 @@ boundary_conditions = (
 r       = 0.7
 re_mech = 5π
 lτ_re_m = minimum(extent(grid)) / re_mech
-vdτ     = minimum(spacing(grid)) / sqrt(5.1)
+vdτ     = minimum(spacing(grid)) / sqrt(8.1)
 θ_dτ    = lτ_re_m * (r + 4 / 3) / vdτ
 dτ_r    = 1.0 ./ (θ_dτ .+ 1.0)
 nudτ    = vdτ * lτ_re_m
@@ -48,8 +48,8 @@ other_fields = (
     A = Field(backend, grid, Center()),
 )
 
-init_incl(x, y, z, x0, y0, z0, r, ηi, ηm) = ifelse((x-x0)^2 + (y-y0)^2 + (z-z0)^2 < r^2, ηi, ηm)
-set!(other_fields.A, grid, init_incl; continuous = true, parameters = (x0 = 0.0, y0 = 0.0, z0 = 0.5, r = 0.2, ηi = 1e-1, ηm = 1.0))
+init_incl(x, y, z, x0, y0, z0, r, Ai, Am) = ifelse((x-x0)^2 + (y-y0)^2 + (z-z0)^2 < r^2, Ai, Am)
+set!(other_fields.A, grid, init_incl; continuous = true, parameters = (x0 = 0.0, y0 = 0.0, z0 = 0.5, r = 0.2, Ai = 1e1, Am = 1.0))
 
 model = IsothermalFullStokesModel(;
     backend,
@@ -60,6 +60,12 @@ model = IsothermalFullStokesModel(;
     other_fields
 )
 
+fig = Figure(resolution=(1000,1000), fontsize=32)
+ax  = Axis(fig[1,1][1,1]; aspect=DataAspect(), xlabel="x", ylabel="y")
+
+plt = heatmap!(ax, xcenters(grid), ycenters(grid), interior(model.fields.Pr)[:, :, size(grid,3)÷2]; colormap=:turbo)
+Colorbar(fig[1,1][1,2], plt)
+
 set!(model.fields.Pr, 0.0)
 foreach(x -> set!(x, 0.0), model.fields.τ)
 Isothermal.apply_bcs!(model.backend, model.grid, model.fields, model.boundary_conditions.stress)
@@ -74,17 +80,12 @@ Isothermal.apply_bcs!(model.backend, model.grid, model.fields, model.boundary_co
 set!(model.fields.η, other_fields.A)
 extrapolate!(data(model.fields.η))
 
-
 for it in 1:100
     advance_iteration!(model, 0.0, 1.0; async = false)
+    if it % 10 == 0
+        plt[3][] = interior(model.fields.Pr)[:, :, size(grid,3)÷2]
+        yield()
+    end
 end
 
-fig = Figure(resolution=(1000,1000), fontsize=32)
-ax  = Axis(fig[1,1][1,1]; aspect=DataAspect(), xlabel="x", ylabel="y")
-
-plt = heatmap!(ax, xcenters(grid), ycenters(grid), interior(model.fields.Pr)[:, :, size(grid,3)÷2]; colormap=:turbo)
-Colorbar(fig[1,1][1,2], plt)
-
-plt[3][] = interior(model.fields.Pr)[:, :, size(grid,3)÷2]
 
-yield()
\ No newline at end of file
diff --git a/src/models/full_stokes/isothermal/isothermal.jl b/src/models/full_stokes/isothermal/isothermal.jl
index da670cbc..82a7f3cf 100644
--- a/src/models/full_stokes/isothermal/isothermal.jl
+++ b/src/models/full_stokes/isothermal/isothermal.jl
@@ -241,7 +241,7 @@ function advance_iteration!(model::IsothermalFullStokesModel, t, Δt; async = tr
     set_bcs!(bcs) = apply_bcs!(model.backend, model.grid, model.fields, bcs)
 
     # stress
-    update_σ!(backend, 256, (nx, ny, nz))(interior(Pr), interior(τ), V, η, Δτ, Δ)
+    update_σ!(backend, 256, (nx+1, ny+1, nz+1))(interior(Pr), interior(τ), V, η, Δτ, Δ)
     set_bcs!(model.boundary_conditions.stress)
     # velocity
     update_V!(backend, 256, (nx + 1, ny + 1, nz + 1))(interior(V), Pr, τ, η, Δτ, Δ)
diff --git a/src/models/full_stokes/isothermal/kernels.jl b/src/models/full_stokes/isothermal/kernels.jl
index 7eecce09..0c422aa1 100644
--- a/src/models/full_stokes/isothermal/kernels.jl
+++ b/src/models/full_stokes/isothermal/kernels.jl
@@ -11,9 +11,9 @@ end
 @kernel function update_σ!(Pr, τ, V, η, Δτ, Δ)
     ix, iy, iz = @index(Global, NTuple)
     @inbounds if checkbounds(Bool, Pr, ix, iy, iz)
-        ε̇xx = @∂_x(V.x) / Δ.x
-        ε̇yy = @∂_y(V.y) / Δ.y
-        ε̇zz = @∂_z(V.z) / Δ.z
+        ε̇xx = @∂_xi(V.x) / Δ.x
+        ε̇yy = @∂_yi(V.y) / Δ.y
+        ε̇zz = @∂_zi(V.z) / Δ.z
         ∇V = ε̇xx + ε̇yy + ε̇zz
         # hydrostatic
         @all(Pr) -= ∇V * @inn(η) * Δτ.Pr
@@ -23,13 +23,19 @@ end
         @all(τ.zz) -= (@all(τ.zz) - 2.0 * @inn(η) * (ε̇zz - ∇V / 3.0)) * Δτ.τ.zz
     end
     @inbounds if checkbounds(Bool, τ.xy, ix, iy, iz)
-        @all(τ.xy) -= (@all(τ.xy) - @av_xy(η) * (@∂_x(V.y) / Δ.x + @∂_y(V.x) / Δ.y)) * Δτ.τ.xy
+        exy = (V.y[ix+1, iy, iz+1] - V.y[ix, iy, iz+1]) / Δ.x +
+              (V.x[ix, iy+1, iz+1] - V.x[ix, iy, iz+1]) / Δ.y
+        @all(τ.xy) -= (@all(τ.xy) - @av_xyi(η) * exy) * Δτ.τ.xy
     end
     @inbounds if checkbounds(Bool, τ.xz, ix, iy, iz)
-        @all(τ.xz) -= (@all(τ.xz) - @av_xz(η) * (@∂_x(V.z) / Δ.x + @∂_z(V.x) / Δ.z)) * Δτ.τ.xz
+        exz = (V.z[ix+1, iy+1, iz  ] - V.z[ix, iy+1, iz]) / Δ.x +
+              (V.x[ix  , iy+1, iz+1] - V.x[ix, iy+1, iz]) / Δ.z
+        @all(τ.xz) -= (@all(τ.xz) - @av_xzi(η) * exz) * Δτ.τ.xz
     end
     @inbounds if checkbounds(Bool, τ.yz, ix, iy, iz)
-        @all(τ.yz) -= (@all(τ.yz) - @av_yz(η) * (@∂_y(V.z) / Δ.y + @∂_z(V.y) / Δ.z)) * Δτ.τ.yz
+        eyz = (V.z[ix+1, iy+1, iz  ] - V.z[ix+1, iy, iz]) / Δ.y +
+              (V.y[ix+1, iy  , iz+1] - V.y[ix+1, iy, iz]) / Δ.z
+        @all(τ.yz) -= (@all(τ.yz) - @av_yzi(η) * eyz) * Δτ.τ.yz
     end
 end