diff --git a/src/Fields.jl b/src/Fields.jl
index 18aa5f1..de244c1 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -223,13 +223,12 @@ See also: [`set!`](@ref), [`compatible`](@ref).
 function op!(f::F, ϕ::Field, ϕs::Field...; check::Bool=true) where {F}
     check && check_compatible(ϕ, ϕs...)
     g = getgrid(ϕ)
-    idxs = CartesianIndices(g)
+    idxs = interior(g)  # padded indices of the interior cells, for `cell_center`
     ϕint = interior(ϕ)
     ϕsint = map(interior, ϕs)
     AK.foreachindex(ϕint; backend=AK.get_backend(g)) do j
-        idx = idxs[j]
-        x = cell_center(g, idx)
-        ϕint[j] = f(x, ϕint[idx], map(φ -> φ[idx], ϕsint)...)
+        x = cell_center(g, idxs[j])
+        ϕint[j] = f(x, ϕint[j], map(φ -> φ[j], ϕsint)...)
     end
     return ϕ
 end
diff --git a/test/fields.jl b/test/fields.jl
index 9c26374..ad355e7 100644
--- a/test/fields.jl
+++ b/test/fields.jl
@@ -28,6 +28,15 @@
         @test eltype(u32) === Float32
     end
 
+    @testset "op! sees the same cell as set!" begin
+        # Non-square 2D grid: a padded/interior index mix-up shifts or scrambles cells.
+        g = CartesianGrid(((0.0, 1.0), (0.0, 2.0)), (3, 5))
+        f(x) = x[1] + 10 * x[2]
+        u = set!(f, scalar_field(g))
+        op!((x, a) -> a - f(x), u)
+        @test maximum(abs, collect(interior(u))) < 1e-12
+    end
+
     @testset "component extraction" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (3, 3))
         v = set!(x -> SVector(x[1], -x[2]), vector_field(g))