diff --git a/DESIGN.md b/DESIGN.md
index 3ee2f06..6a5fbd8 100644
--- a/DESIGN.md
+++ b/DESIGN.md
@@ -603,7 +603,7 @@ and drifts by an ulp — invisible to a stencil, which reads only spacing, but e
 to make a coordinate-assembled RHS and therefore the Krylov iteration count depend
 on `nparts`. It is bit-for-bit a no-op for undistributed grids and forest leaf
 grids, where `first(local_range[d]) == 1`. The user-facing surface is
-`boundary_rhs(P)`, `set!(::MultiDeviceVector, P, fun)`, `assemble_rhs(P, f)`, and
+`boundary_rhs(P)`, `set!(fun, ::MultiDeviceVector, P)`, `assemble_rhs(P, f)`, and
 `local_grids(P)`. Only `prepare_distributed` carries a `distributed` qualifier,
 because only it shadows a single-device function; everything downstream dispatches
 on `P` and is named for what it computes, not for where it runs.
diff --git a/Project.toml b/Project.toml
index 6dcc527..5829ea6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "MatrixFreeOperators"
 uuid = "bf251007-91fc-422b-8727-f8d5215c1c76"
-version = "0.1.0"
+version = "0.2.0"
 authors = ["Kyle Beggs <beggskw@gmail.com>"]
 
 [deps]
diff --git a/README.md b/README.md
index 19e3dae..74425a5 100644
--- a/README.md
+++ b/README.md
@@ -101,13 +101,13 @@ g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (32, 32))
 
 # σ(x,y) = 1 + xy as a coefficient field. Fields are device arrays plus
 # grid metadata; set! fills them from a function of position.
-σ = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+σ = set!(x -> 1 + x[1] * x[2], scalar_field(g))
 
 # The operator is composed symbolically — no matrix is ever assembled.
 K = scaling(σ) - laplacian(g)
 
 # Manufactured right-hand side for u = sin(πx)sin(πy).
-f = set!(scalar_field(g), x -> (2π^2 + 1 + x[1] * x[2]) * sinpi(x[1]) * sinpi(x[2]))
+f = set!(x -> (2π^2 + 1 + x[1] * x[2]) * sinpi(x[1]) * sinpi(x[2]), scalar_field(g))
 
 # prepare walks the operator tree once and allocates all scratch buffers;
 # the result supports mul!/size/eltype with zero steady-state allocations,
@@ -116,7 +116,7 @@ P = prepare(K, scalar_field(g))
 u, stats = cg(P, flatten(f))
 
 # Compare against the exact solution on the interior DOFs.
-u_exact = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
+u_exact = flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
 maximum(abs, u .- u_exact)   # ~1e-3, second-order accurate
 ```
 
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 512b809..f087635 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -6,5 +6,5 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 BenchmarkTools = "1.5"
-MatrixFreeOperators = "0.1"
+MatrixFreeOperators = "0.2"
 julia = "1.10"
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index efbb278..ec95946 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -16,15 +16,15 @@ g3 = CartesianGrid(((0.0, 1.0), (0.0, 1.0), (0.0, 1.0)), (64, 64, 64))
 for (dim, g) in (("2D 256²", g2), ("3D 64³", g3))
     L = laplacian(g)
     P = prepare(L, scalar_field(g))
-    x = flatten(set!(scalar_field(g), p -> sin(4p[1]) + cos(3p[end])))
+    x = flatten(set!(p -> sin(4p[1]) + cos(3p[end]), scalar_field(g)))
     y = similar(x)
     SUITE["grid"][dim]["laplacian prepare"] = @benchmarkable prepare($L, $(scalar_field(g)))
     SUITE["grid"][dim]["laplacian mul!"] = @benchmarkable mul!($y, $P, $x)
 end
 
-xs2 = flatten(set!(scalar_field(g2), p -> sin(4p[1]) + cos(3p[2])))
+xs2 = flatten(set!(p -> sin(4p[1]) + cos(3p[2]), scalar_field(g2)))
 ys2 = similar(xs2)
-xv2 = flatten(set!(vector_field(g2), p -> SVector(sin(p[2]), cos(p[1]))))
+xv2 = flatten(set!(p -> SVector(sin(p[2]), cos(p[1])), vector_field(g2)))
 yv2 = similar(xv2)
 
 PG = prepare(gradient(g2), scalar_field(g2))
@@ -33,11 +33,11 @@ SUITE["grid"]["2D 256²"]["gradient mul!"] = @benchmarkable mul!($yv2, $PG, $xs2
 PD = prepare(divergence(g2), vector_field(g2))
 SUITE["grid"]["2D 256²"]["divergence mul!"] = @benchmarkable mul!($ys2, $PD, $xv2)
 
-vel = set!(vector_field(g2), p -> SVector(sin(p[2]), cos(p[1])))
+vel = set!(p -> SVector(sin(p[2]), cos(p[1])), vector_field(g2))
 PA = prepare(advection(g2, vel), scalar_field(g2))
 SUITE["grid"]["2D 256²"]["advection mul!"] = @benchmarkable mul!($ys2, $PA, $xs2)
 
-κ = set!(scalar_field(g2), p -> 1 + 0.5 * sin(p[1]))
+κ = set!(p -> 1 + 0.5 * sin(p[1]), scalar_field(g2))
 PS = prepare(2.0 * laplacian(g2) + scaling(κ), scalar_field(g2))
 SUITE["grid"]["2D 256²"]["2λ + κ·I mul!"] = @benchmarkable mul!($ys2, $PS, $xs2)
 
@@ -50,15 +50,15 @@ SUITE["grid"]["2D 256²"]["∇·(κ∇u) mul!"] = @benchmarkable mul!($ys2, $PK,
 # benchpkg runs this file on the PR's base branch too, where Diffusion does not exist.
 if isdefined(MatrixFreeOperators, :Diffusion)
     for (dim, g, xf, yf) in (("2D 256²", g2, xs2, ys2),)
-        D = diffusion(g, set!(scalar_field(g), p -> 1 + 0.5 * sin(p[1])))
+        D = diffusion(g, set!(p -> 1 + 0.5 * sin(p[1]), scalar_field(g)))
         PDf = prepare(D, scalar_field(g))
         SUITE["grid"][dim]["diffusion prepare"] = @benchmarkable prepare(
             $D, $(scalar_field(g))
         )
         SUITE["grid"][dim]["diffusion mul!"] = @benchmarkable mul!($yf, $PDf, $xf)
     end
-    D3 = diffusion(g3, set!(scalar_field(g3), p -> 1 + 0.5 * sin(p[1])))
-    x3 = flatten(set!(scalar_field(g3), p -> sin(4p[1]) + cos(3p[3])))
+    D3 = diffusion(g3, set!(p -> 1 + 0.5 * sin(p[1]), scalar_field(g3)))
+    x3 = flatten(set!(p -> sin(4p[1]) + cos(3p[3]), scalar_field(g3)))
     y3 = similar(x3)
     PD3 = prepare(D3, scalar_field(g3))
     SUITE["grid"]["3D 64³"]["diffusion mul!"] = @benchmarkable mul!($y3, $PD3, $x3)
@@ -73,10 +73,10 @@ end
 # revisions without the distributed seam have no `_slab_op`.
 if isdefined(MatrixFreeOperators, :Diffusion) && isdefined(MatrixFreeOperators, :_slab_op)
     for (dim, g) in (("2D 256²", g2), ("3D 64³", g3))
-        D = diffusion(g, set!(scalar_field(g), p -> 1 + 0.5 * sin(p[1])))
+        D = diffusion(g, set!(p -> 1 + 0.5 * sin(p[1]), scalar_field(g)))
         lg = partition_grid(g, 2)[1]
         Ds = MatrixFreeOperators._slab_op(D, lg)
-        ȳs = set!(scalar_field(lg), p -> sin(4p[1]) + cos(3p[end]))
+        ȳs = set!(p -> sin(4p[1]) + cos(3p[end]), scalar_field(lg))
         x̄s = scalar_field(lg)
         SUITE["grid"][dim]["diffusion slab apply!"] =
             @benchmarkable apply!($x̄s, $Ds, $ȳs, $lg, 1.0, 0.0)
@@ -102,8 +102,8 @@ end
 # rank-changers gather unconditionally; a Laplacian would shortcut to its forward
 # action on this all-physical grid and measure nothing.
 
-sadj = set!(scalar_field(g2), p -> sin(4p[1]) + cos(3p[2]))
-vadj = set!(vector_field(g2), p -> SVector(sin(p[2]), cos(p[1])))
+sadj = set!(p -> sin(4p[1]) + cos(3p[2]), scalar_field(g2))
+vadj = set!(p -> SVector(sin(p[2]), cos(p[1])), vector_field(g2))
 sadj_out = scalar_field(g2)
 vadj_out = vector_field(g2)
 Dx, Dy = derivative(g2, 1), derivative(g2, 2)
@@ -132,7 +132,7 @@ SUITE["grid"]["2D 256²"]["adjoint(∂x + ∂y) mul!"] = @benchmarkable mul!($ys
 # 8×8 root tiling of 32² blocks (64 uniform leaves) — same DOFs as the 2D grid above,
 # so the forest overhead (halo exchange + per-leaf dispatch) is directly comparable.
 bf = BlockForest(g2; blocksize=(32, 32), maxlevel=2)
-xb = set!(scalar_field(bf), p -> sin(4p[1]) + cos(3p[2]))
+xb = set!(p -> sin(4p[1]) + cos(3p[2]), scalar_field(bf))
 Lf = laplacian(bf)
 Pf = prepare(Lf, scalar_field(bf))
 xf = flatten(xb)
@@ -159,7 +159,7 @@ end
 # Same 32² blocks, so the extra cost over the uniform leg is coarse–fine work.
 bfr = BlockForest(g2; blocksize=(32, 32), maxlevel=2)
 refine!(bfr, p -> p[1] < 0.5)
-xbr = set!(scalar_field(bfr), p -> sin(4p[1]) + cos(3p[2]))
+xbr = set!(p -> sin(4p[1]) + cos(3p[2]), scalar_field(bfr))
 Lr = laplacian(bfr)
 
 SUITE["forest"]["2D refined"]["halo_update!"] = @benchmarkable halo_update!($xbr, $bfr)
@@ -184,7 +184,7 @@ SUITE["forest"]["2D refined"]["laplacian apply_adjoint!"] =
 g3r = CartesianGrid(((0.0, 1.0), (0.0, 1.0), (0.0, 1.0)), (32, 32, 32))
 bfr3 = BlockForest(g3r; blocksize=(8, 8, 8), maxlevel=2)
 refine!(bfr3, p -> p[1] < 0.5)
-xbr3 = set!(scalar_field(bfr3), p -> sin(4p[1]) + cos(3p[3]))
+xbr3 = set!(p -> sin(4p[1]) + cos(3p[3]), scalar_field(bfr3))
 
 SUITE["forest"]["3D refined"]["halo_update!"] = @benchmarkable halo_update!($xbr3, $bfr3)
 SUITE["forest"]["3D refined"]["halo_update_adjoint!"] =
diff --git a/benchmark/gpu.jl b/benchmark/gpu.jl
index d3bc518..61e4f1f 100644
--- a/benchmark/gpu.jl
+++ b/benchmark/gpu.jl
@@ -33,7 +33,7 @@ println("laplacian mul! — CPU per-leaf / GPU per-leaf / GPU packed (min time)"
     "forest", "nleaves", "DOFs", "CPU/leaf", "GPU/leaf", "GPU packed", "packed speedup", "host alloc")
 for (n, refined) in ((256, false), (512, false), (1024, false), (2048, false), (512, true))
     bf = make_forest(n; refined)
-    u = set!(scalar_field(bf), x -> sin(4x[1]) + cos(3x[2]))
+    u = set!(x -> sin(4x[1]) + cos(3x[2]), scalar_field(bf))
     L = laplacian(bf)
     x = flatten(u)
     y = similar(x)
@@ -70,7 +70,7 @@ println("\nhalo_update! on the device packed field — batched kernels vs per-de
     "forest", "nleaves", "kernels", "loop", "kernel speedup")
 for (n, refined) in ((512, false), (2048, false), (512, true))
     bf = make_forest(n; refined)
-    u = set!(scalar_field(bf), x -> sin(4x[1]) + cos(3x[2]))
+    u = set!(x -> sin(4x[1]) + cos(3x[2]), scalar_field(bf))
     pg = Adapt.adapt(CuArray, pack(u))
     sched = MFO._exchange_schedule(bf)
     halo_update!(pg, pg.grid)                        # warm the _device_schedule cache
diff --git a/docs/index.qmd b/docs/index.qmd
index f997e07..63de2eb 100644
--- a/docs/index.qmd
+++ b/docs/index.qmd
@@ -36,12 +36,12 @@ using MatrixFreeOperators, Krylov, LinearAlgebra
 
 # -Δu + σu = f on (0,1)², homogeneous Dirichlet, manufactured solution
 g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (32, 32))
-σ = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+σ = set!(x -> 1 + x[1] * x[2], scalar_field(g))
 K = scaling(σ) - laplacian(g)
 
-f = set!(scalar_field(g), x -> (2π^2 + 1 + x[1] * x[2]) * sinpi(x[1]) * sinpi(x[2]))
+f = set!(x -> (2π^2 + 1 + x[1] * x[2]) * sinpi(x[1]) * sinpi(x[2]), scalar_field(g))
 P = prepare(K, scalar_field(g))
 u, stats = cg(P, flatten(f))
 
-maximum(abs, u .- flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2]))))
+maximum(abs, u .- flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g))))
 ```
diff --git a/docs/pages/amr.qmd b/docs/pages/amr.qmd
index ed79125..7722f2e 100644
--- a/docs/pages/amr.qmd
+++ b/docs/pages/amr.qmd
@@ -171,8 +171,8 @@ different, block-by-block order), so any order-independent reduction matches:
 g  = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
 bf = BlockForest(g; blocksize = (4, 4), maxlevel = 3)   # level 0 ⇒ same 8×8
 
-u  = set!(scalar_field(g),  x -> sinpi(x[1]) * sinpi(x[2]))
-uf = set!(scalar_field(bf), x -> sinpi(x[1]) * sinpi(x[2]))
+u  = set!(x -> sinpi(x[1]) * sinpi(x[2]),  scalar_field(g))
+uf = set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(bf))
 
 norm(flatten(laplacian(g) * u)) ≈ norm(flatten(laplacian(bf) * uf))
 ```
@@ -188,7 +188,7 @@ conflicts), re-establishes 2:1 balance, and returns freshly allocated fields wit
 the data carried across:
 
 ```{julia}
-bump = set!(scalar_field(forest), x -> exp(-((x[1] - 0.5)^2 + (x[2] - 0.5)^2) / 0.01))
+bump = set!(x -> exp(-((x[1] - 0.5)^2 + (x[2] - 0.5)^2) / 0.01), scalar_field(forest))
 bump = regrid!(bump; refine = b -> maximum(abs, interior(b)) > 0.5)
 length(collect(leaves(forest)))                  # resolution now follows the bump
 ```
diff --git a/docs/pages/autodiff.qmd b/docs/pages/autodiff.qmd
index bdfae8a..7c7ef31 100644
--- a/docs/pages/autodiff.qmd
+++ b/docs/pages/autodiff.qmd
@@ -45,7 +45,7 @@ function diffusion_loss(κdata, udata, w, gg)
 end
 
 rng = MersenneTwister(1)
-u = set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2]))
+u = set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g))
 w = rand(rng, local_size(g)...)
 κ = 1.0 .+ rand(rng, padded_size(g)...)
 
@@ -156,12 +156,12 @@ using Enzyme, StaticArrays, LinearAlgebra
 
 g1 = CartesianGrid(((0.0, 2π),), (32,); bc=((Periodic(), Periodic()),))
 F  = advection(g1, SelfAdvection())          # nonlinear u·∇u
-u0 = set!(vector_field(g1), x -> SVector(2 + sin(x[1])))
+u0 = set!(x -> SVector(2 + sin(x[1])), vector_field(g1))
 
 J  = linearize(F, u0)                        # FiniteDifferenceJVP — the default
 Ja = linearize(F, u0, EnzymeJVP())           # exact, and has a transpose
 
-v = set!(vector_field(g1), x -> SVector(cos(2x[1])))
+v = set!(x -> SVector(cos(2x[1])), vector_field(g1))
 Jv  = collect(interior(apply(J,  copy(v))))
 Jav = collect(interior(apply(Ja, copy(v))))
 
@@ -174,7 +174,7 @@ Only the Enzyme-backed one has a transpose:
 # A ramp rather than a trig mode: on this periodic grid J·v is orthogonal to the
 # low harmonics, so a sin/cos test vector would make both sides vanish and the
 # identity would hold vacuously.
-w = set!(vector_field(g1), x -> SVector(x[1] / 2π))
+w = set!(x -> SVector(x[1] / 2π), vector_field(g1))
 lhs = sum(dot.(collect(interior(apply(Ja, copy(v)))), collect(interior(w))))
 rhs = sum(dot.(collect(interior(v)), collect(interior(apply(adjoint(Ja), copy(w))))))
 (lhs, rhs, abs(lhs - rhs))                   # ⟨Jv, w⟩ = ⟨v, Jᵀw⟩
diff --git a/docs/pages/distributed.qmd b/docs/pages/distributed.qmd
index 5515efc..d9b7bd7 100644
--- a/docs/pages/distributed.qmd
+++ b/docs/pages/distributed.qmd
@@ -229,7 +229,7 @@ assembles the lift of its own slab and the result comes back as a
 b = assemble_rhs(P, x -> sinpi(x[1]) * sinpi(x[2]))   # f - boundary_rhs(P)
 
 # ...or in two steps, if you want the pieces
-src = set!(MultiDeviceVector{Float64}(undef, P.spec), P, x -> sinpi(x[1]))
+src = set!(x -> sinpi(x[1]), MultiDeviceVector{Float64}(undef, P.spec), P)
 b   = src .- boundary_rhs(P)
 ```
 
@@ -247,7 +247,7 @@ Bring your own source data with [`local_grids`](../pages/api.qmd), which hands
 back each partition's slab grid, on the host:
 
 ```julia
-fields = [set!(scalar_field(lg), myfun) for lg in local_grids(P)]
+fields = [set!(myfun, scalar_field(lg)) for lg in local_grids(P)]
 b = assemble_rhs(P, fields)
 ```
 
diff --git a/docs/pages/internals.qmd b/docs/pages/internals.qmd
index 89cdc78..f5beb37 100644
--- a/docs/pages/internals.qmd
+++ b/docs/pages/internals.qmd
@@ -1084,7 +1084,7 @@ g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (256, 256))  # homogeneous Dirichlet
 L = -laplacian(g)
 P = prepare(L)
 M = MultigridPreconditioner(L)
-b = flatten(set!(scalar_field(g), x -> 2π^2 * sinpi(x[1]) * sinpi(x[2])))
+b = flatten(set!(x -> 2π^2 * sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
 u, stats = Krylov.cg(P, b; M)
 ```
 
diff --git a/examples/adaptive_poisson.jl b/examples/adaptive_poisson.jl
index 335441b..848309a 100644
--- a/examples/adaptive_poisson.jl
+++ b/examples/adaptive_poisson.jl
@@ -25,7 +25,7 @@ function l1_error(u, bf)
     e = 0.0
     for i in 1:MFO.nleaves(bf)
         sp = MFO._leaf_spacing(bf, bf.forest.leaves[i].level)
-        ref = set!(similar(MFO.block(u, i)), uexact)
+        ref = set!(uexact, similar(MFO.block(u, i)))
         e += sum(abs, interior(MFO.block(u, i)) .- interior(ref)) * prod(sp)
     end
     return e
@@ -33,7 +33,7 @@ end
 
 for cycle in 1:4
     P = prepare(laplacian(bf), u)
-    rhs = .-flatten(set!(scalar_field(bf), f))          # Δu = -f; Dirichlet lift is zero
+    rhs = .-flatten(set!(f, scalar_field(bf)))          # Δu = -f; Dirichlet lift is zero
     sol, stats = Krylov.gmres(P, rhs; rtol=1e-10)       # nonsymmetric on an adapted forest
     flat_to_interior!(u, sol)
     lo, hi = extrema(k -> k.level, bf.forest.leaves)
diff --git a/examples/heat_equation.jl b/examples/heat_equation.jl
index 962e086..6be38b2 100644
--- a/examples/heat_equation.jl
+++ b/examples/heat_equation.jl
@@ -20,8 +20,8 @@ nsteps = nframes * steps_per_frame
 dt = (tspan[2] - tspan[1]) / nsteps
 # 20 Gaussian blobs at random centers; diameter ≈ 4σ drawn uniform in [0.1, 0.3]
 blobs = [(rand(), rand(), 0.2 * (0.5 + rand())) for _ in 1:20]
-u0 = set!(scalar_field(g), x ->
-    sum(exp(-((x[1] - cx)^2 + (x[2] - cy)^2) / (s^2 / 8)) for (cx, cy, s) in blobs))
+u0 = set!(x ->
+    sum(exp(-((x[1] - cx)^2 + (x[2] - cy)^2) / (s^2 / 8)) for (cx, cy, s) in blobs), scalar_field(g))
 interior(u0) ./= maximum(interior(u0))
 
 function step!(u, du, L, αdt)
diff --git a/examples/inverse_diffusion.jl b/examples/inverse_diffusion.jl
index 0ad6e8d..1a6e990 100644
--- a/examples/inverse_diffusion.jl
+++ b/examples/inverse_diffusion.jl
@@ -62,10 +62,10 @@ end
 
 # This example uses one drive. The data are sensitive to κ through the flux κ∇u,
 # so recovery remains weak where this excitation's gradient is small.
-u = set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2]))
+u = set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g))
 
 rng = MersenneTwister(20260731)
-κ★ = set!(scalar_field(g), κ_true)
+κ★ = set!(κ_true, scalar_field(g))
 clean = collect(response(κ★.data, u.data, g))
 obs = clean .+ 0.01 * maximum(abs, clean) .* randn(rng, size(clean))
 
diff --git a/examples/monodomain_amr.jl b/examples/monodomain_amr.jl
index 1256c85..aff5d01 100644
--- a/examples/monodomain_amr.jl
+++ b/examples/monodomain_amr.jl
@@ -100,7 +100,7 @@ function simulate(; snap_times=Float64[], on_frame=nothing, frame_every=250)
                          bc=ntuple(_ -> (Neumann(), Neumann()), 2))
     bf = BlockForest(base; blocksize=(BS, BS), maxlevel=MAXLEV)
 
-    V = set!(scalar_field(bf), x -> x[1] < 4.0 ? 1.0 : 0.0)   # S1 along the left edge
+    V = set!(x -> x[1] < 4.0 ? 1.0 : 0.0, scalar_field(bf))   # S1 along the left edge
     W = scalar_field(bf)
     dv, dw = zeros(BS, BS), zeros(BS, BS)
 
diff --git a/examples/multigrid_poisson.jl b/examples/multigrid_poisson.jl
index 547534d..d9e8acb 100644
--- a/examples/multigrid_poisson.jl
+++ b/examples/multigrid_poisson.jl
@@ -16,8 +16,8 @@ A = prepare(L)
 uexact(x) = 1 + exp(x[1]) * sinpi(x[1]) * sinpi(x[2])
 f(x) = -exp(x[1]) * sinpi(x[2]) * ((1 - 2 * pi^2) * sinpi(x[1]) + 2 * pi * cospi(x[1]))
 # inhomogeneous boundary data enters the RHS through the affine lift
-b = flatten(set!(scalar_field(g), f)) .- flatten(boundary_rhs(L, g))
-uex = flatten(set!(scalar_field(g), uexact))
+b = flatten(set!(f, scalar_field(g))) .- flatten(boundary_rhs(L, g))
+uex = flatten(set!(uexact, scalar_field(g)))
 
 mg = MultigridPreconditioner(L)                          # Jacobi(2/3), levels=:auto
 println(mg)
diff --git a/examples/niederer_benchmark.jl b/examples/niederer_benchmark.jl
index bce33dc..c546bf0 100644
--- a/examples/niederer_benchmark.jl
+++ b/examples/niederer_benchmark.jl
@@ -403,8 +403,8 @@ function simulate(; adaptive::Bool)
         balance!(bf)
     end
 
-    V = set!(scalar_field(bf), _ -> initial_voltage())
-    S = set!(state_field(bf), _ -> initial_state())
+    V = set!(_ -> initial_voltage(), scalar_field(bf))
+    S = set!(_ -> initial_state(), state_field(bf))
     LV = scalar_field(bf)
     Lop = diffusion_operator(bf)
 
@@ -419,8 +419,8 @@ function simulate(; adaptive::Bool)
     diffuse!(V, LV, Lop, bf)
     react!(V, S, bf, 0.0)
     probe_all(V, bf, probes)
-    V = set!(scalar_field(bf), _ -> initial_voltage())
-    S = set!(state_field(bf), _ -> initial_state())
+    V = set!(_ -> initial_voltage(), scalar_field(bf))
+    S = set!(_ -> initial_state(), state_field(bf))
     LV = scalar_field(bf)
 
     snaps, hist = Any[], Tuple{Float64,Int}[]
diff --git a/ext/MatrixFreeOperatorsMDLAExt.jl b/ext/MatrixFreeOperatorsMDLAExt.jl
index 42e9bbc..6c2c6ed 100644
--- a/ext/MatrixFreeOperatorsMDLAExt.jl
+++ b/ext/MatrixFreeOperatorsMDLAExt.jl
@@ -293,7 +293,7 @@ call it concurrently with a solve on the same prepared operator.
 
 ```julia
 P = prepare_distributed(laplacian(g), 2)
-b = set!(MultiDeviceVector{Float64}(undef, P.spec), P, x -> sin(x[1]))
+b = set!(x -> sin(x[1]), MultiDeviceVector{Float64}(undef, P.spec), P)
 b .-= boundary_rhs(P)
 u, stats = Krylov.cg(P, b)
 ```
@@ -333,12 +333,12 @@ MatrixFreeOperators.local_grids(P::MDLAPreparedOperator) =
     [Adapt.adapt(Array, p.grid) for p in P.parts]
 
 """
-    set!(x::MultiDeviceVector, P::MDLAPreparedOperator, fun) -> x
+    set!(fun, x::MultiDeviceVector, P::MDLAPreparedOperator) -> x
 
 Fill `x` with `fun(coords)` evaluated slab-locally on each partition.
 
-The distributed twin of `set!(::Field, fun)`, and exact: `cell_center` evaluates
-at the global cell index, so this is bit-for-bit `MultiDeviceVector(flatten(set!(scalar_field(g), fun)), P.spec)`
+The distributed twin of `set!(fun, ::Field)`, and exact: `cell_center` evaluates
+at the global cell index, so this is bit-for-bit `MultiDeviceVector(flatten(set!(fun, scalar_field(g))), P.spec)`
 without ever building the global field. Uses `P`'s input scratch, so the same
 concurrency caveat as `mul!` applies.
 
@@ -349,7 +349,7 @@ captured host arrays. For anything heavier, build the fields yourself on
 them.
 """
 function MatrixFreeOperators.set!(
-    x::MultiDeviceVector{T}, P::MDLAPreparedOperator{T}, fun
+    fun, x::MultiDeviceVector{T}, P::MDLAPreparedOperator{T}
 ) where {T}
     _dist_set!(P.xpads, fun, P.ctx)
     _dist_map!(P.ctx) do d
@@ -388,7 +388,7 @@ function MatrixFreeOperators.assemble_rhs(P::MDLAPreparedOperator{T}, f) where {
     return x
 end
 
-_source!(x, P::MDLAPreparedOperator, fun) = MatrixFreeOperators.set!(x, P, fun)
+_source!(x, P::MDLAPreparedOperator, fun) = MatrixFreeOperators.set!(fun, x, P)
 function _source!(x, P::MDLAPreparedOperator, fields::AbstractVector)
     length(fields) == length(P.parts) || throw(
         ArgumentError(
diff --git a/src/BlockForest.jl b/src/BlockForest.jl
index e76c891..2d533ff 100644
--- a/src/BlockForest.jl
+++ b/src/BlockForest.jl
@@ -136,6 +136,41 @@ Leaf grids are all-[`Interface`](@ref); physical BCs live on `bf.bc`.
 """
 leaves(bf::BlockForest) = ((key, leaf_grid(bf, key)) for key in bf.forest.leaves)
 
+@inline function _grid_mismatch(a::BlockForest{N}, b::BlockForest{N}) where {N}
+    a === b && return nothing
+    a.blocksize == b.blocksize || return :blocksize
+    a.halo == b.halo || return :halo
+    a.spacing0 == b.spacing0 || return :spacing0
+    a.extent == b.extent || return :extent
+    a.bc === b.bc || return :bc
+    return _forest_mismatch(a.forest, b.forest)
+end
+
+# Forest topology. `===` covers every Adapt twin (they share one `Forest`), so the
+# O(nleaves) leaf compare is reached only for two independently built forests.
+# `maxlevel` is compared even though equal leaf sets would sort identically under a
+# different `maxlevel`: a field's future regrids depend on it, and "same forest"
+# should mean the same forest going forward, not just today.
+function _forest_mismatch(a::Forest{N}, b::Forest{N}) where {N}
+    a === b && return nothing
+    a.nroot == b.nroot || return :nroot
+    a.periodic == b.periodic || return :periodic
+    a.maxlevel == b.maxlevel || return :maxlevel
+    a.leaves == b.leaves || return :leaves
+    return nothing
+end
+
+# Same block shape and the same number of blocks: block `i` of one field can be
+# broadcast against block `i` of the other. Whether block `i` covers the same
+# region in both forests is a *grid* question, not a layout one.
+@inline function _layout_mismatch(a::BlockForest{N}, b::BlockForest{N}) where {N}
+    a === b && return nothing
+    a.blocksize == b.blocksize || return :blocksize
+    a.halo == b.halo || return :halo
+    nleaves(a) == nleaves(b) || return :nleaves
+    return nothing
+end
+
 #--------------------------------------------------------------------------------# Adaptivity
 
 """
diff --git a/src/Fields.jl b/src/Fields.jl
index c79a9c6..97380aa 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -108,6 +108,8 @@ struct Field{L,A<:AbstractArray,G<:AbstractGrid} <: AbstractField
 end
 Field(data::AbstractArray, grid::AbstractGrid) = Field{Center}(data, grid)
 
+AbstractGrid(ϕ::AbstractField) = ϕ.grid
+
 """
     scalar_field(g::AbstractGrid, T=eltype(spacing(g))) -> Field
 
@@ -160,21 +162,72 @@ View of the interior (owned, non-halo) cells of the field.
 interior(f::Field) = view(f.data, interior(f.grid))
 
 """
-    set!(f::Field, fun) -> f
+    set!(fun, ϕ::AbstractField) -> ϕ
+
+Set the interior of `ϕ` to `fun(x)` evaluated at cell centers, where `x` is the
+`SVector` of physical coordinates. Ghost cells are left untouched. On a forest
+field the sweep runs block by block, so `ϕ` must be current with its forest
+(see [`regrid!`](@ref)).
 
-Set the interior of `f` to `fun(x)` evaluated at cell centers, where `x` is the
-`SVector` of physical coordinates.
+`fun` runs on the field's device, so it must be device-compatible: plain
+arithmetic on the coordinate `SVector`, with no captured host arrays. The
+function comes first so the `do`-block form reads naturally.
 
 ### Examples
 
 ```julia
 g = CartesianGrid(((0.0, 2π),), (64,))
-u = set!(scalar_field(g), x -> sin(x[1]))
+u = set!(x -> sin(x[1]), scalar_field(g))
+v = set!(scalar_field(g)) do x
+    exp(-x[1]^2)
+end
 ```
+
+See also: [`op!`](@ref), [`cell_center`](@ref).
 """
-function set!(f::Field, fun::F) where {F}
-    interior(f) .= fun.(cell_center.(Ref(f.grid), interior(f.grid)))
-    return f
+function set!(f::F, ϕ::Field) where {F}
+    g = AbstractGrid(ϕ)
+    map!(interior(ϕ), interior(g)) do idx
+        x = cell_center(g, idx)
+        f(x)
+    end
+    return ϕ
+end
+
+"""
+    op!(fun, ϕ::AbstractField, ϕs::AbstractField...; check=true) -> ϕ
+
+Pointwise update in place: set every interior cell of `ϕ` to
+`fun(x, ϕ[I], ϕs[1][I], ϕs[2][I], …)`, where `x` is the `SVector` of the cell's
+physical coordinates and the remaining arguments are the current values of `ϕ`
+and of each field in `ϕs` at that cell. Ghost cells are left untouched. Like
+[`set!`](@ref), `fun` runs on the device and the function comes first so the
+`do`-block form reads naturally.
+
+With `check=true` (the default) every field is first verified to live on the
+same grid as `ϕ` via [`check_compatible`](@ref); pass `check=false` only from a
+caller that has already checked. On forest fields the check runs once on the
+whole forest and the sweep then runs block by block.
+
+### Examples
+
+```julia
+g = CartesianGrid(((0.0, 1.0),), (64,))
+u = set!(x -> sin(x[1]), scalar_field(g))
+v = set!(x -> cos(x[1]), scalar_field(g))
+op!((x, a, b) -> a + x[1] * b, u, v)     # u ← u + x⋅v, cell by cell
+```
+
+See also: [`set!`](@ref), [`compatible`](@ref).
+"""
+function op!(f::F, ϕ::Field, ϕs::Field...; check::Bool=true) where {F}
+    check && check_compatible(ϕ, ϕs...)
+    g = AbstractGrid(ϕ)
+    map!(interior(ϕ), interior(g), interior(ϕ), map(interior, ϕs)...) do idx, a, b...
+        x = cell_center(g, idx)
+        f(x, a, b...)
+    end
+    return ϕ
 end
 
 """
@@ -285,3 +338,86 @@ function interior_to_flat!(v::AbstractVector, f::Field, α::Number=true, β::Num
     end
     return v
 end
+
+block(ϕ::Field, ::Integer, lg=ϕ.grid) = ϕ
+_block_array(ϕ::Field, ::Integer) = ϕ.data
+_require_current(::Field) = nothing
+
+# A `Field` is exactly as compatible as its grid.
+@inline _field_mismatch(a::Field, b::Field) = _grid_mismatch(a.grid, b.grid)
+@inline _field_layout_mismatch(a::Field, b::Field) = _layout_mismatch(a.grid, b.grid)
+
+# `Field` vs. forest field, or two forest fields of different storage layout
+# (`BlockField` vs. `PackedBlockField`): the per-block sweep can still pair them
+# through `_block_array`, so layout is *not* refused on the field type — only the
+# grid decides.  (Change these to `:type` if a function needs identical storage.)
+_field_mismatch(::AbstractField, ::AbstractField) = :type
+_field_layout_mismatch(::AbstractField, ::AbstractField) = :type
+
+"""
+    compatible(a::AbstractField, b::AbstractField...) -> Bool
+
+Whether every field lives on a grid `==` to `a`'s grid. Forest fields must
+also be current — allocated on the forest's present leaf set (see
+[`regrid!`](@ref)); a stale field is never compatible with anything.
+
+### Examples
+
+```julia
+g = CartesianGrid(((0.0, 1.0),), (64,))
+u = scalar_field(g); v = scalar_field(g)
+compatible(u, v)                       # true — same grid object, one `===`
+compatible(u, scalar_field(coarsen(g)))  # false
+```
+
+See also: [`check_compatible`](@ref), [`same_layout`](@ref).
+"""
+@inline compatible(::AbstractField) = true
+@inline compatible(a::AbstractField, b::AbstractField, rest::AbstractField...) =
+    isnothing(_field_mismatch(a, b)) && compatible(a, rest...)
+
+"""
+    check_compatible(a::AbstractField, b::AbstractField...) -> nothing
+    check_compatible(a::AbstractGrid, b::AbstractGrid...) -> nothing
+
+Throw an `ArgumentError` naming the first property on which any argument differs
+from `a` (grid `==`, plus regrid currency for forest fields).
+Return `nothing` otherwise. The guard for a multi-field function:
+
+```julia
+function fma!(y::AbstractField, α, x::AbstractField, z::AbstractField)
+    check_compatible(y, x, z)
+    ...
+end
+```
+
+The identity fast path makes this free when all fields share one grid object;
+the throw is out of line, so the check inlines into the caller.
+
+See also: [`check_layout`](@ref), [`compatible`](@ref).
+"""
+# Recursion over the argument tuple rather than `foreach` with a closure: the
+# closure form allocates ~1 KB per call from three fields up under
+# `--check-bounds=yes`; the recursive form measures 0 B for Field and BlockField.
+@inline check_compatible(::AbstractField) = nothing
+@inline function check_compatible(a::AbstractField, b::AbstractField, rest::AbstractField...)
+    _check_pair(_field_mismatch(a, b), a, b, "grid")
+    return check_compatible(a, rest...)
+end
+
+"""
+    check_layout(a::AbstractField, b::AbstractField...) -> nothing
+    check_layout(a::AbstractGrid, b::AbstractGrid...) -> nothing
+
+The [`same_layout`](@ref) counterpart of [`check_compatible`](@ref): throw an
+`ArgumentError` unless every argument has `a`'s padded storage shape (and, for
+forest fields, is current). For pointwise kernels that never read spacing or
+boundary conditions.
+"""
+@inline check_layout(::AbstractField) = nothing
+@inline function check_layout(a::AbstractField, b::AbstractField, rest::AbstractField...)
+    _check_pair(_field_layout_mismatch(a, b), a, b, "layout")
+    return check_layout(a, rest...)
+end
+
+nleaves(ϕ::Field) = nleaves(AbstractGrid(ϕ))
diff --git a/src/Grids.jl b/src/Grids.jl
index e04ef50..c8a130d 100644
--- a/src/Grids.jl
+++ b/src/Grids.jl
@@ -84,6 +84,8 @@ function _validate_bc(bc::Tuple, ::Val{N}) where {N}
     return nothing
 end
 
+AbstractGrid(g::AbstractGrid) = g
+
 #--------------------------------------------------------------------------------# Grid interface
 
 """
@@ -205,9 +207,124 @@ call this before any stencil that reads neighbor cells.
 """
 halo_update!(x, ::AbstractGrid) = x
 
+nleaves(::CartesianGrid) = 1
+function leaf_grid(g::CartesianGrid, j::Integer)
+    j == 1 || throw(ArgumentError("CartesianGrid only has a single leaf"))
+    g
+end
+
 function Adapt.adapt_structure(to, g::CartesianGrid{N}) where {N}
     device = KernelAbstractions.get_backend(Adapt.adapt(to, similar(Vector{Bool}, 0)))
     return CartesianGrid{N,eltype(g.spacing),typeof(g.bc),typeof(device),typeof(g.topology)}(
         g.extent, g.spacing, g.size, g.halo, g.bc, device, g.local_range, g.topology
     )
 end
+
+
+@inline function _grid_mismatch(a::CartesianGrid{N}, b::CartesianGrid{N}) where {N}
+    a === b && return nothing
+    a.size == b.size || return :size
+    a.halo == b.halo || return :halo
+    a.spacing == b.spacing || return :spacing
+    a.extent == b.extent || return :extent
+    a.bc === b.bc || return :bc
+    a.local_range == b.local_range || return :local_range
+    a.topology === b.topology || return :topology
+    return nothing
+end
+
+# Different grid types, or the same type in different dimensions: never compatible.
+_grid_mismatch(::AbstractGrid, ::AbstractGrid) = :type
+
+@inline function _layout_mismatch(a::CartesianGrid{N}, b::CartesianGrid{N}) where {N}
+    a === b && return nothing
+    a.size == b.size || return :size
+    a.halo == b.halo || return :halo
+    return nothing
+end
+
+_layout_mismatch(::AbstractGrid, ::AbstractGrid) = :type
+
+"""
+    a::AbstractGrid == b::AbstractGrid -> Bool
+
+Whether `a` and `b` describe the same discretization: equal cell counts, halo,
+spacing, extent, boundary conditions, ownership (`local_range`) and distributed
+topology for [`CartesianGrid`](@ref)s; equal block size, halo, root spacing,
+extent, physical boundary conditions and forest topology (root tiling,
+periodicity, `maxlevel`, leaf set) for [`BlockForest`](@ref)s. Grids of
+different types or dimensions are never equal. The `device` is deliberately
+**not** compared: a host grid and its `Adapt`-ed twin are the same
+discretization.
+
+`a === b` short-circuits, so the common case — two fields allocated from one
+grid object — costs a single comparison. `hash` is not specialized to match, so
+grids that are `==` but distinct objects are not interchangeable `Dict` keys.
+
+See also: [`same_layout`](@ref), [`check_compatible`](@ref).
+"""
+Base.:(==)(a::AbstractGrid, b::AbstractGrid) = isnothing(_grid_mismatch(a, b))
+
+@inline function _check_pair(what, a, b, level::String)
+    what === nothing || _throw_incompatible(what, a, b, level)
+    return nothing
+end
+
+# Grid-level twins of the field checks; documented with them in Fields.jl.
+# Tuple recursion, not `foreach` + closure — see the note there.
+@inline check_compatible(::AbstractGrid) = nothing
+@inline function check_compatible(a::AbstractGrid, b::AbstractGrid, rest::AbstractGrid...)
+    _check_pair(_grid_mismatch(a, b), a, b, "grid")
+    return check_compatible(a, rest...)
+end
+@inline check_layout(::AbstractGrid) = nothing
+@inline function check_layout(a::AbstractGrid, b::AbstractGrid, rest::AbstractGrid...)
+    _check_pair(_layout_mismatch(a, b), a, b, "layout")
+    return check_layout(a, rest...)
+end
+
+"""
+    same_layout(a::AbstractGrid, b::AbstractGrid) -> Bool
+
+Whether fields on `a` and `b` have the same padded storage shape, so that a
+pointwise broadcast between them (block by block for forests) is legal. Weaker
+than grid `==`: says nothing about geometry or boundary conditions.
+The right check for kernels that only ever combine values cell by cell.
+
+See also: [`check_layout`](@ref), [`check_compatible`](@ref).
+"""
+same_layout(a::AbstractGrid, b::AbstractGrid) = isnothing(_layout_mismatch(a, b))
+
+
+# Why each property matters, for an error that names the consequence rather than
+# just the field name.
+const _MISMATCH_REASONS = Dict{Symbol,String}(
+    :type => "they are different grid types or dimensions",
+    :size => "their interior cell counts differ, so the padded arrays cannot be broadcast together",
+    :halo => "their ghost-layer widths differ, so the same padded index means different cells",
+    :spacing => "their cell spacings differ",
+    :extent => "they cover different physical domains",
+    :bc => "their boundary conditions differ, so ghost cells would mean different things",
+    :local_range => "they own different global index ranges (a slab vs. the grid it was cut from?)",
+    :topology => "they belong to different distributed topologies",
+    :blocksize => "their block sizes differ",
+    :spacing0 => "their root-level spacings differ",
+    :nroot => "their root tilings differ",
+    :periodic => "their periodicity differs",
+    :maxlevel => "their maximum refinement levels differ",
+    :leaves => "their leaf sets differ (same domain, different refinement)",
+    :nleaves => "they have different numbers of leaf blocks",
+    :stale => "a forest field was allocated before its forest was regridded; allocate a fresh field on the current forest",
+)
+
+@noinline function _throw_incompatible(what::Symbol, a, b, level::String)
+    reason = get(_MISMATCH_REASONS, what, "they differ in $what")
+    throw(
+        ArgumentError(
+            "incompatible $level: $reason. Got " *
+            "$(summary(AbstractGrid(a))) vs $(summary(AbstractGrid(b))) [$what]",
+        ),
+    )
+end
+
+
diff --git a/src/MatrixFreeOperators.jl b/src/MatrixFreeOperators.jl
index 97c6a42..9d82098 100644
--- a/src/MatrixFreeOperators.jl
+++ b/src/MatrixFreeOperators.jl
@@ -15,6 +15,7 @@ export refine!, coarsen!, balance!, leaves, regrid!
 export Interpolated, Conservative, SlopeLimited, with_transfer
 export AbstractField, Field, Center, scalar_field, vector_field, set!, ncomponents, component
 export flatten, flat_to_interior!, interior_to_flat!
+export op!, compatible, check_compatible, check_layout, same_layout
 export AbstractOperator, apply, apply!, apply_adjoint!, AdjointOp
 export islinear, isconstant, isselfadjoint, isdiagonal, shares_exchange, operator_diagonal
 export Laplacian, laplacian, laplacian_stencil, laplacian_7pt_noflux
diff --git a/src/blockfield.jl b/src/blockfield.jl
index 119c2b1..ec09c96 100644
--- a/src/blockfield.jl
+++ b/src/blockfield.jl
@@ -187,11 +187,19 @@ end
 Base.copy(f::BlockField{L,P}) where {L,P} =
     BlockField{L,P}([copy(b) for b in f.blocks], f.grid, f.generation)
 
-function set!(f::AbstractBlockField, fun::F) where {F}
-    for i in 1:nleaves(f.grid)
-        set!(block(f, i), fun)
+function set!(f::F, ϕ::AbstractBlockField) where {F}
+    for i in 1:nleaves(AbstractGrid(ϕ))
+        set!(f, block(ϕ, i))
     end
-    return f
+    return ϕ
+end
+
+function op!(f::F, ϕ::AbstractBlockField, ϕs::AbstractBlockField...; check::Bool=true) where {F}
+    check && check_compatible(ϕ, ϕs...)
+    for i ∈ 1:nleaves(AbstractGrid(ϕ))
+        op!(f, block(ϕ, i), map(ζ -> block(ζ, i), ϕs)...; check=false)
+    end
+    return ϕ
 end
 
 function zero_ghosts!(f::AbstractBlockField)
@@ -251,3 +259,24 @@ function _interior_to_flat_leaves!(v::AbstractVector, f::AbstractBlockField, α:
     end
     return v
 end
+
+# A forest field additionally carries the regrid generation its storage was
+# allocated on.  Two fields that are both current on compatible forests agree on
+# the leaf set by construction, so currency is the whole check — comparing the
+# generation *numbers* would be wrong across two independently built forests
+# (one regridded three times, one built fresh from the same leaf set).
+@inline _is_current(f::AbstractBlockField) = f.generation == f.grid.forest.generation[]
+
+@inline function _field_mismatch(a::AbstractBlockField, b::AbstractBlockField)
+    _is_current(a) || return :stale
+    _is_current(b) || return :stale
+    return _grid_mismatch(a.grid, b.grid)
+end
+
+@inline function _field_layout_mismatch(a::AbstractBlockField, b::AbstractBlockField)
+    _is_current(a) || return :stale
+    _is_current(b) || return :stale
+    return _layout_mismatch(a.grid, b.grid)
+end
+
+
diff --git a/src/distributed.jl b/src/distributed.jl
index 4b9b9ec..f355189 100644
--- a/src/distributed.jl
+++ b/src/distributed.jl
@@ -734,7 +734,7 @@ end
 # right-hand side independent of the partition count.
 function _dist_set!(fields, fun::F, ctx) where {F}
     _dist_map!(ctx) do p
-        set!(fields[p], fun)
+        set!(fun, fields[p])
     end
     return fields
 end
diff --git a/src/linalg.jl b/src/linalg.jl
index 0a79959..203ee8f 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -105,11 +105,11 @@ step at field level with [`apply!`](@ref) instead.
 ```julia
 g = CartesianGrid(((0.0, 1.0),), (64,))
 A = prepare(laplacian(g))
-b = flatten(set!(scalar_field(g), x -> sin(π * x[1])))
+b = flatten(set!(x -> sin(π * x[1]), scalar_field(g)))
 u, stats = Krylov.minres(A, b)
 
 # explicit stepping stays at field level:
-uf = set!(scalar_field(g), x -> sin(π * x[1]))
+uf = set!(x -> sin(π * x[1]), scalar_field(g))
 du = similar(uf)
 dt = 0.4 * spacing(g)[1]^2              # forward-Euler bound
 apply!(du, A, uf)                       # or apply!(du, laplacian(g), uf)
diff --git a/src/multigrid.jl b/src/multigrid.jl
index a148885..a789204 100644
--- a/src/multigrid.jl
+++ b/src/multigrid.jl
@@ -345,7 +345,7 @@ g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
 L = -laplacian(g)
 A = prepare(L)
 M = MultigridPreconditioner(L)
-b = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
+b = flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
 u, stats = Krylov.cg(A, b; M)
 ```
 
diff --git a/src/operators/abstract.jl b/src/operators/abstract.jl
index eb5254c..c350c02 100644
--- a/src/operators/abstract.jl
+++ b/src/operators/abstract.jl
@@ -136,7 +136,7 @@ use in-place [`apply!`](@ref) instead (explicit stepping), or [`prepare`](@ref)
 
 ```julia
 g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
-u = set!(scalar_field(g), x -> sin(x[1]))
+u = set!(x -> sin(x[1]), scalar_field(g))
 Δu = apply(laplacian(g), u)        # equivalently laplacian(g)(u) or laplacian(g) * u
 ```
 """
diff --git a/src/operators/advection.jl b/src/operators/advection.jl
index 497845e..6b71871 100644
--- a/src/operators/advection.jl
+++ b/src/operators/advection.jl
@@ -50,7 +50,7 @@ nonlinear — see [`linearize`](@ref) for its Jacobian.
 ```julia
 g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                   bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
-v = set!(vector_field(g), x -> SVector(1.0, 0.0))
+v = set!(x -> SVector(1.0, 0.0), vector_field(g))
 A = advection(g, v)                  # linear: passive transport by v
 B = advection(g, SelfAdvection())    # nonlinear: u·∇u
 ```
diff --git a/src/operators/derivative.jl b/src/operators/derivative.jl
index 6f49ea1..90e3699 100644
--- a/src/operators/derivative.jl
+++ b/src/operators/derivative.jl
@@ -67,7 +67,7 @@ componentwise on any element type.
 
 ```julia
 g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
-u = set!(scalar_field(g), x -> sin(x[1]))
+u = set!(x -> sin(x[1]), scalar_field(g))
 ∂u = derivative(g, 1) * u        # ≈ cos
 ```
 
diff --git a/src/operators/diffusion.jl b/src/operators/diffusion.jl
index ad87e61..dd8f090 100644
--- a/src/operators/diffusion.jl
+++ b/src/operators/diffusion.jl
@@ -518,7 +518,7 @@ object.
 
 ```julia
 g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
-κ = set!(scalar_field(g), x -> 1 + x[1]^2)
+κ = set!(x -> 1 + x[1]^2, scalar_field(g))
 L = diffusion(g, κ)                                  # ∇·(κ∇u), compact 5-point
 Lj = diffusion(g, κ; averaging=HarmonicMean())       # flux-conserving across jumps
 ```
@@ -598,7 +598,7 @@ is unavailable on forest leaves, matching [`laplacian`](@ref).
 
 ```julia
 bf = BlockForest(CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16)); blocksize=(8, 8), maxlevel=2)
-κ = set!(scalar_field(bf), x -> 1 + x[1]^2)
+κ = set!(x -> 1 + x[1]^2, scalar_field(bf))
 L = diffusion(bf, κ)
 ```
 """
diff --git a/src/operators/divergence.jl b/src/operators/divergence.jl
index f2dab6f..08eef20 100644
--- a/src/operators/divergence.jl
+++ b/src/operators/divergence.jl
@@ -48,7 +48,7 @@ Variable-coefficient diffusion can be assembled as
 ```julia
 g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                   bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
-v = set!(vector_field(g), x -> SVector(sin(x[1]), cos(x[2])))
+v = set!(x -> SVector(sin(x[1]), cos(x[2])), vector_field(g))
 divv = divergence(g) * v
 ```
 
diff --git a/src/operators/gradient.jl b/src/operators/gradient.jl
index 9f20866..18a4830 100644
--- a/src/operators/gradient.jl
+++ b/src/operators/gradient.jl
@@ -49,7 +49,7 @@ vector field — one of the two rank-changing leaves (with [`divergence`](@ref))
 ```julia
 g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                   bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
-u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
+u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
 ∇u = gradient(g) * u             # SVector{2}-valued Field
 ```
 
diff --git a/src/operators/laplacian.jl b/src/operators/laplacian.jl
index afdc031..59dc793 100644
--- a/src/operators/laplacian.jl
+++ b/src/operators/laplacian.jl
@@ -84,7 +84,7 @@ under the homogeneous ghost fills of all built-in boundary conditions.
 
 ```julia
 g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
-u = set!(scalar_field(g), x -> sin(x[1]))
+u = set!(x -> sin(x[1]), scalar_field(g))
 Δu = laplacian(g) * u            # ≈ -u
 ```
 
diff --git a/src/operators/linearize.jl b/src/operators/linearize.jl
index 7f5b001..274eb6b 100644
--- a/src/operators/linearize.jl
+++ b/src/operators/linearize.jl
@@ -70,7 +70,7 @@ The default [`FiniteDifferenceJVP`](@ref) keeps the core dependency-free;
 ```julia
 g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
 F = advection(g, SelfAdvection())            # nonlinear u·∇u
-u0 = set!(vector_field(g), x -> SVector(sin(x[1])))
+u0 = set!(x -> SVector(sin(x[1])), vector_field(g))
 J = linearize(F, u0)                         # linear: v ↦ (∂F/∂u)|_{u0} · v
 P = prepare(J, u0)                           # Krylov-ready JFNK Jacobian
 
diff --git a/src/operators/scaling.jl b/src/operators/scaling.jl
index 3fdda95..7253053 100644
--- a/src/operators/scaling.jl
+++ b/src/operators/scaling.jl
@@ -36,7 +36,7 @@ the natural Jacobi-smoother target.
 
 ```julia
 g = CartesianGrid(((0.0, 1.0),), (64,))
-κ = set!(scalar_field(g), x -> 1 + x[1]^2)
+κ = set!(x -> 1 + x[1]^2, scalar_field(g))
 H = laplacian(g) - scaling(κ)                    # Helmholtz-type: Δu − κu
 ```
 
diff --git a/src/packedfield.jl b/src/packedfield.jl
index cac1507..e852a68 100644
--- a/src/packedfield.jl
+++ b/src/packedfield.jl
@@ -156,7 +156,7 @@ remaining layout mismatch degrades to the per-leaf reference sweep.
 ### Examples
 
 ```julia
-u  = set!(scalar_field(bf), x -> sin(π * x[1]))
+u  = set!(x -> sin(π * x[1]), scalar_field(bf))
 P  = prepare(laplacian(bf), pack(u))
 ```
 
diff --git a/test/algebra.jl b/test/algebra.jl
index f49cf52..7a0d88c 100644
--- a/test/algebra.jl
+++ b/test/algebra.jl
@@ -1,7 +1,7 @@
 @testset "Operator algebra" begin
     g = periodic_grid_2d(24)
-    u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
-    v = set!(vector_field(g), x -> SVector(sin(x[2]), cos(x[1])))
+    u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
+    v = set!(x -> SVector(sin(x[2]), cos(x[1])), vector_field(g))
     A = laplacian(g)
     B = advection(g, v)
     ints(f) = collect(interior(f))
@@ -55,15 +55,15 @@
     @testset "variable-coefficient diffusion ∇·(κ∇u) stress test" begin
         function vc_error(n)
             gn = periodic_grid_2d(n)
-            κ = set!(scalar_field(gn), x -> 2 + cos(x[1]))
-            un = set!(scalar_field(gn), x -> sin(x[1]) * sin(x[2]))
+            κ = set!(x -> 2 + cos(x[1]), scalar_field(gn))
+            un = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(gn))
             K = divergence(gn) * scaling(κ) * MatrixFreeOperators.gradient(gn)
             y = K * un
             ref = set!(
-                scalar_field(gn),
                 x ->
                     -2 * (2 + cos(x[1])) * sin(x[1]) * sin(x[2]) -
                     sin(x[1]) * cos(x[1]) * sin(x[2]),
+                scalar_field(gn),
             )
             return maximum(abs, collect(interior(y)) .- collect(interior(ref)))
         end
@@ -75,7 +75,7 @@
             ((0.0, 1.0), (0.0, 1.0)), (4, 3);
             bc=((Dirichlet(), Dirichlet()), (Neumann(), Neumann())),
         )
-        κ = set!(scalar_field(gn), x -> 1 + x[1] * x[2])
+        κ = set!(x -> 1 + x[1] * x[2], scalar_field(gn))
         K = divergence(gn) * scaling(κ) * MatrixFreeOperators.gradient(gn)
         M = materialize(prepare(K, scalar_field(gn)))
         Mt = materialize(prepare(adjoint(K), scalar_field(gn)))
@@ -87,7 +87,7 @@
             ((0.0, 1.0), (0.0, 1.0)), (4, 4);
             bc=((Periodic(), Periodic()), (Dirichlet(), Neumann())),
         )
-        vn = set!(vector_field(gn), x -> SVector(1 + x[1], x[2] - 2))
+        vn = set!(x -> SVector(1 + x[1], x[2] - 2), vector_field(gn))
         Adv = advection(gn, vn)
         M = materialize(prepare(Adv, scalar_field(gn)))
         Mt = materialize(prepare(adjoint(Adv), scalar_field(gn)))
diff --git a/test/amr_driver.jl b/test/amr_driver.jl
index 1241168..e3210e2 100644
--- a/test/amr_driver.jl
+++ b/test/amr_driver.jl
@@ -14,7 +14,7 @@
             key = bf.forest.leaves[i]
             sp = MFO._leaf_spacing(bf, key.level)
             lg = MFO.leaf_grid(bf, i)
-            ref = set!(scalar_field(lg), fun)
+            ref = set!(fun, scalar_field(lg))
             r = collect(interior(MFO.block(u, i))) .- collect(interior(ref))
             e += sum(abs, r) * prod(sp)
         end
@@ -25,7 +25,7 @@
     function max_error(u, bf, fun)
         maximum(1:MFO.nleaves(bf)) do i
             lg = MFO.leaf_grid(bf, i)
-            ref = set!(similar(MFO.block(u, i)), fun)
+            ref = set!(fun, similar(MFO.block(u, i)))
             maximum(abs, interior(MFO.block(u, i)) .- interior(ref))
         end
     end
@@ -33,7 +33,7 @@
     @testset "marking and topology semantics" begin
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=dirbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)   # 4×4 roots
-        u = set!(scalar_field(bf), x -> x[1])
+        u = set!(x -> x[1], scalar_field(bf))
 
         # refine the right column (blocks whose data exceeds the threshold)
         u = regrid!(u; refine=b -> maximum(interior(b)) > 0.75)
@@ -69,7 +69,7 @@
     @testset "no-op regrid preserves fields and prepared operators" begin
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(scalar_field(bf), lin)
+        u = set!(lin, scalar_field(bf))
         w = scalar_field(bf)
         P = prepare(laplacian(bf), u)
         gen = bf.forest.generation[]
@@ -84,7 +84,7 @@
         for bc in (dirbc, mixbc, perbc)
             base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-            u = set!(scalar_field(bf), lin)
+            u = set!(lin, scalar_field(bf))
 
             # refine a band touching the physical boundary: one-sided edge taps and
             # boundary-adjacent blocks are exercised for every BC kind
@@ -102,7 +102,7 @@
     @testset "balance-cascade leaves are transferred (key-based, not mark-based)" begin
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=dirbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-        u = set!(scalar_field(bf), lin)
+        u = set!(lin, scalar_field(bf))
         u = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.25 && b.grid.extent[2][1] < 0.25)
         gen1 = bf.forest.generation[]
         # corner → level 2 forces balance! to refine unmarked neighbors to level 1
@@ -118,20 +118,20 @@
         vfun = x -> SVector(x[1] - 2x[2], 1 + x[2])
 
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(scalar_field(bf), lin)
-        w = set!(vector_field(bf), vfun)
+        u = set!(lin, scalar_field(bf))
+        w = set!(vfun, vector_field(bf))
         band = b -> b.grid.extent[1][1] < 0.5
         u2, w2 = regrid!(u, w; refine=band)
         @test u2 isa BlockField && w2 isa BlockField
         werr = maximum(1:MFO.nleaves(bf)) do i
-            ref = set!(similar(MFO.block(w2, i)), vfun)
+            ref = set!(vfun, similar(MFO.block(w2, i)))
             maximum(norm, interior(MFO.block(w2, i)) .- interior(ref))
         end
         @test werr < 1e-13
 
         # the same regrid of a lone field gives bit-identical blocks
         bf3 = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u3 = set!(scalar_field(bf3), lin)
+        u3 = set!(lin, scalar_field(bf3))
         u3 = regrid!(u3; refine=band)
         @test all(i -> u3.blocks[i] == u2.blocks[i], 1:MFO.nleaves(bf3))
     end
@@ -141,7 +141,7 @@
         function refine_transfer_error(ncells)
             base = CartesianGrid(((0.0, 2π), (0.0, 2π)), (ncells, ncells); bc=perbc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-            u = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(bf))
             u = regrid!(u; refine=b -> b.grid.extent[1][1] < 1.6)
             @test !bf.forest.uniform[]
             return l1_error(u, bf, fun)
@@ -149,9 +149,9 @@
         function coarsen_transfer_error(ncells)
             base = CartesianGrid(((0.0, 2π), (0.0, 2π)), (ncells, ncells); bc=perbc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-            ind = set!(scalar_field(bf), x -> x[1] < 1.6 ? 1.0 : 0.0)
+            ind = set!(x -> x[1] < 1.6 ? 1.0 : 0.0, scalar_field(bf))
             regrid!(ind; refine=b -> maximum(interior(b)) > 0.5)
-            u = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(bf))
             u = regrid!(u; refine=Returns(false), coarsen=Returns(true))
             @test bf.forest.uniform[]
             return l1_error(u, bf, fun)
@@ -167,7 +167,7 @@
     @testset "staleness and argument errors" begin
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=dirbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(scalar_field(bf), lin)
+        u = set!(lin, scalar_field(bf))
         P = prepare(laplacian(bf), u)
         v = flatten(u)
         u2 = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.5)
@@ -188,9 +188,9 @@
         rng = Random.MersenneTwister(29)
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=mixbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
-        ind = set!(scalar_field(bf), x -> x[1] < 0.5 && x[2] < 0.5 ? 1.0 : 0.0)
+        ind = set!(x -> x[1] < 0.5 && x[2] < 0.5 ? 1.0 : 0.0, scalar_field(bf))
         ind = regrid!(ind; refine=b -> maximum(interior(b)) > 0.5)
-        ind2 = set!(scalar_field(bf), x -> x[1] < 0.2 && x[2] < 0.2 ? 1.0 : 0.0)
+        ind2 = set!(x -> x[1] < 0.2 && x[2] < 0.2 ? 1.0 : 0.0, scalar_field(bf))
         regrid!(ind2; refine=b -> maximum(interior(b)) > 0.5)    # levels 0–2
         @test maximum(k -> k.level, bf.forest.leaves) == 2
 
@@ -245,10 +245,10 @@
         function l1_action_error(ncells)
             base = CartesianGrid(((0.0, 2π), (0.0, 2π)), (ncells, ncells); bc=perbc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-            ind = set!(scalar_field(bf), x -> x[1] < 1.6 ? 1.0 : 0.0)
+            ind = set!(x -> x[1] < 1.6 ? 1.0 : 0.0, scalar_field(bf))
             regrid!(ind; refine=b -> maximum(interior(b)) > 0.5)
             @test !bf.forest.uniform[]
-            u = set!(scalar_field(bf), x -> sin(x[1]) * sin(x[2]))
+            u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(bf))
             Lu = laplacian(bf) * u
             e = 0.0
             for i in 1:MFO.nleaves(bf)
@@ -282,7 +282,7 @@
         errs = Float64[]
         for cycle in 1:3
             P = prepare(laplacian(bf), u)
-            rhs = .-flatten(set!(scalar_field(bf), f))   # Δu = −f
+            rhs = .-flatten(set!(f, scalar_field(bf)))   # Δu = −f
             sol, stats = Krylov.gmres(P, rhs; rtol=1e-10)
             @test stats.solved
             flat_to_interior!(u, sol)
@@ -310,7 +310,7 @@
     @testset "re-prepared mul! stays within the steady-state allocation budget" begin
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=mixbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(scalar_field(bf), lin)
+        u = set!(lin, scalar_field(bf))
         u = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.5)
         P = prepare(laplacian(bf), u)
         v = flatten(u)
@@ -326,7 +326,7 @@
     @testset "invalid transfer policies fail before regridding" begin
         base = CartesianGrid(((0.0, 1.0),), (8,))
         bf = BlockForest(base; blocksize=(4,), maxlevel=1)
-        u = set!(scalar_field(bf), x -> x[1])
+        u = set!(x -> x[1], scalar_field(bf))
         p = pack(u)
         generation = bf.forest.generation[]
         keys = copy(bf.forest.leaves)
@@ -354,7 +354,7 @@
     )
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(vector_field(bf; transfer=pol), x -> SVector(x[1], 2x[2]))
+        u = set!(x -> SVector(x[1], 2x[2]), vector_field(bf; transfer=pol))
         p = pack(u)
         generation = u.generation
 
@@ -405,8 +405,8 @@
         # Changing field units must not turn a sloped reconstruction into injection.
         base = CartesianGrid(((zero(T), one(T)),), (8,))
         bf = BlockForest(base; blocksize=(4,), maxlevel=1)
-        u = set!(scalar_field(bf, T; transfer=SlopeLimited()), x -> x[1])
-        v = set!(scalar_field(bf, T; transfer=SlopeLimited()), x -> tiny * x[1])
+        u = set!(x -> x[1], scalar_field(bf, T; transfer=SlopeLimited()))
+        v = set!(x -> tiny * x[1], scalar_field(bf, T; transfer=SlopeLimited()))
         u, v = regrid!(u, v; refine=Returns(true))
         @test isapprox(flatten(v) ./ tiny, flatten(u); rtol=8eps(T))
     end
@@ -433,7 +433,7 @@
         for (name, pol) in (("Conservative", Conservative()), ("SlopeLimited", SlopeLimited()))
             base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=mixbc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-            u = set!(scalar_field(bf; transfer=pol), curved)
+            u = set!(curved, scalar_field(bf; transfer=pol))
             m0, s0 = field_mass(u, bf)
             # refine a corner, then its inside — the second regrid forces balance!
             # to refine leaves the criteria never marked (issue #59's acceptance:
@@ -464,7 +464,7 @@
         # has teeth
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=mixbc)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-        u = set!(scalar_field(bf), curved)
+        u = set!(curved, scalar_field(bf))
         m0, s0 = field_mass(u, bf)
         u = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.25 && b.grid.extent[2][1] < 0.25)
         m1, _ = field_mass(u, bf)
@@ -496,7 +496,7 @@
         inner(b) = all(e -> e[2] < 0.26, b.grid.extent)
         for (base, bs, fun) in cases, pol in (Conservative(), SlopeLimited())
             bf = BlockForest(base; blocksize=bs, maxlevel=2)
-            u = set!(scalar_field(bf; transfer=pol), fun)
+            u = set!(fun, scalar_field(bf; transfer=pol))
             nl0 = MFO.nleaves(bf)
             m0, s0 = field_mass(u, bf)
             u = regrid!(u; refine=corner)
@@ -516,7 +516,7 @@
 
             # negative control in the same dimension
             bfi = BlockForest(base; blocksize=bs, maxlevel=2)
-            ui = set!(scalar_field(bfi), fun)
+            ui = set!(fun, scalar_field(bfi))
             mi0, si0 = field_mass(ui, bfi)
             ui = regrid!(ui; refine=corner)
             mi1, _ = field_mass(ui, bfi)
@@ -530,7 +530,7 @@
         # Conservative is exact on linears everywhere (centered and one-sided
         # slopes both reproduce a linear), like the default
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        u = set!(scalar_field(bf; transfer=Conservative()), lin)
+        u = set!(lin, scalar_field(bf; transfer=Conservative()))
         u = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.5)
         @test max_error(u, bf, lin) < 1e-13
 
@@ -538,7 +538,7 @@
         # source's global bounds (minmod interior, zero slope at block edges)
         bf2 = BlockForest(base; blocksize=(4, 4), maxlevel=1)
         step = x -> x[1] < 0.4 ? 0.0 : 1.0
-        v = set!(scalar_field(bf2; transfer=SlopeLimited()), step)
+        v = set!(step, scalar_field(bf2; transfer=SlopeLimited()))
         lo, hi = extrema(
             reduce(vcat, [vec(collect(interior(MFO.block(v, i)))) for i in 1:MFO.nleaves(bf2)])
         )
@@ -549,15 +549,15 @@
         m_step, s_step = field_mass(v, bf2)   # and it still conserved
         # (mass computed post-refine equals pre-refine: recompute the reference)
         bf2b = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-        v0 = set!(scalar_field(bf2b), step)
+        v0 = set!(step, scalar_field(bf2b))
         m0_step, _ = field_mass(v0, bf2b)
         @test abs(m_step - m0_step) ≤ 1e3 * eps() * s_step
 
         # SVector state conserves componentwise
         bf3 = BlockForest(base; blocksize=(4, 4), maxlevel=1)
         w = set!(
-            vector_field(bf3; transfer=Conservative()),
             x -> SVector(sin(3 * x[1]) + x[2]^2, cos(2 * x[2]) - x[1]^2),
+            vector_field(bf3; transfer=Conservative()),
         )
         msum(f, bfx) = sum(
             i -> prod(MFO._leaf_spacing(bfx, bfx.forest.leaves[i].level)) .*
@@ -586,7 +586,7 @@
         function refine_l1(ncells)
             base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (ncells, ncells); bc=perbc)
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=1)
-            u = set!(scalar_field(bf; transfer=Conservative()), smooth)
+            u = set!(smooth, scalar_field(bf; transfer=Conservative()))
             u = regrid!(u; refine=b -> b.grid.extent[1][1] < 0.5)
             return l1_error(u, bf, smooth)
         end
diff --git a/test/autodiff.jl b/test/autodiff.jl
index e3d7118..b5282a1 100644
--- a/test/autodiff.jl
+++ b/test/autodiff.jl
@@ -47,7 +47,7 @@ end
     )
     rng = Random.MersenneTwister(61)
     w = rand(rng, local_size(g)...)
-    κ0 = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+    κ0 = set!(x -> 1 + x[1] * x[2], scalar_field(g))
 
     @testset "field gradient of $(name)" for (name, L) in (
         ("Laplacian", laplacian(g)),
@@ -87,7 +87,7 @@ end
     end
 
     @testset "parameter gradient w.r.t. coefficient field κ (Decision B)" begin
-        u = set!(scalar_field(g), x -> sin(3 * x[1]) * x[2])
+        u = set!(x -> sin(3 * x[1]) * x[2], scalar_field(g))
         κ = 1.0 .+ rand(rng, padded_size(g)...)
         dκ = zero(κ)
         Enzyme.autodiff(
@@ -112,7 +112,7 @@ end
                                                                                           (
         ArithmeticMean(), HarmonicMean()
     )
-        u = set!(scalar_field(g), x -> sin(3 * x[1]) * x[2])
+        u = set!(x -> sin(3 * x[1]) * x[2], scalar_field(g))
         κ = 1.0 .+ rand(rng, padded_size(g)...)
         dκ = zero(κ)
         # set_runtime_activity, mirror image of the field-gradient case above: the leaf
@@ -201,7 +201,7 @@ end
         n = length(flatten(scalar_field(bf)))
         v = rand(rng, n)
         wf = rand(rng, n)
-        D = diffusion(bf, set!(scalar_field(bf), x -> 1.2 + 0.8 * x[1]^2 + 0.5 * x[2]);
+        D = diffusion(bf, set!(x -> 1.2 + 0.8 * x[1]^2 + 0.5 * x[2], scalar_field(bf));
                       averaging=avg)
 
         # Field gradient: Enzyme's taped derivative of the coarse-ghost flux rewrite
@@ -256,7 +256,7 @@ end
     @testset "gradient through the nonlinear leaf u·∇u" begin
         g1 = CartesianGrid(((0.0, 2π),), (16,); bc=((Periodic(), Periodic()),))
         F = advection(g1, SelfAdvection())
-        u = set!(vector_field(g1), x -> SVector(2 + sin(x[1])))
+        u = set!(x -> SVector(2 + sin(x[1])), vector_field(g1))
         wv = [SVector(rand(rng)) for _ in 1:16]
 
         du = zero(u.data)
diff --git a/test/blockfield.jl b/test/blockfield.jl
index ce725b7..56e94e1 100644
--- a/test/blockfield.jl
+++ b/test/blockfield.jl
@@ -23,9 +23,9 @@
     @testset "set! matches per-leaf standalone field" begin
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
         f = scalar_field(bf)
-        set!(f, x -> x[1] + 2x[2])
+        set!(x -> x[1] + 2x[2], f)
         for i in 1:MFO.nleaves(bf)
-            ref = set!(scalar_field(MFO.leaf_grid(bf, i)), x -> x[1] + 2x[2])
+            ref = set!(x -> x[1] + 2x[2], scalar_field(MFO.leaf_grid(bf, i)))
             @test collect(interior(MFO.block(f, i))) ≈ collect(interior(ref))
         end
     end
@@ -33,7 +33,7 @@
     @testset "flat ↔ interior round-trip (uniform)" begin
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
         f = scalar_field(bf)
-        set!(f, x -> sin(3x[1]) * cos(2x[2]))
+        set!(x -> sin(3x[1]) * cos(2x[2]), f)
         v = flatten(f)
         @test length(v) == MFO.flat_length(f) == 4 * 16   # nleaves·prod(blocksize)
         g = scalar_field(bf)
@@ -46,7 +46,7 @@
         refine!(bf, x -> x[1] < 0.4 && x[2] < 0.4)         # refine a corner region
         refine!(bf, x -> x[1] < 0.2 && x[2] < 0.2)
         f = scalar_field(bf)
-        set!(f, x -> x[1]^2 - x[2])
+        set!(x -> x[1]^2 - x[2], f)
         v = flatten(f)
         @test length(v) == MFO.flat_length(f) == MFO.nleaves(bf) * 16
         g = scalar_field(bf)
@@ -57,7 +57,7 @@
     @testset "interior_to_flat! axpby" begin
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
         f = scalar_field(bf)
-        set!(f, x -> x[1] - x[2])
+        set!(x -> x[1] - x[2], f)
         v = flatten(f)
         v2 = copy(v)
         interior_to_flat!(v2, f, 2.0, 3.0)                 # 2·interior + 3·v2 = 5v
@@ -67,7 +67,7 @@
     @testset "vector field round-trip" begin
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
         vf = vector_field(bf)
-        set!(vf, x -> SVector(x[1], x[2]))
+        set!(x -> SVector(x[1], x[2]), vf)
         w = flatten(vf)
         @test length(w) == MFO.flat_length(vf) == 4 * 16 * 2
         g = vector_field(bf)
diff --git a/test/blockforest.jl b/test/blockforest.jl
index 3867f7c..b6a4a2c 100644
--- a/test/blockforest.jl
+++ b/test/blockforest.jl
@@ -94,7 +94,7 @@ struct UnsupportedBC <: MatrixFreeOperators.AbstractBC end
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
         refine!(bf, x -> x[1] < 0.5 && x[2] < 0.5)     # one corner root → mixed levels
         @test !bf.forest.uniform[]
-        uf = set!(scalar_field(bf), x -> x[1]^2 + x[2]^2)   # Δu = 4 exactly
+        uf = set!(x -> x[1]^2 + x[2]^2, scalar_field(bf))   # Δu = 4 exactly
         Lu = laplacian(bf) * uf
         # Every stencil fed only by interior/interface ghosts must be exact; skip
         # the one-cell layer whose stencil reads a homogeneous physical-BC ghost.
@@ -118,19 +118,19 @@ struct UnsupportedBC <: MatrixFreeOperators.AbstractBC end
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
         uf = scalar_field(bf)
         refine!(bf, _ -> false)                         # no-op regrid keeps fields valid
-        @test set!(uf, x -> x[1]) isa BlockField
+        @test set!(x -> x[1], uf) isa BlockField
         refine!(bf, _ -> true)                          # uniform level 1: new leaf set
         @test_throws ArgumentError laplacian(bf) * uf
         @test_throws ArgumentError flatten(uf)
-        @test_throws ArgumentError set!(uf, x -> x[1])
+        @test_throws ArgumentError set!(x -> x[1], uf)
         # copies and similars of a stale field are equally stale
-        @test_throws ArgumentError set!(copy(uf), x -> x[1])
+        @test_throws ArgumentError set!(x -> x[1], copy(uf))
         uf2 = scalar_field(bf)                          # fresh allocation works
-        @test flatten(set!(uf2, x -> x[1])) isa Vector
+        @test flatten(set!(x -> x[1], uf2)) isa Vector
         # coarsening back to the original leaf count is still a different generation
         coarsen!(bf, _ -> true)
         @test MFO.nleaves(bf) == 4
-        @test_throws ArgumentError set!(uf, x -> x[1])
-        @test_throws ArgumentError set!(uf2, x -> x[1])
+        @test_throws ArgumentError set!(x -> x[1], uf)
+        @test_throws ArgumentError set!(x -> x[1], uf2)
     end
 end
diff --git a/test/derivative.jl b/test/derivative.jl
index 03e1f01..d9e1dbe 100644
--- a/test/derivative.jl
+++ b/test/derivative.jl
@@ -1,8 +1,8 @@
 function deriv_periodic_error(n::Int, order::Int)
     g = CartesianGrid(((0.0, 2π),), (n,); bc=((Periodic(), Periodic()),))
-    u = set!(scalar_field(g), x -> sin(x[1]))
+    u = set!(x -> sin(x[1]), scalar_field(g))
     y = derivative(g, 1; order) * u
-    ref = set!(scalar_field(g), order == 1 ? (x -> cos(x[1])) : (x -> -sin(x[1])))
+    ref = set!(order == 1 ? (x -> cos(x[1])) : (x -> -sin(x[1])), scalar_field(g))
     return maximum(abs, collect(interior(y)) .- collect(interior(ref)))
 end
 
@@ -25,9 +25,9 @@ end
             ((0.0, 2π), (0.0, 2π)), (32, 32);
             bc=((Periodic(), Periodic()), (Periodic(), Periodic())),
         )
-        u = set!(scalar_field(g), x -> sin(x[2]))
+        u = set!(x -> sin(x[2]), scalar_field(g))
         dy = derivative(g, 2) * u
-        ref = set!(scalar_field(g), x -> cos(x[2]))
+        ref = set!(x -> cos(x[2]), scalar_field(g))
         @test maximum(abs, collect(interior(dy)) .- collect(interior(ref))) < 0.01
         dx = derivative(g, 1) * u
         @test maximum(abs, collect(interior(dx))) < 1e-12
diff --git a/test/device.jl b/test/device.jl
index 15ed554..de485d5 100644
--- a/test/device.jl
+++ b/test/device.jl
@@ -4,9 +4,9 @@
             ((0.0, 1.0), (0.0, 1.0)), (6, 5);
             bc=((Dirichlet(), Neumann()), (Periodic(), Periodic())),
         )
-        κ = set!(scalar_field(g), x -> 1 + x[1])
-        v = set!(vector_field(g), x -> SVector(x[1], 1.0))
-        u = set!(scalar_field(g), x -> sin(x[1]) * x[2])
+        κ = set!(x -> 1 + x[1], scalar_field(g))
+        v = set!(x -> SVector(x[1], 1.0), vector_field(g))
+        u = set!(x -> sin(x[1]) * x[2], scalar_field(g))
 
         K = divergence(g) * scaling(κ) * MatrixFreeOperators.gradient(g)
         ops = (
@@ -45,7 +45,7 @@
             bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic())),
         )
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
-        uf = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        uf = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         L = laplacian(bf)
 
         L2 = Adapt.adapt(Array, L)
diff --git a/test/device_gpu.jl b/test/device_gpu.jl
index d974d0c..de567a8 100644
--- a/test/device_gpu.jl
+++ b/test/device_gpu.jl
@@ -8,8 +8,8 @@ CUDA.allowscalar(false)
         ((0.0, 2π), (0.0, 1.0)), (32, 24);
         bc=((Periodic(), Periodic()), (Dirichlet(), Neumann())),
     )
-    κ = set!(scalar_field(g), x -> 1 + x[1] / 7)
-    u = set!(scalar_field(g), x -> sin(x[1]) * x[2])
+    κ = set!(x -> 1 + x[1] / 7, scalar_field(g))
+    u = set!(x -> sin(x[1]) * x[2], scalar_field(g))
 
     @testset "operator action parity" begin
         for L in (laplacian(g), derivative(g, 2), scaling(κ) - laplacian(g))
@@ -36,7 +36,7 @@ CUDA.allowscalar(false)
             bc=((Periodic(), Periodic()), (Dirichlet(), Dirichlet())),
         )
         bf = BlockForest(base; blocksize=(8, 8), maxlevel=2)
-        uf = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        uf = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         for L in (laplacian(bf), derivative(bf, 1; order=1))
             y_cpu = [
                 collect(interior(MFO.block(apply(L, copy(uf)), i))) for i in 1:MFO.nleaves(bf)
@@ -59,7 +59,7 @@ CUDA.allowscalar(false)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
         refine!(bf, x -> x[1] < π && x[2] < 0.5)       # mixed levels: CF ghost fills run
         @test !bf.forest.uniform[]
-        uf = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        uf = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         for L in (laplacian(bf), derivative(bf, 1; order=1))
             y_cpu = [
                 collect(interior(MFO.block(apply(L, copy(uf)), i))) for i in 1:MFO.nleaves(bf)
@@ -80,7 +80,7 @@ CUDA.allowscalar(false)
             bc=((Periodic(), Periodic()), (Dirichlet(), Dirichlet())),
         )
         bf = BlockForest(base; blocksize=(8, 8), maxlevel=2)
-        uf = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        uf = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         ug = Adapt.adapt(CuArray, uf)
         vg = flatten(ug)                       # device vector, no scalar indexing
         @test vg isa CuArray
@@ -102,7 +102,7 @@ CUDA.allowscalar(false)
             bc=((Periodic(), Periodic()), (Dirichlet(), Dirichlet())),
         )
         bf = BlockForest(base; blocksize=(8, 8), maxlevel=2)
-        uf = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        uf = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         p = pack(uf)
         pg = Adapt.adapt(CuArray, p)
         @test pg.data isa CuArray
@@ -127,7 +127,7 @@ CUDA.allowscalar(false)
 
         # refined forest: per-leaf levels SoA feeds the kernel on device
         refine!(bf, x -> x[1] < π)
-        ur = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        ur = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         pr = pack(ur)
         yr = apply(laplacian(bf), copy(pr))
         prg = Adapt.adapt(CuArray, pr)
@@ -146,10 +146,10 @@ CUDA.allowscalar(false)
         for refined in (false, true)
             bf = BlockForest(base; blocksize=(8, 8), maxlevel=2)
             refined && refine!(bf, x -> x[1] < π)
-            u = set!(scalar_field(bf), sfun)
-            w = set!(vector_field(bf), wfun)
-            κp = pack(set!(scalar_field(bf), x -> 1 + x[2]^2))
-            velp = pack(set!(vector_field(bf), wfun))
+            u = set!(sfun, scalar_field(bf))
+            w = set!(wfun, vector_field(bf))
+            κp = pack(set!(x -> 1 + x[2]^2, scalar_field(bf)))
+            velp = pack(set!(wfun, vector_field(bf)))
             p = pack(u)
             pw = pack(w)
             pg = Adapt.adapt(CuArray, p)
@@ -186,9 +186,9 @@ CUDA.allowscalar(false)
             @test Array(outg) ≈ out
 
             # declared adjoint transpose-gather kernels + fold path on device
-            ys = pack(set!(scalar_field(bf), x -> cos(x[1]) + x[2]^2))
+            ys = pack(set!(x -> cos(x[1]) + x[2]^2, scalar_field(bf)))
             for L in (derivative(bf, 1; order=1), MFO.gradient(bf))
-                ȳ = L isa MFO.Gradient ? pack(set!(vector_field(bf), wfun)) : ys
+                ȳ = L isa MFO.Gradient ? pack(set!(wfun, vector_field(bf))) : ys
                 x̄ = apply_adjoint!(MFO.allocate_input(L, ȳ), L, copy(ȳ), bf)
                 x̄g = apply_adjoint!(
                     MFO.allocate_input(Adapt.adapt(CuArray, L), Adapt.adapt(CuArray, ȳ)),
@@ -203,7 +203,7 @@ CUDA.allowscalar(false)
             # both averaging policies; on the refined forest the coarse–fine flux
             # rewrite runs as device-view broadcasts ahead of the launch
             for avg in (ArithmeticMean(), HarmonicMean())
-                D = diffusion(bf, set!(scalar_field(bf), x -> 1 + x[2]^2); averaging=avg)
+                D = diffusion(bf, set!(x -> 1 + x[2]^2, scalar_field(bf)); averaging=avg)
                 Dp = MFO.Diffusion(bf, pack(D.κ), D.avg)
                 Dg = Adapt.adapt(CuArray, Dp)
                 y = apply(Dp, copy(p))
@@ -229,7 +229,7 @@ CUDA.allowscalar(false)
         )
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
         refine!(bf, x -> x[1] < π)
-        u = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        u = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         p = pack(u)
         # host reference exchange + BC
         pr = copy(p)
@@ -263,7 +263,7 @@ CUDA.allowscalar(false)
         @test Array(outg2) ≈ out2
         # regrid invalidates the device schedule: fresh generation rebuilds
         refine!(bf, x -> x[2] < 0.5)
-        u2 = set!(scalar_field(bf), x -> sin(x[1]) * x[2])
+        u2 = set!(x -> sin(x[1]) * x[2], scalar_field(bf))
         pr2 = pack(u2)
         MFO.halo_update!(pr2, bf)
         pg3 = Adapt.adapt(CuArray, pack(u2))
@@ -277,9 +277,9 @@ CUDA.allowscalar(false)
     end
 
     @testset "Krylov cg parity" begin
-        σ = set!(scalar_field(g), x -> 1 + x[2])
+        σ = set!(x -> 1 + x[2], scalar_field(g))
         K = scaling(σ) - laplacian(g)
-        f = set!(scalar_field(g), x -> sin(x[1]))
+        f = set!(x -> sin(x[1]), scalar_field(g))
 
         P_cpu = prepare(K, scalar_field(g))
         b_cpu = flatten(f)
@@ -312,8 +312,8 @@ CUDA.allowscalar(false)
         # parity needs two independent forests; the GPU side adapts a twin whose
         # CPU original is used only to build the initial data.
         bf_cpu = mk()
-        u_cpu = set!(scalar_field(bf_cpu), bump)
-        u_gpu = Adapt.adapt(CuArray, set!(scalar_field(mk()), bump))
+        u_cpu = set!(bump, scalar_field(bf_cpu))
+        u_gpu = Adapt.adapt(CuArray, set!(bump, scalar_field(mk())))
         bf_gpu = u_gpu.grid
         @test first(u_gpu.blocks) isa CuArray
 
@@ -334,14 +334,14 @@ CUDA.allowscalar(false)
             flat_to_interior!(u, sol)
             return sol
         end
-        rhs = .-flatten(set!(scalar_field(bf_cpu), rhsf))    # identical topology ⇒ same layout
+        rhs = .-flatten(set!(rhsf, scalar_field(bf_cpu)))    # identical topology ⇒ same layout
         @test Array(cycle!(u_gpu, bf_gpu, CuArray(rhs))) ≈ cycle!(u_cpu, bf_cpu, rhs) rtol = 1e-6
 
         crit2 = b -> maximum(abs, interior(b)) > 0.5
         u_cpu = regrid!(u_cpu; refine=crit2)
         u_gpu = regrid!(u_gpu; refine=crit2)                 # prolongs solved data on device
         @test bf_gpu.forest.leaves == bf_cpu.forest.leaves
-        rhs2 = .-flatten(set!(scalar_field(bf_cpu), rhsf))
+        rhs2 = .-flatten(set!(rhsf, scalar_field(bf_cpu)))
         @test Array(cycle!(u_gpu, bf_gpu, CuArray(rhs2))) ≈ cycle!(u_cpu, bf_cpu, rhs2) rtol = 1e-6
 
         # coarsen everything back: the conservative child-mean path on device
diff --git a/test/differentiation_interface.jl b/test/differentiation_interface.jl
index b3e5603..1e5bf03 100644
--- a/test/differentiation_interface.jl
+++ b/test/differentiation_interface.jl
@@ -32,7 +32,7 @@ end
     )
     rng = Random.MersenneTwister(1907)
     w = rand(rng, local_size(g)...)
-    κ0 = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+    κ0 = set!(x -> 1 + x[1] * x[2], scalar_field(g))
 
     @testset "field gradient matches the declared adjoint: $(name)" for (name, L) in (
         ("Laplacian", laplacian(g)),
@@ -67,7 +67,7 @@ end
     end
 
     @testset "parameter gradient w.r.t. the coefficient field κ" begin
-        u = set!(scalar_field(g), x -> sin(3 * x[1]) * x[2])
+        u = set!(x -> sin(3 * x[1]) * x[2], scalar_field(g))
         κ = 1.0 .+ rand(rng, padded_size(g)...)
         dκ = DI.gradient(
             di_kappa_loss,
diff --git a/test/diffusion.jl b/test/diffusion.jl
index bcb8b45..4388bbe 100644
--- a/test/diffusion.jl
+++ b/test/diffusion.jl
@@ -40,10 +40,10 @@ end
                 ((0.0, 1.0), (0.0, 1.0)), (n, n);
                 bc=((Periodic(), Periodic()), (Periodic(), Periodic())),
             )
-            y = diffusion(g, set!(scalar_field(g), κf); averaging=avg) *
-                set!(scalar_field(g), uf)
+            y = diffusion(g, set!(κf, scalar_field(g)); averaging=avg) *
+                set!(uf, scalar_field(g))
             return maximum(
-                abs, collect(interior(y)) .- collect(interior(set!(scalar_field(g), exactf)))
+                abs, collect(interior(y)) .- collect(interior(set!(exactf, scalar_field(g))))
             )
         end
         @testset "$(nameof(typeof(avg)))" for avg in DIFF_AVGS
@@ -57,8 +57,8 @@ end
     @testset "1-D and 3-D action" begin
         # κ(x) = 1 + x, u(x) = x²  ⇒  ∇·(κ∇u) = (κu')' = (2x + 2x²)' = 2 + 4x
         g1 = CartesianGrid(((0.0, 1.0),), (128,); bc=((Neumann(), Neumann()),))
-        y1 = diffusion(g1, set!(scalar_field(g1), x -> 1 + x[1])) *
-             set!(scalar_field(g1), x -> x[1]^2)
+        y1 = diffusion(g1, set!(x -> 1 + x[1], scalar_field(g1))) *
+             set!(x -> x[1]^2, scalar_field(g1))
         # Interior only: the wall rows carry the homogeneous-Neumann flux, not (κu')'.
         yi = collect(interior(y1))[2:(end - 1)]
         xi = [cell_center(g1, I)[1] for I in interior(g1)][2:(end - 1)]
@@ -71,7 +71,7 @@ end
         )
         κ3 = scalar_field(g3)
         interior(κ3) .= 2.0
-        u3 = set!(scalar_field(g3), x -> prod(sin, x))
+        u3 = set!(x -> prod(sin, x), scalar_field(g3))
         @test collect(interior(diffusion(g3, κ3) * u3)) ≈
             2 .* collect(interior(laplacian(g3) * u3))
     end
@@ -315,7 +315,7 @@ end
         # so ∂(Lu)_I/∂κ_I ≡ 0. The compact form includes κ_I. This local
         # sensitivity does not prove that the full κ-Jacobian has no null space.
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (7, 7))
-        u = set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(2 * x[2]) + x[1] * x[2])
+        u = set!(x -> sinpi(x[1]) * sinpi(2 * x[2]) + x[1] * x[2], scalar_field(g))
         I0 = CartesianIndex(4, 4)                        # a strictly interior cell
 
         function dLu_dkappa(build, κdata, I)
@@ -452,9 +452,9 @@ end
         κ = scalar_field(g, Float32)
         interior(κ) .= 1.0f0
         D = diffusion(g, κ)
-        y = D * set!(scalar_field(g, Float32), x -> sinpi(x[1]))
+        y = D * set!(x -> sinpi(x[1]), scalar_field(g, Float32))
         @test eltype(y.data) === Float32
-        @test collect(interior(y)) ≈ collect(interior(laplacian(g) * set!(scalar_field(g, Float32), x -> sinpi(x[1])))) rtol =
+        @test collect(interior(y)) ≈ collect(interior(laplacian(g) * set!(x -> sinpi(x[1]), scalar_field(g, Float32)))) rtol =
             1.0f-5
         @test Adapt.adapt(Array, D) isa Diffusion
     end
@@ -480,8 +480,8 @@ end
     @testset "coefficient grid compatibility" begin
         g = CartesianGrid(((0.0, 1.0),), (6,))
         equivalent = CartesianGrid(((0.0, 1.0),), (6,))
-        κ = set!(scalar_field(equivalent), _ -> 2)
-        u = set!(scalar_field(g), x -> sinpi(x[1]))
+        κ = set!(_ -> 2, scalar_field(equivalent))
+        u = set!(x -> sinpi(x[1]), scalar_field(g))
         @testset "$(nameof(typeof(avg))) check=$check" for avg in DIFF_AVGS,
             check in (true, false)
 
@@ -502,7 +502,7 @@ end
                 @test padded_size(other) == padded_size(g)
                 # Positive interiors pass the harmonic check on `other`, while
                 # its zero ghosts must never become interior coefficients on g.
-                foreign = set!(scalar_field(other), _ -> 2)
+                foreign = set!(_ -> 2, scalar_field(other))
                 @test all(>(0), interior(foreign))
                 @test_throws ArgumentError diffusion(
                     g, foreign; averaging=avg, check=check
@@ -518,7 +518,7 @@ end
         )
         D = diffusion(g, diff_kappa(g))
         L = laplacian(g)
-        u = set!(scalar_field(g), x -> sinpi(2 * x[1]) * x[2])
+        u = set!(x -> sinpi(2 * x[1]) * x[2], scalar_field(g))
         ints(f) = collect(interior(f))
 
         @test ints((D + L) * copy(u)) ≈ ints(D * copy(u)) .+ ints(L * copy(u))
@@ -533,7 +533,7 @@ end
 
     @testset "multigrid rediscretization" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16))
-        κ = set!(scalar_field(g), x -> 1 + x[1]^2 + x[2])
+        κ = set!(x -> 1 + x[1]^2 + x[2], scalar_field(g))
         D = diffusion(g, κ; averaging=HarmonicMean())
         gc = coarsen(g)
         Dc = MatrixFreeOperators._rediscretize(D, gc)
diff --git a/test/enzyme_rules.jl b/test/enzyme_rules.jl
index 09e2d04..18734ef 100644
--- a/test/enzyme_rules.jl
+++ b/test/enzyme_rules.jl
@@ -140,7 +140,7 @@ end
         n = length(flatten(scalar_field(bf)))
         v = rand(rng, n)
         wf = rand(rng, n)
-        D = diffusion(bf, set!(scalar_field(bf), x -> 1.2 + 0.8 * x[1]^2 + 0.5 * x[2]))
+        D = diffusion(bf, set!(x -> 1.2 + 0.8 * x[1]^2 + 0.5 * x[2], scalar_field(bf)))
 
         # The measurable form of "no additional halo exchange per application": one
         # exchange-rule invocation per apply, exactly the Laplacian's count — the
@@ -253,7 +253,7 @@ end
             bc=((Dirichlet(), Dirichlet()), (Neumann(), Neumann())),
         )
         w = rand(rng, local_size(g)...)
-        κ0 = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+        κ0 = set!(x -> 1 + x[1] * x[2], scalar_field(g))
 
         @testset "$(name)" for (name, L) in (
             ("Laplacian", laplacian(g)),
@@ -291,7 +291,7 @@ end
             bc=((Dirichlet(), Dirichlet()), (Neumann(), Neumann())),
         )
         w = rand(rng, local_size(g)...)
-        u = set!(scalar_field(g), x -> sin(3 * x[1]) * x[2])
+        u = set!(x -> sin(3 * x[1]) * x[2], scalar_field(g))
         κ = 1.0 .+ rand(rng, padded_size(g)...)
         dκ = zero(κ)
         Enzyme.autodiff(
diff --git a/test/exchange_kernels.jl b/test/exchange_kernels.jl
index 11d1514..2f23a9f 100644
--- a/test/exchange_kernels.jl
+++ b/test/exchange_kernels.jl
@@ -114,10 +114,10 @@
             refined && refine!(bf, x -> x[1] < 0.5)
             sched = MFO._exchange_schedule(bf)
             ds = MFO._flatten_schedule(sched, bf, cpu)
-            xr = pack(set!(scalar_field(bf), fun))
+            xr = pack(set!(fun, scalar_field(bf)))
             MFO._run_exchange_host!(xr, sched)
             MFO._run_bc_host!(xr, bf, sched)
-            xk = run_device!(pack(set!(scalar_field(bf), fun)), bf, ds)
+            xk = run_device!(pack(set!(fun, scalar_field(bf))), bf, ds)
             @test noncorner_parity(xk, xr)
             # operator action reads no corners ⇒ exact equality through the sweep
             L = laplacian(bf)
@@ -128,10 +128,10 @@
                 1:MFO.nleaves(bf),
             )
             # SVector fields ride the same descriptors
-            wr = pack(set!(vector_field(bf), vfun))
+            wr = pack(set!(vfun, vector_field(bf)))
             MFO._run_exchange_host!(wr, sched)
             MFO._run_bc_host!(wr, bf, sched)
-            wk = run_device!(pack(set!(vector_field(bf), vfun)), bf, ds)
+            wk = run_device!(pack(set!(vfun, vector_field(bf))), bf, ds)
             @test noncorner_parity(wk, wr)
         end
     end
@@ -145,7 +145,7 @@
         refine!(bf, x -> x[2] < 0.5)
         dims = (bf.blocksize..., MFO.nleaves(bf))
         for mk in (scalar_field, vector_field)
-            f = pack(set!(mk(bf), mk === scalar_field ? fun : vfun))
+            f = pack(set!(mk === scalar_field ? fun : vfun, mk(bf)))
             # gather (interior_to_flat!): broadcast body vs loops, α/β combos
             vref = flatten(f)
             for (α, β) in ((1.0, 0.0), (2.0, 0.5))
@@ -199,7 +199,7 @@
             refine!(bfn, x -> x[1] < 0.5)
             sched = MFO._exchange_schedule(bfn)
             ds = MFO._flatten_schedule(sched, bfn, cpu)
-            x = pack(set!(scalar_field(bfn), fun))
+            x = pack(set!(fun, scalar_field(bfn)))
             a, s = alloc_device(x, bfn, ds)
             @test isfinite(s)
             (a, MFO.nleaves(bfn))
diff --git a/test/exchange_schedule.jl b/test/exchange_schedule.jl
index 0dd2be2..b127982 100644
--- a/test/exchange_schedule.jl
+++ b/test/exchange_schedule.jl
@@ -272,7 +272,7 @@ end
         interior_t = (h[1] + 1):(h[1] + n[1])             # transverse-interior band 2:5
         for bc in (periodic, mixed)
             bf = make_bf(bc)
-            f = set!(scalar_field(bf), x -> sinpi(x[1]) * cospi(2x[2]) + 0.5x[1])
+            f = set!(x -> sinpi(x[1]) * cospi(2x[2]) + 0.5x[1], scalar_field(bf))
             halo_update!(f, bf)
             for (i, K) in enumerate(bf.forest.leaves), d in 1:2, s in (-1, 1)
                 nbr = MFO.face_neighbor(bf.forest, K, d, s)
@@ -356,7 +356,7 @@ end
             lhs ≈ rhs || @info "adjoint identity" N T lhs rhs
             # constant reproduction (to roundoff of the weighted sum) on every
             # fill-owned ghost cell
-            c = set!(scalar_field(bf), _ -> T(0.75))
+            c = set!(_ -> T(0.75), scalar_field(bf))
             halo_update!(c, bf)
             tol = 50 * eps(T)
             for (phase, fills) in (("interp", sched.interp), ("restrict", sched.restrict))
@@ -656,14 +656,14 @@ end
         @test s2 !== s1
         @test s2.generation == bf.forest.generation[]
         @test length(s2.copies) == count_faces(bf)      # 4× the leaves ⇒ more descriptors
-        f = set!(scalar_field(bf), x -> x[1] - 2x[2])   # fresh field exchanges on the new forest
+        f = set!(x -> x[1] - 2x[2], scalar_field(bf))   # fresh field exchanges on the new forest
         @test halo_update!(f, bf) === f
     end
 
     @testset "inference and zero-allocation" begin
         bf = make_bf(((Dirichlet(), Dirichlet()), (Neumann(), Neumann())))
         @inferred MFO._exchange_schedule(bf)
-        uf = set!(scalar_field(bf), x -> sinpi(x[1]) * x[2])
+        uf = set!(x -> sinpi(x[1]) * x[2], scalar_field(bf))
         @inferred apply_bc!(uf, bf)
         @inferred MFO.fold_bc!(uf, bf)
         function alloc_halo(f, g)
diff --git a/test/fields.jl b/test/fields.jl
index c1d0dcc..9c26374 100644
--- a/test/fields.jl
+++ b/test/fields.jl
@@ -19,18 +19,18 @@
 
     @testset "set! at cell centers" begin
         g = CartesianGrid(((0.0, 1.0),), (4,))
-        u = set!(scalar_field(g), x -> 2 * x[1])
+        u = set!(x -> 2 * x[1], scalar_field(g))
         @test vec(collect(interior(u))) ≈ [0.25, 0.75, 1.25, 1.75]
         @test u.data[1] == 0.0 && u.data[end] == 0.0
 
         g32 = CartesianGrid(((0.0f0, 1.0f0),), (4,))
-        u32 = set!(scalar_field(g32), x -> x[1]^2)
+        u32 = set!(x -> x[1]^2, scalar_field(g32))
         @test eltype(u32) === Float32
     end
 
     @testset "component extraction" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (3, 3))
-        v = set!(vector_field(g), x -> SVector(x[1], -x[2]))
+        v = set!(x -> SVector(x[1], -x[2]), vector_field(g))
         vx = component(v, 1)
         vy = component(v, 2)
         @test eltype(vx) === Float64
@@ -38,14 +38,14 @@
         @test collect(interior(vy)) ≈ getindex.(collect(interior(v)), 2)
         @test_throws ArgumentError component(v, 3)
 
-        u = set!(scalar_field(g), x -> x[1])
+        u = set!(x -> x[1], scalar_field(g))
         @test collect(interior(component(u, 1))) == collect(interior(u))
         @test_throws ArgumentError component(u, 2)
     end
 
     @testset "flat round trip (scalar)" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (3, 4))
-        u = set!(scalar_field(g), x -> x[1] + 10 * x[2])
+        u = set!(x -> x[1] + 10 * x[2], scalar_field(g))
         flat = flatten(u)
         @test length(flat) == 12
         w = scalar_field(g)
@@ -56,7 +56,7 @@
 
     @testset "flat round trip (SVector)" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (3, 4))
-        v = set!(vector_field(g), x -> SVector(x[1], -x[2]))
+        v = set!(x -> SVector(x[1], -x[2]), vector_field(g))
         flat = flatten(v)
         @test length(flat) == 24
         @test eltype(flat) === Float64
@@ -67,14 +67,14 @@
 
     @testset "interior_to_flat! axpby fusion" begin
         g = CartesianGrid(((0.0, 1.0),), (5,))
-        u = set!(scalar_field(g), x -> x[1])
+        u = set!(x -> x[1], scalar_field(g))
         v = ones(5)
         interior_to_flat!(v, u, 2.0, 3.0)
         @test v ≈ 2 .* vec(collect(interior(u))) .+ 3
         interior_to_flat!(v, u)
         @test v ≈ vec(collect(interior(u)))
 
-        vf = set!(vector_field(CartesianGrid(((0.0, 1.0),), (3,))), x -> SVector(x[1]))
+        vf = set!(x -> SVector(x[1]), vector_field(CartesianGrid(((0.0, 1.0),), (3,))))
         fv = fill(0.5, 3)
         interior_to_flat!(fv, vf, 1.0, -1.0)
         @test fv ≈ getindex.(vec(collect(interior(vf))), 1) .- 0.5
@@ -82,7 +82,7 @@
 
     @testset "field-level BC and copy/similar" begin
         g = CartesianGrid(((0.0, 1.0),), (4,); bc=((Neumann(), Neumann()),))
-        u = set!(scalar_field(g), x -> x[1])
+        u = set!(x -> x[1], scalar_field(g))
         apply_bc!(u)
         @test u.data[1] == u.data[2]
         @test u.data[end] == u.data[end - 1]
diff --git a/test/forest_amr.jl b/test/forest_amr.jl
index 7920a13..4b5062e 100644
--- a/test/forest_amr.jl
+++ b/test/forest_amr.jl
@@ -44,7 +44,7 @@
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
             refine!(bf, predicate)
             @test !bf.forest.uniform[]
-            u = set!(scalar_field(bf), quad)
+            u = set!(quad, scalar_field(bf))
             halo_update!(u, bf)
             nchecked = 0
             foreach_cf_ghost(u, bf) do val, center
@@ -65,7 +65,7 @@
         refine!(bf, x -> x[1] < 0.2 && x[2] < 0.2)   # levels 0–2, balance! keeps 2:1
         @test minimum(k -> k.level, bf.forest.leaves) == 0
         @test maximum(k -> k.level, bf.forest.leaves) == 2
-        u = set!(scalar_field(bf), quad)
+        u = set!(quad, scalar_field(bf))
         halo_update!(u, bf)
         foreach_cf_ghost(u, bf) do val, center
             @test val ≈ quad(center) atol = 1e-12
@@ -87,7 +87,7 @@
             # resolutions): CF interface at x = π/2 and a wrapped one at x = 0
             refine!(bf, x -> x[1] < 1.6)
             @test !bf.forest.uniform[]
-            u = set!(scalar_field(bf), x -> sin(x[1]) * sin(x[2]))
+            u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(bf))
             Lu = laplacian(bf) * u
             e = 0.0
             for i in 1:MFO.nleaves(bf)
@@ -270,6 +270,6 @@
         @test_throws ArgumentError halo_update!(u2, bf2)
         # uniform forests keep the looser v1 constraints
         bf2u = BlockForest(base4; blocksize=(2, 2), maxlevel=2)
-        @test laplacian(bf2u) * set!(scalar_field(bf2u), x -> x[1]) isa BlockField
+        @test laplacian(bf2u) * set!(x -> x[1], scalar_field(bf2u)) isa BlockField
     end
 end
diff --git a/test/forest_diffusion.jl b/test/forest_diffusion.jl
index fef15a4..64b4d33 100644
--- a/test/forest_diffusion.jl
+++ b/test/forest_diffusion.jl
@@ -81,7 +81,7 @@
 
             @testset "$name: Diffusion, varying κ, $(nameof(typeof(avg)))" for avg in
                                                                                FOREST_DIFF_AVGS
-                κv = set!(scalar_field(bf), κ_varying)
+                κv = set!(κ_varying, scalar_field(bf))
                 y = diffusion(bf, κv; averaging=avg) * u
                 defect, scale = conservation_defect(y, bf)
                 @info "conservation Σ V·(Lu)" case = "$name varying-κ $(nameof(typeof(avg)))" defect scale
@@ -131,10 +131,10 @@
         for bc in PARITY_BCS, avg in FOREST_DIFF_AVGS
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)   # level 0 ⇒ same (8, 8)
-            Dg = diffusion(g, set!(scalar_field(g), κ_fun); averaging=avg)
-            Df = diffusion(bf, set!(scalar_field(bf), κ_fun); averaging=avg)
-            ref = collect(interior(Dg * set!(scalar_field(g), u_fun)))
-            rec = reconstruct(Df * set!(scalar_field(bf), u_fun), (8, 8))
+            Dg = diffusion(g, set!(κ_fun, scalar_field(g)); averaging=avg)
+            Df = diffusion(bf, set!(κ_fun, scalar_field(bf)); averaging=avg)
+            ref = collect(interior(Dg * set!(u_fun, scalar_field(g))))
+            rec = reconstruct(Df * set!(u_fun, scalar_field(bf)), (8, 8))
             @test rec == ref                                    # bit-identical
         end
     end
@@ -146,10 +146,10 @@
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
             refine!(bf, _ -> true)                              # uniform level 1 = 16×16
             @test all(k -> k.level == 1, bf.forest.leaves)
-            Dg = diffusion(g16, set!(scalar_field(g16), κ_fun))
-            Df = diffusion(bf, set!(scalar_field(bf), κ_fun))
-            ref = collect(interior(Dg * set!(scalar_field(g16), u_fun)))
-            rec = reconstruct(Df * set!(scalar_field(bf), u_fun), (16, 16))
+            Dg = diffusion(g16, set!(κ_fun, scalar_field(g16)))
+            Df = diffusion(bf, set!(κ_fun, scalar_field(bf)))
+            ref = collect(interior(Dg * set!(u_fun, scalar_field(g16))))
+            rec = reconstruct(Df * set!(u_fun, scalar_field(bf)), (16, 16))
             @test rec == ref
         end
     end
@@ -158,8 +158,8 @@
         bci = ((Dirichlet(2.0), Dirichlet(-1.0)), (Neumann(0.5), Dirichlet(3.0)))
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bci)
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        Dg = diffusion(g, set!(scalar_field(g), κ_fun))
-        Df = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        Dg = diffusion(g, set!(κ_fun, scalar_field(g)))
+        Df = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         bfor = boundary_rhs(Df, scalar_field(bf))
         @test bfor isa BlockField
         @test reconstruct(bfor, (8, 8)) == collect(interior(boundary_rhs(Dg, scalar_field(g))))
@@ -194,7 +194,7 @@
             avg in FOREST_DIFF_AVGS
 
             bf = small_refined(bc)
-            D = diffusion(bf, set!(scalar_field(bf), κ_fun); averaging=avg)
+            D = diffusion(bf, set!(κ_fun, scalar_field(bf)); averaging=avg)
             @test islinear(D) && isconstant(D) && !isdiagonal(D)
             @test !isselfadjoint(D)                  # CF coupling breaks the symmetry
             @test adjoint(D) isa AdjointOp           # must not fold to D
@@ -248,7 +248,7 @@
         bf = BlockForest(
             CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8)); blocksize=(4, 4), maxlevel=2
         )
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         @test isselfadjoint(D)
         @test adjoint(D) === D
         A = materialize(prepare(D))
@@ -263,7 +263,7 @@
         # skipping them.
         rng = Random.MersenneTwister(47)
         bf = small_refined(((Dirichlet(), Dirichlet()), (Neumann(), Neumann())))
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         ȳ = block_rand!(scalar_field(bf), bf, rng)
         x̄ = scalar_field(bf)
         for i in 1:MFO.nleaves(bf)
@@ -289,7 +289,7 @@
         base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bci)
         bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
         refine!(bf, x -> x[1] < 0.5 && x[2] < 0.5)
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         rng = Random.MersenneTwister(53)
         u = block_rand!(scalar_field(bf), bf, rng)
 
@@ -335,7 +335,7 @@
         # the reflected point.
         κlin(x) = 0.7 + 2.0 * x[1] - 1.3 * x[2]
         bf = small_refined(((Dirichlet(), Dirichlet()), (Neumann(), Neumann())))
-        κx = MFO.fill_coefficient_ghosts!(copy(set!(scalar_field(bf), κlin)), bf)
+        κx = MFO.fill_coefficient_ghosts!(copy(set!(κlin, scalar_field(bf))), bf)
         n = bf.blocksize
         lo = ntuple(k -> bf.extent[k][1], 2)
         counts = Dict("mirror" => 0, "copy" => 0, "inject" => 0, "average" => 0)
@@ -380,7 +380,7 @@
     @testset "the operator reads its exchanged κ ghosts (negative control)" begin
         rng = Random.MersenneTwister(59)
         bf = small_refined(((Periodic(), Periodic()), (Periodic(), Periodic())))
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         u = block_rand!(scalar_field(bf), bf, rng)
         clean = D * copy(u)
         κz = copy(D.κ)
@@ -400,7 +400,7 @@
 
     @testset "mutating κ after construction is inert" begin
         bf = small_refined(((Dirichlet(), Dirichlet()), (Dirichlet(), Dirichlet())))
-        κ = set!(scalar_field(bf), κ_fun)
+        κ = set!(κ_fun, scalar_field(bf))
         D = diffusion(bf, κ)
         rng = Random.MersenneTwister(61)
         u = block_rand!(scalar_field(bf), bf, rng)
@@ -439,11 +439,11 @@
         @test diffusion(bf, κpm; check=false) isa Diffusion
 
         # operator_diagonal stays unavailable on forests, matching Laplacian
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
         @test_throws ArgumentError operator_diagonal(D)
 
         # a regrid invalidates the operator through its κ generation stamp
-        κold = set!(scalar_field(bf), κ_fun)
+        κold = set!(κ_fun, scalar_field(bf))
         Dold = diffusion(bf, κold)
         nl0 = MFO.nleaves(bf)
         refine!(bf, x -> x[1] > 0.7 && x[2] > 0.7)   # matches the top-right leaf CENTER
@@ -478,7 +478,7 @@
             return a, sum(out)                          # DCE-proof: consume the output
         end
         function forest_alloc(bf, adj)
-            L = diffusion(bf, set!(scalar_field(bf), κ_fun))
+            L = diffusion(bf, set!(κ_fun, scalar_field(bf)))
             P = prepare(adj ? adjoint(L) : L)
             v = rand(Random.MersenneTwister(3), size(P, 2))
             a, s = alloc_mul(P, similar(v), v)
@@ -501,8 +501,8 @@
 
     @testset "packed coefficient path (prepared prototype packs κ)" begin
         bf = small_refined(((Dirichlet(), Dirichlet()), (Neumann(), Neumann())))
-        D = diffusion(bf, set!(scalar_field(bf), κ_fun))
-        uf = set!(scalar_field(bf), u_fun)
+        D = diffusion(bf, set!(κ_fun, scalar_field(bf)))
+        uf = set!(u_fun, scalar_field(bf))
         P = prepare(D, pack(uf))
         @test P.xpad isa PackedBlockField            # packed prototype ⇒ packed scratch
         @test P.op.κ isa PackedBlockField            # _prepare_tree packed the coefficient
@@ -531,8 +531,8 @@
             )
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
-            uf = set!(scalar_field(bf), u_fun)
-            D = diffusion(bf, set!(scalar_field(bf), κ_fun); averaging=avg)
+            uf = set!(u_fun, scalar_field(bf))
+            D = diffusion(bf, set!(κ_fun, scalar_field(bf)); averaging=avg)
             Dp = MFO.Diffusion(bf, pack(D.κ), D.avg)  # inner ctor: ghosts ride pack
 
             # public packed path: packed κ on packed x hits the new override; a
@@ -579,7 +579,7 @@
             # direct adjoint launch: full padded equality pre-fold (the ghost
             # cotangents are the point of the padded ndrange)
             ndp = (bf.blocksize .+ 2 .* bf.halo..., MFO.nleaves(bf))
-            ȳ = pack(set!(scalar_field(bf), x -> cospi(x[1]) + x[2]^2))
+            ȳ = pack(set!(x -> cospi(x[1]) + x[2]^2, scalar_field(bf)))
             ȳk = copy(ȳ)
             ȳr = copy(ȳ)
             x̄k = MFO.allocate_input(Dp, ȳk)
diff --git a/test/forest_packed.jl b/test/forest_packed.jl
index 1ca1927..2d29dec 100644
--- a/test/forest_packed.jl
+++ b/test/forest_packed.jl
@@ -19,13 +19,13 @@
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            uf = set!(scalar_field(bf), fun)
+            uf = set!(fun, scalar_field(bf))
             for makeL in (laplacian, g -> derivative(g, 1; order=1))
                 L = makeL(bf)
                 @test interiors_equal(L * pack(uf), L * uf)
             end
             refine!(bf, x -> x[1] < 0.5)
-            ur = set!(scalar_field(bf), fun)
+            ur = set!(fun, scalar_field(bf))
             @test interiors_equal(laplacian(bf) * pack(ur), laplacian(bf) * ur)
         end
     end
@@ -41,7 +41,7 @@
             )
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)   # kernel reads per-leaf levels
-            u = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(bf))
             x = pack(u)
             MFO.halo_update!(x, bf)
             MFO.apply_bc!(x, bf)
@@ -73,7 +73,7 @@
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            uf = set!(scalar_field(bf), fun)
+            uf = set!(fun, scalar_field(bf))
             # Uniform forest: isselfadjoint(Laplacian) is live-true, so the adjoint
             # action IS the kernel sweep.
             At = apply_adjoint!(similar(pack(uf)), laplacian(bf), pack(uf), bf)
@@ -113,14 +113,14 @@
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            uf = set!(scalar_field(bf), fun)
+            uf = set!(fun, scalar_field(bf))
             p = pack(uf)
             S = 2.0 * laplacian(bf) + adjoint(derivative(bf, 1; order=1))
             @test interiors_equal(S * copy(p), S * copy(uf))
             DG = divergence(bf) * MFO.gradient(bf)   # packed SVector intermediate
             @test interiors_equal(DG * copy(p), DG * copy(uf))
             @test interiors_equal(MFO.gradient(bf) * p, MFO.gradient(bf) * uf)
-            w = set!(vector_field(bf), vfun)
+            w = set!(vfun, vector_field(bf))
             @test interiors_equal(divergence(bf) * pack(w), divergence(bf) * w)
         end
     end
@@ -142,7 +142,7 @@
         bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        uf = set!(scalar_field(bf), fun)
+        uf = set!(fun, scalar_field(bf))
         v = flatten(uf)
         out = similar(v)
         A = prepare(laplacian(bf), pack(uf))
@@ -190,7 +190,7 @@
         for n in (16, 32)
             gn = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (n, n); bc=bc)
             bfn = BlockForest(gn; blocksize=(4, 4), maxlevel=2)
-            un = set!(scalar_field(bfn), fun)
+            un = set!(fun, scalar_field(bfn))
             vn = flatten(un)
             P = prepare(laplacian(bfn), pack(un))
             a, s = alloc_mul(P, similar(vn), vn)
@@ -216,7 +216,7 @@
         allocs = map((16, 32)) do n
             gn = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (n, n); bc=bc)
             bfn = BlockForest(gn; blocksize=(4, 4), maxlevel=2)
-            x = pack(set!(scalar_field(bfn), fun))
+            x = pack(set!(fun, scalar_field(bfn)))
             MFO.halo_update!(x, bfn)
             MFO.apply_bc!(x, bfn)
             a, s = alloc_launch(kernel!, MFO._zero_all!(similar(x)), x, bfn)
@@ -232,9 +232,9 @@
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
-            uf = set!(scalar_field(bf), fun)
-            κ = set!(scalar_field(bf), x -> 1 + x[1]^2 + 0.5 * x[2])
-            vel = set!(vector_field(bf), vfun)
+            uf = set!(fun, scalar_field(bf))
+            κ = set!(x -> 1 + x[1]^2 + 0.5 * x[2], scalar_field(bf))
+            vel = set!(vfun, vector_field(bf))
             for makeL in (
                 g -> derivative(g, 2; order=2),
                 MFO.gradient,
@@ -248,7 +248,7 @@
                 L = makeL(bf)
                 @test interiors_equal(L * pack(uf), L * uf)
             end
-            w = set!(vector_field(bf), vfun)
+            w = set!(vfun, vector_field(bf))
             @test interiors_equal(divergence(bf) * pack(w), divergence(bf) * w)
             A = advection(bf, SelfAdvection())
             @test interiors_equal(apply(A, pack(w)), apply(A, w))
@@ -265,14 +265,14 @@
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
             nd = (bf.blocksize..., MFO.nleaves(bf))
-            x = pack(set!(scalar_field(bf), fun))
+            x = pack(set!(fun, scalar_field(bf)))
             MFO.halo_update!(x, bf)
             MFO.apply_bc!(x, bf)
-            xw = pack(set!(vector_field(bf), vfun))
+            xw = pack(set!(vfun, vector_field(bf)))
             MFO.halo_update!(xw, bf)
             MFO.apply_bc!(xw, bf)
-            κp = pack(set!(scalar_field(bf), x -> 1 + x[1]^2))
-            velp = pack(set!(vector_field(bf), vfun))
+            κp = pack(set!(x -> 1 + x[1]^2, scalar_field(bf)))
+            velp = pack(set!(vfun, vector_field(bf)))
             cases = (
                 (
                     derivative(bf, 1; order=1), x,
@@ -360,8 +360,8 @@
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
             ndp = (bf.blocksize .+ 2 .* bf.halo..., MFO.nleaves(bf))
-            ȳs = pack(set!(scalar_field(bf), fun))     # scalar cotangent
-            ȳv = pack(set!(vector_field(bf), vfun))    # vector cotangent
+            ȳs = pack(set!(fun, scalar_field(bf)))     # scalar cotangent
+            ȳv = pack(set!(vfun, vector_field(bf)))    # vector cotangent
             cases = (
                 (
                     laplacian(bf), ȳs,
@@ -413,11 +413,11 @@
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
-            vel = set!(vector_field(bf), vfun)
-            xs = pack(set!(scalar_field(bf), fun))
-            ys = pack(set!(scalar_field(bf), gfun))
-            yv = pack(set!(vector_field(bf), x -> SVector(gfun(x), fun(x))))
-            xv = pack(set!(vector_field(bf), vfun))
+            vel = set!(vfun, vector_field(bf))
+            xs = pack(set!(fun, scalar_field(bf)))
+            ys = pack(set!(gfun, scalar_field(bf)))
+            yv = pack(set!(x -> SVector(gfun(x), fun(x)), vector_field(bf)))
+            xv = pack(set!(vfun, vector_field(bf)))
             for L in (derivative(bf, 1; order=1), derivative(bf, 2; order=2))
                 Lx = apply(L, copy(xs))
                 Lty = apply_adjoint!(similar(xs), L, copy(ys), bf)
@@ -448,9 +448,9 @@
             )
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
             refined && refine!(bf, x -> x[1] < 0.5)
-            κc = set!(scalar_field(bf, ComplexF64), x -> (1 + x[1]) + im * x[2])
-            xs = pack(set!(scalar_field(bf, ComplexF64), x -> fun(x) + 0.5im * x[1]))
-            ys = pack(set!(scalar_field(bf, ComplexF64), x -> gfun(x) - im * x[2]))
+            κc = set!(x -> (1 + x[1]) + im * x[2], scalar_field(bf, ComplexF64))
+            xs = pack(set!(x -> fun(x) + 0.5im * x[1], scalar_field(bf, ComplexF64)))
+            ys = pack(set!(x -> gfun(x) - im * x[2], scalar_field(bf, ComplexF64)))
             for S in (scaling(κc), scaling(pack(κc)), scaling(1.5 + 2.0im))
                 Sx = apply(S, copy(xs))
                 Sty = apply_adjoint!(similar(xs), S, copy(ys), bf)
@@ -463,9 +463,9 @@
         bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        uf = set!(scalar_field(bf), fun)
-        κ = set!(scalar_field(bf), x -> 1 + x[1]^2 + 0.5 * x[2])
-        vel = set!(vector_field(bf), vfun)
+        uf = set!(fun, scalar_field(bf))
+        κ = set!(x -> 1 + x[1]^2 + 0.5 * x[2], scalar_field(bf))
+        vel = set!(vfun, vector_field(bf))
         v = flatten(uf)
 
         K = divergence(bf) * scaling(κ) * MFO.gradient(bf)
@@ -511,18 +511,18 @@
         bc = ((Dirichlet(), Dirichlet()), (Dirichlet(), Dirichlet()))
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        κp = pack(set!(scalar_field(bf), x -> 1 + x[1]))
-        velp = pack(set!(vector_field(bf), vfun))
+        κp = pack(set!(x -> 1 + x[1], scalar_field(bf)))
+        velp = pack(set!(vfun, vector_field(bf)))
         S = scaling(κp)
         A = advection(bf, velp)
         refine!(bf, _ -> true)
-        xfresh = pack(set!(scalar_field(bf), fun))
+        xfresh = pack(set!(fun, scalar_field(bf)))
         @test_throws ArgumentError apply(S, xfresh)   # stale packed coefficient
         @test_throws ArgumentError apply(A, xfresh)   # stale packed velocity
 
         bf2 = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        κ2 = set!(scalar_field(bf2), x -> 1 + x[1])
-        u2 = set!(scalar_field(bf2), fun)
+        κ2 = set!(x -> 1 + x[1], scalar_field(bf2))
+        u2 = set!(fun, scalar_field(bf2))
         P = prepare(divergence(bf2) * scaling(κ2) * MFO.gradient(bf2), pack(u2))
         v2 = flatten(u2)
         refine!(bf2, _ -> true)
@@ -554,10 +554,10 @@
         allocs = map((16, 32)) do n
             gn = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (n, n); bc=bc)
             bfn = BlockForest(gn; blocksize=(4, 4), maxlevel=2)
-            x = pack(set!(scalar_field(bfn), fun))
+            x = pack(set!(fun, scalar_field(bfn)))
             MFO.halo_update!(x, bfn)
             MFO.apply_bc!(x, bfn)
-            velp = pack(set!(vector_field(bfn), vfun))
+            velp = pack(set!(vfun, vector_field(bfn)))
             a1, s1 = alloc_adv(adv!, MFO._zero_all!(similar(x)), x, velp, bfn)
             @test isfinite(s1)
             ȳ = copy(x)
diff --git a/test/forest_parity.jl b/test/forest_parity.jl
index 2d59c90..0c8d753 100644
--- a/test/forest_parity.jl
+++ b/test/forest_parity.jl
@@ -36,8 +36,8 @@ struct NoTraitOp85 <: AbstractOperator end
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)   # level 0 ⇒ same (8, 8)
-            u = set!(scalar_field(g), fun)
-            uf = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(g))
+            uf = set!(fun, scalar_field(bf))
             for makeL in leaf_ops
                 ref = collect(interior(makeL(g) * u))
                 rec = reconstruct(makeL(bf) * uf, (8, 8))
@@ -53,8 +53,8 @@ struct NoTraitOp85 <: AbstractOperator end
             bf = BlockForest(base; blocksize=(4, 4), maxlevel=3)
             refine!(bf, _ -> true)                              # uniform level 1 = 16×16
             @test all(k -> k.level == 1, bf.forest.leaves)
-            u = set!(scalar_field(g16), fun)
-            uf = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(g16))
+            uf = set!(fun, scalar_field(bf))
             for makeL in leaf_ops
                 ref = collect(interior(makeL(g16) * u))
                 rec = reconstruct(makeL(bf) * uf, (16, 16))
@@ -109,7 +109,7 @@ struct NoTraitOp85 <: AbstractOperator end
         @test size(A) == (64, 64)                   # nleaves·prod(blocksize) = 4·16
         M = materialize(A)
         @test M ≈ M'                                # forest Laplacian is self-adjoint
-        uf = set!(scalar_field(bf), fun)
+        uf = set!(fun, scalar_field(bf))
         v = flatten(uf)
         out = similar(v)
         mul!(out, A, v)
@@ -120,8 +120,8 @@ struct NoTraitOp85 <: AbstractOperator end
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            u = set!(scalar_field(g), fun)
-            uf = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(g))
+            uf = set!(fun, scalar_field(bf))
             D = derivative(g, 1; order=1)            # order 1 ⇒ adjoint is an AdjointOp
             Df = derivative(bf, 1; order=1)
             @test Df' isa AdjointOp
@@ -156,8 +156,8 @@ struct NoTraitOp85 <: AbstractOperator end
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            u = set!(scalar_field(g), fun)
-            uf = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(g))
+            uf = set!(fun, scalar_field(bf))
             L2g = laplacian(g) * laplacian(g)
             L2f = laplacian(bf) * laplacian(bf)
             @test reconstruct(L2f * uf, (8, 8)) == collect(interior(L2g * u))
@@ -170,7 +170,7 @@ struct NoTraitOp85 <: AbstractOperator end
         bf = BlockForest(
             CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8)); blocksize=(4, 4), maxlevel=2
         )
-        uf = set!(scalar_field(bf), fun)
+        uf = set!(fun, scalar_field(bf))
         L2 = laplacian(bf) * laplacian(bf)
         A = prepare(L2, uf)
         v = flatten(uf)
@@ -186,7 +186,7 @@ struct NoTraitOp85 <: AbstractOperator end
         # rank-changing composition: the intermediate is a vector BlockField whose
         # inter-block exchange must reproduce the single-grid wide Laplacian exactly
         g8 = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
-        u8 = set!(scalar_field(g8), fun)
+        u8 = set!(fun, scalar_field(g8))
         wide = apply(divergence(g8), apply(MFO.gradient(g8), u8))
         DG = divergence(bf) * MFO.gradient(bf)
         @test reconstruct(DG * copy(uf), (8, 8)) == collect(interior(wide))
@@ -197,12 +197,12 @@ struct NoTraitOp85 <: AbstractOperator end
         for bc in bcs
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8); bc=bc)
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            u = set!(scalar_field(g), fun)
-            uf = set!(scalar_field(bf), fun)
+            u = set!(fun, scalar_field(g))
+            uf = set!(fun, scalar_field(bf))
             @test reconstruct(MFO.gradient(bf) * uf, (8, 8)) ==
                 collect(interior(MFO.gradient(g) * u))
-            w = set!(vector_field(g), vfun)
-            wf = set!(vector_field(bf), vfun)
+            w = set!(vfun, vector_field(g))
+            wf = set!(vfun, vector_field(bf))
             @test reconstruct(divergence(bf) * wf, (8, 8)) ==
                 collect(interior(divergence(g) * w))
         end
@@ -261,8 +261,8 @@ struct NoTraitOp85 <: AbstractOperator end
         @testset "trait table" begin
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
             bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-            κ = set!(scalar_field(bf), x -> 1 + x[1])
-            vel = set!(vector_field(bf), x -> SVector(1.0, 0.5))
+            κ = set!(x -> 1 + x[1], scalar_field(bf))
+            vel = set!(x -> SVector(1.0, 0.5), vector_field(bf))
             for L in (
                 laplacian(bf), derivative(bf, 1; order=1), derivative(bf, 2; order=2),
                 MFO.gradient(bf), divergence(bf), scaling(κ), scaling(2.0), identity_op(),
@@ -317,7 +317,7 @@ struct NoTraitOp85 <: AbstractOperator end
                 balance!(bf)
                 @test !bf.forest.uniform[]
             end
-            u = set!(scalar_field(bf), f)
+            u = set!(f, scalar_field(bf))
             T = eltype(u)
             @test T === eltype(groot.spacing)
             N = length(bs)
@@ -352,7 +352,7 @@ struct NoTraitOp85 <: AbstractOperator end
                 @test count_exchanges!(ys, T(3) * aniso, x, bf) == (1, 1)
                 @test interiors_equal(ys, per_operand!(similar(x), T(3) * aniso, x, bf, true, false))
                 # a diagonal operand rides along on the shared exchange
-                κ = set!(scalar_field(bf), z -> 1 + z[1] * z[2])
+                κ = set!(z -> 1 + z[1] * z[2], scalar_field(bf))
                 κx = x isa PackedBlockField ? pack(κ) : κ
                 mixed = laplacian(bf) + scaling(κx)
                 ym = similar(x)
@@ -384,8 +384,8 @@ struct NoTraitOp85 <: AbstractOperator end
             κ_fun = z -> 1 + z[1] + 0.5 * z[2]^2      # variable κ: rewrite ≠ restriction
             # uniform: cfflux is empty, so the sum shares one exchange, bit-exactly
             bfu = BlockForest(g; blocksize=(4, 4), maxlevel=3)
-            u = set!(scalar_field(bfu), fun)
-            Du = diffusion(bfu, set!(scalar_field(bfu), κ_fun)) + laplacian(bfu)
+            u = set!(fun, scalar_field(bfu))
+            Du = diffusion(bfu, set!(κ_fun, scalar_field(bfu))) + laplacian(bfu)
             @test shares_exchange(Du)
             yu = similar(u)
             @test count_exchanges!(yu, Du, u, bfu) == (1, 1)
@@ -394,8 +394,8 @@ struct NoTraitOp85 <: AbstractOperator end
             bfr = BlockForest(g; blocksize=(4, 4), maxlevel=3)
             refine!(bfr, x -> x[1] < 0.5 && x[2] < 0.5)
             balance!(bfr)
-            ur = set!(scalar_field(bfr), fun)
-            Dr = diffusion(bfr, set!(scalar_field(bfr), κ_fun))
+            ur = set!(fun, scalar_field(bfr))
+            Dr = diffusion(bfr, set!(κ_fun, scalar_field(bfr)))
             Lr = Dr + laplacian(bfr)
             @test !shares_exchange(Lr)
             yr = similar(ur)
diff --git a/test/forest_prepare.jl b/test/forest_prepare.jl
index 544ffc9..9d4070e 100644
--- a/test/forest_prepare.jl
+++ b/test/forest_prepare.jl
@@ -8,7 +8,7 @@
     bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
     g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
     bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)   # 16 leaves
-    uf = set!(scalar_field(bf), fun)
+    uf = set!(fun, scalar_field(bf))
     v = flatten(uf)
     out = similar(v)
 
@@ -68,7 +68,7 @@
         @test out ≈ materialize(prepare(inner))' * v
         @test alloc_mul(At, out, v) ≤ alloc_bound(MFO.nleaves(bf))
 
-        κ = set!(scalar_field(bf), x -> 1 + x[1] * x[2])
+        κ = set!(x -> 1 + x[1] * x[2], scalar_field(bf))
         K = divergence(bf) * scaling(κ)                          # vector → scalar
         Kt = prepare(MFO.AdjointOp(K), scalar_field(bf))
         B = materialize(prepare(K, vector_field(bf)))
@@ -122,7 +122,7 @@
         mul!(out2, A, v, 2.0, 3.0)
         @test out2 ≈ 2.0 .* per_operand_flat(aniso, uf) .+ 3.0 .* v
         # nested sums and a diagonal operand share too
-        κ = set!(scalar_field(bf), x -> 1 + x[1] * x[2])
+        κ = set!(x -> 1 + x[1] * x[2], scalar_field(bf))
         three = (aniso + scaling(κ)) + identity_op()
         A3 = prepare(three)
         @test prepared_exchanges(A3, v) == (1, 1)
@@ -162,7 +162,7 @@
         bfr = BlockForest(g; blocksize=(4, 4), maxlevel=3)
         refine!(bfr, x -> x[1] < 0.5 && x[2] < 0.5)
         balance!(bfr)
-        ur = set!(scalar_field(bfr), fun)
+        ur = set!(fun, scalar_field(bfr))
         vr = flatten(ur)
         anisor = 0.13 * laplacian(bfr) + 0.7 * derivative(bfr, 1; order=2)
         @test !isselfadjoint(anisor)
diff --git a/test/gradient_divergence.jl b/test/gradient_divergence.jl
index fb7326f..144ccf6 100644
--- a/test/gradient_divergence.jl
+++ b/test/gradient_divergence.jl
@@ -8,28 +8,28 @@ end
 @testset "Gradient and Divergence (rank-changers)" begin
     @testset "gradient analytic action" begin
         g = periodic_grid_2d(32)
-        u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
+        u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
         ∇u = MatrixFreeOperators.gradient(g) * u
         @test eltype(∇u) === SVector{2,Float64}
-        ref = set!(vector_field(g), x -> SVector(cos(x[1]) * sin(x[2]), sin(x[1]) * cos(x[2])))
+        ref = set!(x -> SVector(cos(x[1]) * sin(x[2]), sin(x[1]) * cos(x[2])), vector_field(g))
         @test maximum(norm.(collect(interior(∇u)) .- collect(interior(ref)))) < 0.01
     end
 
     @testset "divergence analytic action" begin
         g = periodic_grid_2d(32)
-        v = set!(vector_field(g), x -> SVector(sin(x[1]) * cos(x[2]), cos(x[1]) * sin(x[2])))
+        v = set!(x -> SVector(sin(x[1]) * cos(x[2]), cos(x[1]) * sin(x[2])), vector_field(g))
         divv = divergence(g) * v
         @test eltype(divv) === Float64
-        ref = set!(scalar_field(g), x -> 2 * cos(x[1]) * cos(x[2]))
+        ref = set!(x -> 2 * cos(x[1]) * cos(x[2]), scalar_field(g))
         @test maximum(abs, collect(interior(divv)) .- collect(interior(ref))) < 0.02
     end
 
     @testset "div ∘ grad agrees with laplacian (both vs analytic)" begin
         g = periodic_grid_2d(48)
-        u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
+        u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
         wide = apply(divergence(g), apply(MatrixFreeOperators.gradient(g), u))
         compact = laplacian(g) * u
-        ref = -2 .* collect(interior(set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))))
+        ref = -2 .* collect(interior(set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))))
         @test maximum(abs, collect(interior(wide)) .- ref) < 0.02
         @test maximum(abs, collect(interior(compact)) .- ref) < 0.02
     end
diff --git a/test/laplacian.jl b/test/laplacian.jl
index 26c63a9..80af7fc 100644
--- a/test/laplacian.jl
+++ b/test/laplacian.jl
@@ -4,7 +4,7 @@ function lap_periodic_error(n::Int, ::Val{D}) where {D}
         ntuple(_ -> n, Val(D));
         bc=ntuple(_ -> (Periodic(), Periodic()), Val(D)),
     )
-    u = set!(scalar_field(g), x -> prod(sin, x))
+    u = set!(x -> prod(sin, x), scalar_field(g))
     y = laplacian(g) * u
     return maximum(abs, collect(interior(y)) .+ D .* collect(interior(u)))
 end
@@ -62,9 +62,9 @@ end
             ((0.0, 2π), (0.0, 2π)), (16, 16);
             bc=((Periodic(), Periodic()), (Periodic(), Periodic())),
         )
-        v = set!(vector_field(g), x -> SVector(sin(x[1]) * sin(x[2]), cos(x[1])))
-        u1 = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
-        u2 = set!(scalar_field(g), x -> cos(x[1]))
+        v = set!(x -> SVector(sin(x[1]) * sin(x[2]), cos(x[1])), vector_field(g))
+        u1 = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
+        u2 = set!(x -> cos(x[1]), scalar_field(g))
         L = laplacian(g)
         Lv = L * v
         @test getindex.(collect(interior(Lv)), 1) ≈ collect(interior(L * u1))
@@ -73,7 +73,7 @@ end
 
     @testset "stencil primitive returns center value and Laplacian" begin
         g = CartesianGrid(((0.0, 1.0),), (4,); bc=((Neumann(), Neumann()),))
-        u = set!(scalar_field(g), x -> x[1]^2)
+        u = set!(x -> x[1]^2, scalar_field(g))
         apply_bc!(u)
         inv_h2 = inv.(spacing(g) .^ 2)
         uc, lap = laplacian_stencil(u.data, CartesianIndex(3), inv_h2)
diff --git a/test/linearize.jl b/test/linearize.jl
index f0c0db1..6cd85e5 100644
--- a/test/linearize.jl
+++ b/test/linearize.jl
@@ -5,7 +5,7 @@
             bc=((Dirichlet(), Dirichlet()), (Neumann(), Neumann())),
         )
         rng = Random.MersenneTwister(41)
-        for F in (laplacian(g), advection(g, set!(vector_field(g), x -> SVector(x[1], 1.0))))
+        for F in (laplacian(g), advection(g, set!(x -> SVector(x[1], 1.0), vector_field(g))))
             u0 = scalar_field(g)
             v = scalar_field(g)
             interior(u0) .= rand(rng, local_size(g)...)
@@ -21,15 +21,15 @@
     @testset "self-advection JVP matches directional finite difference" begin
         g = CartesianGrid(((0.0, 2π),), (32,); bc=((Periodic(), Periodic()),))
         F = advection(g, SelfAdvection())
-        u0 = set!(vector_field(g), x -> SVector(sin(x[1])))
-        v = set!(vector_field(g), x -> SVector(cos(2 * x[1])))
+        u0 = set!(x -> SVector(sin(x[1])), vector_field(g))
+        v = set!(x -> SVector(cos(2 * x[1])), vector_field(g))
         J = linearize(F, u0)
         Jv = apply(J, copy(v))
 
         # analytic: ∂/∂ε (u+εv)·∇(u+εv)|₀ = v·∇u + u·∇v
         ref = set!(
-            vector_field(g),
             x -> SVector(cos(2 * x[1]) * cos(x[1]) - 2 * sin(x[1]) * sin(2 * x[1])),
+            vector_field(g),
         )
         @test maximum(norm.(collect(interior(Jv)) .- collect(interior(ref)))) < 0.05
 
@@ -45,9 +45,9 @@
     @testset "linearize! refresh equals fresh linearize" begin
         g = CartesianGrid(((0.0, 2π),), (16,); bc=((Periodic(), Periodic()),))
         F = advection(g, SelfAdvection())
-        u0 = set!(vector_field(g), x -> SVector(sin(x[1])))
-        u1 = set!(vector_field(g), x -> SVector(cos(x[1])))
-        v = set!(vector_field(g), x -> SVector(sin(2 * x[1])))
+        u0 = set!(x -> SVector(sin(x[1])), vector_field(g))
+        u1 = set!(x -> SVector(cos(x[1])), vector_field(g))
+        v = set!(x -> SVector(sin(2 * x[1])), vector_field(g))
 
         J = linearize(F, u0)
         linearize!(J, u1)
@@ -59,12 +59,12 @@
     @testset "prepared Jacobian drives Krylov (implicit-Euler JFNK system)" begin
         g = CartesianGrid(((0.0, 2π),), (24,); bc=((Periodic(), Periodic()),))
         F = advection(g, SelfAdvection())
-        u0 = set!(vector_field(g), x -> SVector(2 + sin(x[1]) / 4))
+        u0 = set!(x -> SVector(2 + sin(x[1]) / 4), vector_field(g))
         J = linearize(F, u0)
         @test_throws ArgumentError adjoint(J)
         @test size(J) == (24, 24)
 
-        v = set!(vector_field(g), x -> SVector(cos(x[1])))
+        v = set!(x -> SVector(cos(x[1])), vector_field(g))
         Pj = prepare(J, u0)
         jv = similar(flatten(v))
         mul!(jv, Pj, flatten(v))
@@ -88,7 +88,7 @@
             bc=((Dirichlet(), Dirichlet()), (Neumann(), Neumann())),
         )
         rng = Random.MersenneTwister(93)
-        vel = set!(vector_field(g), x -> SVector(x[1], 1.0))
+        vel = set!(x -> SVector(x[1], 1.0), vector_field(g))
 
         @testset "Jacobian of the linear $(name) is the operator, exactly" for (name, F) in
                                                                               (
@@ -116,7 +116,7 @@
         @testset "nonlinear self-advection: JVP vs FD, and the adjoint identity" begin
             g1 = CartesianGrid(((0.0, 2π),), (16,); bc=((Periodic(), Periodic()),))
             F = advection(g1, SelfAdvection())
-            u0 = set!(vector_field(g1), x -> SVector(2 + sin(x[1])))
+            u0 = set!(x -> SVector(2 + sin(x[1])), vector_field(g1))
             J = linearize(F, u0, EnzymeJVP())
 
             x = vector_field(g1)
@@ -149,9 +149,9 @@
         @testset "linearize! refresh, and the unloaded-backend error" begin
             g1 = CartesianGrid(((0.0, 2π),), (16,); bc=((Periodic(), Periodic()),))
             F = advection(g1, SelfAdvection())
-            u0 = set!(vector_field(g1), x -> SVector(sin(x[1])))
-            u1 = set!(vector_field(g1), x -> SVector(cos(x[1])))
-            v = set!(vector_field(g1), x -> SVector(sin(2 * x[1])))
+            u0 = set!(x -> SVector(sin(x[1])), vector_field(g1))
+            u1 = set!(x -> SVector(cos(x[1])), vector_field(g1))
+            v = set!(x -> SVector(sin(2 * x[1])), vector_field(g1))
             J = linearize(F, u0, EnzymeJVP())
             linearize!(J, u1)
             fresh = linearize(F, u1, EnzymeJVP())
diff --git a/test/mdla_gpu.jl b/test/mdla_gpu.jl
index 1597f19..5de56d6 100644
--- a/test/mdla_gpu.jl
+++ b/test/mdla_gpu.jl
@@ -26,7 +26,7 @@ end
 @testset "guards" begin
     g = mdla_grid((Dirichlet(), Dirichlet()))
     gc = coarsen(g)
-    v = set!(vector_field(g), x -> SVector(1.0, 0.0))
+    v = set!(x -> SVector(1.0, 0.0), vector_field(g))
     @test_throws ArgumentError prepare_distributed(MatrixFreeOperators.gradient(g), 1)
     @test_throws ArgumentError prepare_distributed(MatrixFreeOperators.divergence(g), 1)
     @test_throws ArgumentError prepare_distributed(advection(g, v) + laplacian(g), 1)
@@ -35,13 +35,13 @@ end
     @test_throws ArgumentError prepare_distributed(laplacian(gc) * restriction(g, gc), 1)
     @test_throws ArgumentError prepare_distributed(laplacian(g), NGPUS_MDLA + 1)
     # a coefficient on some OTHER grid: each leaf is fine alone, only the tree shows it
-    κc = set!(scalar_field(gc), x -> 1 + x[1] / 7)
+    κc = set!(x -> 1 + x[1] / 7, scalar_field(gc))
     @test_throws ArgumentError prepare_distributed(laplacian(g) + scaling(κc), 1)
 end
 
 @testset "newly distributable operators are accepted" begin
     g = mdla_grid((Dirichlet(), Dirichlet()))
-    κ = set!(scalar_field(g), x -> 1 + x[1] / 7)
+    κ = set!(x -> 1 + x[1] / 7, scalar_field(g))
     for L in (
         laplacian(g) * laplacian(g),
         adjoint(derivative(g, 1)),
@@ -68,7 +68,7 @@ end
 # Cases with a *false* trait catch a forwarding that hard-codes `true`.
 @testset "traits forward to the global tree" begin
     g = mdla_grid((Dirichlet(), Neumann()))
-    κ = set!(scalar_field(g), x -> 1 + x[1] / 7)
+    κ = set!(x -> 1 + x[1] / 7, scalar_field(g))
     D1 = derivative(g, 1; order=1)
     for L in (
         laplacian(g),                                  # self-adjoint, not diagonal
@@ -351,7 +351,7 @@ end
         bc=((Dirichlet(), Dirichlet()), (Dirichlet(), Dirichlet())),
     )
     L = adjoint(derivative(g, 2)) * derivative(g, 2) + 1.0 * identity_op()
-    f = set!(scalar_field(g), x -> sin(π * x[1]) * sin(π * x[2]))
+    f = set!(x -> sin(π * x[1]) * sin(π * x[2]), scalar_field(g))
     bflat = flatten(f)
 
     u_cpu, stats_cpu = Krylov.cg(
@@ -376,7 +376,7 @@ end
 @testset "Krylov.cg on distributed Poisson" begin
     g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (24, 26))
     L = -1.0 * laplacian(g)   # SPD under homogeneous Dirichlet
-    f = set!(scalar_field(g), x -> sin(π * x[1]) * sin(π * x[2]))
+    f = set!(x -> sin(π * x[1]) * sin(π * x[2]), scalar_field(g))
     bflat = flatten(f)
     n = length(bflat)
 
@@ -411,7 +411,7 @@ mdla_coeff(x) = 1.5 + x[2] + 0.3 * x[1] * x[2] + 0.2 * x[2]^2
 @testset "the coefficient is sliced and uploaded per partition" begin
     if NGPUS_MDLA >= 2
         g = mdla_grid((Dirichlet(), Neumann()))
-        κ = set!(scalar_field(g), mdla_coeff)
+        κ = set!(mdla_coeff, scalar_field(g))
         P = prepare_distributed(laplacian(g) * scaling(κ), 2)
         locals = partition_grid(g, 2)
         for d in 1:2
@@ -435,7 +435,7 @@ end
 @testset "the diffusion coefficient is uploaded with its cut-plane ghosts" begin
     if NGPUS_MDLA >= 2
         g = mdla_grid((Dirichlet(), Neumann()))
-        κ = set!(scalar_field(g), mdla_coeff)
+        κ = set!(mdla_coeff, scalar_field(g))
         Dg = diffusion(g, κ)
         P = prepare_distributed(laplacian(g) * Dg, 2)
         locals = partition_grid(g, 2)
@@ -465,7 +465,7 @@ end
             g = mdla_grid(cutbc)
             n = prod(local_size(g))
             xflat = rand(rng, n)
-            κ = set!(scalar_field(g), mdla_coeff)
+            κ = set!(mdla_coeff, scalar_field(g))
             for L in (
                 scaling(κ),
                 laplacian(g) * scaling(κ),
@@ -499,7 +499,7 @@ end
         rng = Random.MersenneTwister(67)
         g = mdla_grid((Dirichlet(), Neumann()))
         n = prod(local_size(g))
-        κ = set!(scalar_field(g), mdla_coeff)
+        κ = set!(mdla_coeff, scalar_field(g))
         x, y = rand(rng, n), rand(rng, n)
         for L in (
             scaling(κ) * laplacian(g),
@@ -523,7 +523,7 @@ end
         gs = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (3, 6);
             bc=((Dirichlet(), Dirichlet()), (Dirichlet(), Neumann())))
         ns = prod(local_size(gs))
-        κs = set!(scalar_field(gs), mdla_coeff)
+        κs = set!(mdla_coeff, scalar_field(gs))
         for Lt in (
             scaling(κs) * laplacian(gs),
             adjoint(derivative(gs, 1) * scaling(κs)),
@@ -550,7 +550,7 @@ end
             ((0.0, 2π), (0.0, 1.0)), (16, 18);
             bc=((Dirichlet(0.75), Neumann(-1.25)), cut),
         )
-        κ = set!(scalar_field(g), mdla_coeff)
+        κ = set!(mdla_coeff, scalar_field(g))
         D1 = derivative(g, 1)
         for L in (
             laplacian(g),
@@ -585,7 +585,7 @@ end
         ((0.0, 1.0), (0.0, 1.0)), (24, 26);
         bc=((Dirichlet(0.5), Dirichlet(-0.25)), (Dirichlet(1.0), Dirichlet(-0.5))),
     )
-    κ = set!(scalar_field(g), mdla_coeff)
+    κ = set!(mdla_coeff, scalar_field(g))
     fun = x -> sin(π * x[1]) * sin(π * x[2]) + 0.3x[2]
     # `-diffusion` is the compact-form counterpart of `-(scaling(κ)*laplacian)`:
     # SPD, so cg applies, and exactly symmetric rather than merely close — an
@@ -597,7 +597,7 @@ end
         -1.0 * diffusion(g, κ),
         -1.0 * diffusion(g, κ; averaging=HarmonicMean()),
     )
-        bflat = flatten(set!(scalar_field(g), fun)) .- flatten(boundary_rhs(L, g))
+        bflat = flatten(set!(fun, scalar_field(g))) .- flatten(boundary_rhs(L, g))
         u_cpu, stats_cpu = Krylov.cg(prepare(L), bflat; atol=1e-10, rtol=1e-10)
         @test stats_cpu.solved
 
@@ -628,11 +628,11 @@ end
         )
         L = laplacian(g)
         fun = x -> sin(3x[1]) * exp(-x[2]) + 0.25x[1] * x[2]
-        ref = flatten(set!(scalar_field(g), fun)) .- flatten(boundary_rhs(L, g))
+        ref = flatten(set!(fun, scalar_field(g))) .- flatten(boundary_rhs(L, g))
         P = prepare_distributed(L, 2)
         grids = local_grids(P)
         @test length(grids) == 2
-        fields = [set!(scalar_field(lg), fun) for lg in grids]
+        fields = [set!(fun, scalar_field(lg)) for lg in grids]
         @test gather(assemble_rhs(P, fields)) == ref
         @test_throws ArgumentError assemble_rhs(P, fields[1:1])
     else
diff --git a/test/multigpu/mdla_3partition.jl b/test/multigpu/mdla_3partition.jl
index e6cbb12..bf87550 100644
--- a/test/multigpu/mdla_3partition.jl
+++ b/test/multigpu/mdla_3partition.jl
@@ -138,7 +138,7 @@ mp_grid(cutbc) = CartesianGrid(
         @testset "3-partition Krylov.cg parity" begin
             g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (24, 26))
             L = -1.0 * laplacian(g)   # SPD under homogeneous Dirichlet
-            f = set!(scalar_field(g), x -> sin(π * x[1]) * sin(π * x[2]))
+            f = set!(x -> sin(π * x[1]) * sin(π * x[2]), scalar_field(g))
             bflat = flatten(f)
             n = length(bflat)
 
@@ -164,7 +164,7 @@ mp_grid(cutbc) = CartesianGrid(
                 g = mp_grid(cutbc)
                 n = prod(local_size(g))
                 xflat = rand(rng, n)
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 for L in (
                     scaling(κ),
                     laplacian(g) * scaling(κ),
@@ -188,7 +188,7 @@ mp_grid(cutbc) = CartesianGrid(
             for cutbc in ((Dirichlet(), Neumann()), (Periodic(), Periodic()))
                 g = mp_grid(cutbc)
                 n = prod(local_size(g))
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 for L in (scaling(κ) * laplacian(g), laplacian(g) * scaling(κ))
                     P = prepare_distributed(L, 3)
                     x = MultiDeviceVector(rand(rng, n), P.spec)
@@ -211,7 +211,7 @@ mp_grid(cutbc) = CartesianGrid(
                     ((0.0, 2π), (0.0, 1.0)), (16, 18);
                     bc=((Dirichlet(0.75), Neumann(-1.25)), cut),
                 )
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 D1 = derivative(g, 1)
                 for L in (
                     laplacian(g),
@@ -233,7 +233,7 @@ mp_grid(cutbc) = CartesianGrid(
             )
             L = -1.0 * laplacian(g)
             fun = x -> sin(3x[1]) * exp(-x[2]) + 0.25x[1] * x[2]
-            bflat = flatten(set!(scalar_field(g), fun)) .- flatten(boundary_rhs(L, g))
+            bflat = flatten(set!(fun, scalar_field(g))) .- flatten(boundary_rhs(L, g))
             u_cpu, stats_cpu = Krylov.cg(prepare(L), bflat; atol=1e-10, rtol=1e-10)
             @test stats_cpu.solved
 
@@ -258,7 +258,7 @@ mp_grid(cutbc) = CartesianGrid(
         @testset "3-partition diffusion coefficient upload" begin
             for cutbc in ((Dirichlet(), Neumann()), (Periodic(), Periodic()))
                 g = mp_grid(cutbc)
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 Dg = diffusion(g, κ)
                 P = prepare_distributed(laplacian(g) * Dg, 3)
                 locals = partition_grid(g, 3)
@@ -282,7 +282,7 @@ mp_grid(cutbc) = CartesianGrid(
                 g = mp_grid(cutbc)
                 n = prod(local_size(g))
                 xflat = rand(rng, n)
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 for L in (
                     diffusion(g, κ),
                     diffusion(g, κ; averaging=HarmonicMean()),
@@ -306,7 +306,7 @@ mp_grid(cutbc) = CartesianGrid(
             for cutbc in ((Dirichlet(), Neumann()), (Periodic(), Periodic()))
                 g = mp_grid(cutbc)
                 n = prod(local_size(g))
-                κ = set!(scalar_field(g), mp_coeff)
+                κ = set!(mp_coeff, scalar_field(g))
                 for L in (diffusion(g, κ), laplacian(g) * diffusion(g, κ))
                     P = prepare_distributed(L, 3)
                     x = MultiDeviceVector(rand(rng, n), P.spec)
diff --git a/test/multigrid.jl b/test/multigrid.jl
index 14b77ea..9737aff 100644
--- a/test/multigrid.jl
+++ b/test/multigrid.jl
@@ -70,7 +70,7 @@
 
         # ScalingOp / IdentityOp leaves
         @test operator_diagonal(scaling(2.5)) == 2.5
-        κ = set!(scalar_field(g), x -> 1 + x[1]^2)
+        κ = set!(x -> 1 + x[1]^2, scalar_field(g))
         @test operator_diagonal(scaling(κ)) === κ
         @test operator_diagonal(identity_op()) === true
 
@@ -79,7 +79,7 @@
             ((0.0, 1.0),), (6,); bc=((Periodic(), Periodic()),)
         )
         @test operator_diagonal(3 * laplacian(gp)) ≈ 3 * operator_diagonal(laplacian(gp))
-        σ = set!(scalar_field(g), x -> x[1] + x[2])
+        σ = set!(x -> x[1] + x[2], scalar_field(g))
         M = scaling(σ) - laplacian(g)
         @test flatten(operator_diagonal(M)) ≈ diag(materialize(prepare(M)))
         @test operator_diagonal(scaling(2.0) * scaling(κ)) isa Field
@@ -97,7 +97,7 @@
         gf = CartesianGrid(((0.0, 1.0),), (8,))
         gc = coarsen(gf)
         P = prolongation(gc, gf)
-        u = set!(scalar_field(gc), x -> 2 * x[1] + 1)
+        u = set!(x -> 2 * x[1] + 1, scalar_field(gc))
         y = P * u
         uc = collect(interior(u))
         yf = collect(interior(y))
@@ -110,7 +110,7 @@
 
         # Neumann wall child mirrors: (3/4)u₁ + (1/4)u₁ = u₁; constants preserved
         gfn = CartesianGrid(((0.0, 1.0),), (8,); bc=((Neumann(), Neumann()),))
-        un = set!(scalar_field(coarsen(gfn)), x -> 3.5 + 0 * x[1])
+        un = set!(x -> 3.5 + 0 * x[1], scalar_field(coarsen(gfn)))
         @test all(collect(interior(prolongation(coarsen(gfn), gfn) * un)) .≈ 3.5)
 
         # P·1 = 1 under Periodic and Neumann, D ∈ 1:3
@@ -121,7 +121,7 @@
                 bc=ntuple(_ -> (bc, bc), D),
             )
             gcD = coarsen(gD)
-            ones_c = set!(scalar_field(gcD), _ -> 1.0)
+            ones_c = set!(_ -> 1.0, scalar_field(gcD))
             @test all(collect(interior(prolongation(gcD, gD) * ones_c)) .≈ 1.0)
         end
     end
@@ -147,7 +147,7 @@
                 ntuple(_ -> 8, D);
                 bc=ntuple(_ -> (bc, bc), D),
             )
-            ones_f = set!(scalar_field(gD), _ -> 1.0)
+            ones_f = set!(_ -> 1.0, scalar_field(gD))
             @test all(collect(interior(restriction(gD) * ones_f)) .≈ 1.0)
         end
     end
@@ -309,18 +309,18 @@
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
         L = -1 * laplacian(g)
         A = prepare(L)
-        b = flatten(set!(scalar_field(g), x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2])))
+        b = flatten(set!(x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
         u, _ = Krylov.cg(A, b; M=MultigridPreconditioner(L), rtol=1e-10)
-        uex = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
+        uex = flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
         @test maximum(abs, u .- uex) < 1e-3
     end
 
     @testset "variable-coefficient SPD system" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16))
-        σ = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+        σ = set!(x -> 1 + x[1] * x[2], scalar_field(g))
         L = scaling(σ) - laplacian(g)
         A = prepare(L)
-        f = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
+        f = flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
         u, stats = Krylov.cg(A, f; M=MultigridPreconditioner(L; levels=2), rtol=1e-10)
         r = similar(f)
         mul!(r, A, u)
@@ -333,13 +333,13 @@
     @testset "standalone MultigridSolver" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
         L = -1 * laplacian(g)
-        b = flatten(set!(scalar_field(g), x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2])))
+        b = flatten(set!(x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
         u = solve(MultigridSolver(L), b; rtol=1e-10)
         A = prepare(L)
         r = similar(b)
         mul!(r, A, u)
         @test norm(b .- r) <= 1e-10 * norm(b)
-        uex = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
+        uex = flatten(set!(x -> sinpi(x[1]) * sinpi(x[2]), scalar_field(g)))
         @test maximum(abs, u .- uex) < 1e-3
 
         # inhomogeneous Dirichlet folded through boundary_rhs: u = 1 on ∂Ω
@@ -349,11 +349,11 @@
             bc=ntuple(_ -> (Dirichlet(1.0), Dirichlet(1.0)), 2),
         )
         Li = -1 * laplacian(gi)
-        fi = flatten(set!(scalar_field(gi), x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2])))
+        fi = flatten(set!(x -> 2 * pi^2 * sinpi(x[1]) * sinpi(x[2]), scalar_field(gi)))
         rhs = fi .- flatten(boundary_rhs(Li, gi))
         ui = solve(MultigridSolver(Li), rhs; rtol=1e-10)
         uexi = flatten(
-            set!(scalar_field(gi), x -> 1 + sinpi(x[1]) * sinpi(x[2]))
+            set!(x -> 1 + sinpi(x[1]) * sinpi(x[2]), scalar_field(gi))
         )
         @test maximum(abs, ui .- uexi) < 4e-3
     end
diff --git a/test/ode_rhs.jl b/test/ode_rhs.jl
index 7c86141..4331a73 100644
--- a/test/ode_rhs.jl
+++ b/test/ode_rhs.jl
@@ -67,13 +67,13 @@ _grid_1d(n) = CartesianGrid(((0.0, 2π),), (n,); bc=((Periodic(), Periodic()),))
 
 function build_1d_leaf(n)
     g = _grid_1d(n)
-    u = set!(scalar_field(g), x -> sin(x[1]))
+    u = set!(x -> sin(x[1]), scalar_field(g))
     return similar(u), laplacian(g), u
 end
 
 function build_1d_added(n)
     g = _grid_1d(n)
-    u = set!(scalar_field(g), x -> sin(x[1]))
+    u = set!(x -> sin(x[1]), scalar_field(g))
     A = laplacian(g) + 2.0 * derivative(g, 1; order=1)          # Added(leaf, Scaled)
     return similar(u), A, u
 end
@@ -84,7 +84,7 @@ function build_2d(s::Int, prepared::Bool)
         ((0.0, 1.0), (0.0, 1.0)), (8s, 6s);
         bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic())),
     )
-    κ = set!(scalar_field(g), x -> 1 + x[1])
+    κ = set!(x -> 1 + x[1], scalar_field(g))
     K = divergence(g) * scaling(κ) * MatrixFreeOperators.gradient(g)
     u = scalar_field(g)
     interior(u) .= rand(Random.MersenneTwister(87 + s), local_size(g)...)
@@ -100,7 +100,7 @@ function build_3d_prepared(s::Int)
         ((0.0f0, 1.0f0), (0.0f0, 2.0f0), (0.0f0, 1.0f0)), (6s, 5s, 4s);
         bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic()), (Neumann(), Neumann())),
     )
-    κ = set!(scalar_field(g), x -> 1.0f0 + x[1] * x[3])
+    κ = set!(x -> 1.0f0 + x[1] * x[3], scalar_field(g))
     K = divergence(g) * scaling(κ) * MatrixFreeOperators.gradient(g) + 0.5f0 * laplacian(g)
     u = scalar_field(g)
     interior(u) .= rand(Random.MersenneTwister(87 + s), Float32, local_size(g)...)
@@ -114,7 +114,7 @@ function forest_step_alloc(n::Int, bs::Int)
     bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
     g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (n, n); bc=bc)
     bf = BlockForest(g; blocksize=(bs, bs), maxlevel=2)
-    uf = set!(scalar_field(bf), x -> sinpi(x[1]) * cospi(2x[2]) + 0.3 * x[1])
+    uf = set!(x -> sinpi(x[1]) * cospi(2x[2]) + 0.3 * x[1], scalar_field(bf))
     Sf = laplacian(bf) + adjoint(derivative(bf, 1; order=1))
     P = prepare(Sf)
     up = pack(uf)
@@ -130,13 +130,13 @@ end
         L = laplacian(g)
         f!(du, u) = apply!(du, L, u)
 
-        u = set!(scalar_field(g), x -> sin(x[1]))
+        u = set!(x -> sin(x[1]), scalar_field(g))
         dt = 0.005
         nsteps = 200
         rk4_field!(f!, u, dt, nsteps)
 
         t = dt * nsteps
-        uref = exp(-t) .* interior(set!(scalar_field(g), x -> sin(x[1])))
+        uref = exp(-t) .* interior(set!(x -> sin(x[1]), scalar_field(g)))
         err = maximum(abs, interior(u) .- uref)
         @info "field-level RK4 heat equation" t err
         @test err < 1e-3
@@ -160,7 +160,7 @@ end
         @test any(!iszero, interior(b))
         dt = 0.4 * spacing(g)[1]^2
         nsteps = 4000                                 # t = 6.25 ≫ 1/π², fully relaxed
-        exact = interior(set!(scalar_field(g), x -> 1 + x[1]))
+        exact = interior(set!(x -> 1 + x[1], scalar_field(g)))
 
         u = scalar_field(g)
         du = similar(u)
@@ -197,7 +197,7 @@ end
             ((0.0, 1.0), (0.0, 1.0)), (8, 6);
             bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic())),
         )
-        κ = set!(scalar_field(g), x -> 1 + x[1])
+        κ = set!(x -> 1 + x[1], scalar_field(g))
         K = divergence(g) * scaling(κ) * MatrixFreeOperators.gradient(g)   # Composed
         P = prepare(K, scalar_field(g))
 
@@ -240,7 +240,7 @@ end
             ((0.0f0, 1.0f0), (0.0f0, 2.0f0), (0.0f0, 1.0f0)), (6, 5, 4);
             bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic()), (Neumann(), Neumann())),
         )
-        κ3 = set!(scalar_field(g3), x -> 1.0f0 + x[1] * x[3])
+        κ3 = set!(x -> 1.0f0 + x[1] * x[3], scalar_field(g3))
         K3 = divergence(g3) * scaling(κ3) * MatrixFreeOperators.gradient(g3) + 0.5f0 * laplacian(g3)
         P3 = prepare(K3, scalar_field(g3))
         u3 = scalar_field(g3)
@@ -276,7 +276,7 @@ end
         bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)
-        uf = set!(scalar_field(bf), x -> sinpi(x[1]) * cospi(2x[2]) + 0.3 * x[1])
+        uf = set!(x -> sinpi(x[1]) * cospi(2x[2]) + 0.3 * x[1], scalar_field(bf))
         Df = derivative(bf, 1; order=1)
         Sf = laplacian(bf) + adjoint(Df)           # Added(leaf, PreparedAdjoint) after prepare
         P = prepare(Sf)
@@ -341,7 +341,7 @@ end
             bc=((Neumann(), Neumann()), (Neumann(), Neumann())),
         )
         P = prepare(laplacian(gd))
-        u = set!(scalar_field(gd), x -> x[1] + x[2]^2)
+        u = set!(x -> x[1] + x[2]^2, scalar_field(gd))
         du = similar(u)
         @test apply!(du, P, u) === du
 
@@ -372,9 +372,9 @@ end
         bf = BlockForest(g; blocksize=(4, 4), maxlevel=1)
         bf2 = BlockForest(g; blocksize=(4, 4), maxlevel=1)
         Pf = prepare(laplacian(bf))
-        uf = set!(scalar_field(bf), x -> x[1])
+        uf = set!(x -> x[1], scalar_field(bf))
         @test apply!(similar(uf), Pf, uf) isa BlockField
-        uf2 = set!(scalar_field(bf2), x -> x[1])
+        uf2 = set!(x -> x[1], scalar_field(bf2))
         @test_throws ArgumentError apply!(similar(uf2), Pf, uf2)
         @test_throws ArgumentError apply!(similar(uf2), Pf, uf)
     end
@@ -384,13 +384,13 @@ end
         P = prepare(laplacian(g))
         f!(du, u) = mul!(du, P, u)
 
-        u = flatten(set!(scalar_field(g), x -> sin(x[1])))
+        u = flatten(set!(x -> sin(x[1]), scalar_field(g)))
         dt = 0.005
         nsteps = 200
         rk4!(f!, u, dt, nsteps)
 
         t = dt * nsteps
-        uref = exp(-t) .* flatten(set!(scalar_field(g), x -> sin(x[1])))
+        uref = exp(-t) .* flatten(set!(x -> sin(x[1]), scalar_field(g)))
         @test maximum(abs, u .- uref) < 1e-3
     end
 end
diff --git a/test/operators_abstract.jl b/test/operators_abstract.jl
index 90f131c..ceffb9f 100644
--- a/test/operators_abstract.jl
+++ b/test/operators_abstract.jl
@@ -37,7 +37,7 @@ end
 
 @testset "Operator abstraction" begin
     g = CartesianGrid(((0.0, 1.0),), (4,))
-    u = set!(scalar_field(g), x -> x[1])
+    u = set!(x -> x[1], scalar_field(g))
 
     @testset "trait defaults make the weak claim" begin
         L = DummyNonlinearOp()
@@ -66,7 +66,7 @@ end
         @test collect(interior(L(u))) == collect(interior(y))
         @test collect(interior(L * u)) == collect(interior(y))
 
-        z = set!(scalar_field(g), x -> 1.0)
+        z = set!(x -> 1.0, scalar_field(g))
         MatrixFreeOperators.apply!(z, L, u, g, 3.0, 2.0)
         @test collect(interior(z)) ≈ 6 .* collect(interior(u)) .+ 2
     end
@@ -191,7 +191,7 @@ end
             a = @allocated apply_adjoint!(x̄, L, ȳ, g, 1.5, 2.0)
             return a, sum(interior(x̄))   # DCE-proof: consume the output
         end
-        κa = set!(scalar_field(ga), x -> 1 + x[1] * x[2])
+        κa = set!(x -> 1 + x[1] * x[2], scalar_field(ga))
         κa.data .= ifelse.(iszero.(κa.data), one(eltype(κa.data)), κa.data)   # incl. ghosts
         for L in (
             laplacian(ga),
@@ -199,7 +199,7 @@ end
             laplacian(ga) + derivative(ga, 2; order=2),
             MatrixFreeOperators.Diffusion(ga, κa, ArithmeticMean()),
         )
-            ȳ = set!(scalar_field(ga), x -> sinpi(x[1]) * exp(-x[2]))
+            ȳ = set!(x -> sinpi(x[1]) * exp(-x[2]), scalar_field(ga))
             a, s = alloc_adjoint(L, scalar_field(ga), ȳ, ga)
             @test isfinite(s) && !iszero(s)
             @test a ≤ 512
diff --git a/test/packedfield.jl b/test/packedfield.jl
index 62d3a58..8f8ad40 100644
--- a/test/packedfield.jl
+++ b/test/packedfield.jl
@@ -6,7 +6,7 @@
     bc = ((Dirichlet(), Dirichlet()), (Neumann(), Neumann()))
     g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16); bc=bc)
     bf = BlockForest(g; blocksize=(4, 4), maxlevel=2)   # 16 leaves
-    uf = set!(scalar_field(bf), fun)
+    uf = set!(fun, scalar_field(bf))
 
     @testset "pack/unpack round trip is bit-exact" begin
         p = pack(uf)
@@ -23,7 +23,7 @@
             @test up.blocks[i] == uf.blocks[i]
         end
 
-        w = set!(vector_field(bf), vfun)                # SVector eltype packs too
+        w = set!(vfun, vector_field(bf))                # SVector eltype packs too
         pw = pack(w)
         @test eltype(pw) === SVector{2,Float64}
         for i in 1:MFO.nleaves(bf)
@@ -32,7 +32,7 @@
     end
 
     @testset "set! on packed matches pack ∘ set!" begin
-        p = set!(MFO._zero_all!(pack(scalar_field(bf))), fun)
+        p = set!(fun, MFO._zero_all!(pack(scalar_field(bf))))
         @test p.data == pack(uf).data
     end
 
@@ -74,7 +74,7 @@
         bf2 = BlockForest(
             CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8)); blocksize=(4, 4), maxlevel=2
         )
-        u2 = set!(scalar_field(bf2), fun)
+        u2 = set!(fun, scalar_field(bf2))
         p2 = pack(u2)
         @test_throws ArgumentError regrid!(p2; refine=Returns(true))         # packed never regrids
         @test_throws ArgumentError regrid!(u2, p2; refine=Returns(true))     # nor mixed in varargs
diff --git a/test/partitioning.jl b/test/partitioning.jl
index 3620ffb..addccc1 100644
--- a/test/partitioning.jl
+++ b/test/partitioning.jl
@@ -209,7 +209,7 @@ function dist_boundary_rhs(D, T=Float64)
     return b
 end
 
-# CPU twin of the extension's set!(::MultiDeviceVector, P, fun).
+# CPU twin of the extension's set!(fun, ::MultiDeviceVector, P).
 function dist_set(D, fun, T=Float64)
     MatrixFreeOperators._dist_set!(D.xs, fun, D.ctx)
     x = Vector{T}(undef, sum(length, D.owned))
@@ -293,10 +293,10 @@ end
         # global one bit for bit, so a distributed RHS is partition-independent.
         g = CartesianGrid(((0.3, 1.7), (-1.1, 2.9)), (5, 9))
         fun = x -> sin(3x[1]) * exp(-x[2]) + 0.25x[1] * x[2]
-        ref = set!(scalar_field(g), fun)
+        ref = set!(fun, scalar_field(g))
         for np in (2, 3)
             for lp in partition_grid(g, np)
-                loc = set!(scalar_field(lp), fun)
+                loc = set!(fun, scalar_field(lp))
                 @test collect(interior(loc)) ==
                     collect(view(interior(ref), lp.local_range...))
             end
@@ -751,7 +751,7 @@ end
         g = gridof((Dirichlet(), Neumann()))
         n = prod(local_size(g))
         x = rand(MersenneTwister(11), n)
-        κ = set!(scalar_field(g), coeff_fun)
+        κ = set!(coeff_fun, scalar_field(g))
         for L in (scaling(κ), laplacian(g) * scaling(κ), scaling(κ) * laplacian(g))
             D = dist_prepare(L, g, 2)
             clean = dist_mul(D, x)
@@ -772,7 +772,7 @@ end
             g = gridof(cut)
             n = prod(local_size(g))
             x = rand(rng, n)
-            κ = set!(scalar_field(g), coeff_fun)
+            κ = set!(coeff_fun, scalar_field(g))
             ops = (
                 scaling(κ),
                 laplacian(g) * scaling(κ),
@@ -793,7 +793,7 @@ end
         g = gridof((Dirichlet(), Neumann()))
         n = prod(local_size(g))
         x = rand(MersenneTwister(13), n)
-        κ = set!(scalar_field(g), coeff_fun)
+        κ = set!(coeff_fun, scalar_field(g))
         D = dist_prepare(laplacian(g) * scaling(κ), g, 2)
         good = dist_mul(D, x)
         @test good == flatten(apply(laplacian(g) * scaling(κ), _field(g, x)))
@@ -813,7 +813,7 @@ end
         for cut in cutbcs, np in (2, 3)
             g = gridof(cut)
             n = prod(local_size(g))
-            κ = set!(scalar_field(g), coeff_fun)
+            κ = set!(coeff_fun, scalar_field(g))
             x, y = rand(rng, n), rand(rng, n)
             for L in (
                 scaling(κ) * laplacian(g),
@@ -870,7 +870,7 @@ end
     @testset "boundary_rhs parity" begin
         for cut in inhom_cuts, np in (1, 2, 3)
             g = inhom_grid(cut)
-            κ = set!(scalar_field(g), coeff_fun)
+            κ = set!(coeff_fun, scalar_field(g))
             D1 = derivative(g, 1)
             ops = (
                 laplacian(g),
@@ -962,11 +962,11 @@ end
     @testset "a full inhomogeneous RHS assembles slab-locally" begin
         for cut in inhom_cuts, np in (2, 3)
             g = inhom_grid(cut, (8, 12))
-            κ = set!(scalar_field(g), coeff_fun)
+            κ = set!(coeff_fun, scalar_field(g))
             fun = x -> sin(3x[1]) * exp(-x[2]) + 0.25x[1] * x[2]
             for L in (laplacian(g), scaling(κ) * laplacian(g))
                 D = dist_prepare(L, g, np)
-                ref = flatten(set!(scalar_field(g), fun)) .- flatten(boundary_rhs(L, g))
+                ref = flatten(set!(fun, scalar_field(g))) .- flatten(boundary_rhs(L, g))
                 @test dist_set(D, fun) .- dist_boundary_rhs(D) == ref
             end
         end
@@ -1002,7 +1002,7 @@ end
         for (ext, sz) in diff_exts, cut in ((Dirichlet(), Neumann()), (Periodic(), Periodic()))
             g = diff_grid(ext, sz, cut)
             N = length(sz)
-            plain = set!(scalar_field(g), diff_coeff_fun)
+            plain = set!(diff_coeff_fun, scalar_field(g))
             extended = diffusion(g, plain).κ
             for (κ, ghosts_read) in ((plain, false), (extended, true)), np in (1, 2, 3)
                 parts = partition_grid(g, np)
@@ -1043,7 +1043,7 @@ end
 
             g = diff_grid(ext, sz, cut)
             n = prod(local_size(g))
-            κ = set!(scalar_field(g), diff_coeff_fun)
+            κ = set!(diff_coeff_fun, scalar_field(g))
             Dop = diffusion(g, κ; averaging=avg)
             x, y = rand(rng, n), rand(rng, n)
             # `laplacian(g) + Dop` puts the leaf SECOND under the Added, so the slab
@@ -1081,7 +1081,7 @@ end
             g = diff_grid(ext, sz, cut)
             n = prod(local_size(g))
             x = rand(MersenneTwister(23), n)
-            κ = set!(scalar_field(g), diff_coeff_fun)
+            κ = set!(diff_coeff_fun, scalar_field(g))
             for L in (diffusion(g, κ), laplacian(g) * diffusion(g, κ))
                 D = dist_prepare(L, g, 2)
                 clean = dist_mul(D, x)
@@ -1107,7 +1107,7 @@ end
     # therefore be bit-for-bit the Laplacian's — same gating rule, same nodes.
     @testset "Diffusion adds no per-apply exchange" begin
         g = gridof((Dirichlet(), Neumann()))
-        κ = set!(scalar_field(g), diff_coeff_fun)
+        κ = set!(diff_coeff_fun, scalar_field(g))
         Dop = diffusion(g, κ)
         # A bare stencil leaf: no mid-tree node exists to hold an exchange.
         @test dist_prepare(Dop, g, 2).tree isa MatrixFreeOperators.DistLeaf
@@ -1142,10 +1142,10 @@ end
         for mk in (
             g -> laplacian(g) * laplacian(g),
             # A localized leaf dispatches dynamically; that must stay O(1), not O(cells).
-            g -> laplacian(g) * scaling(set!(scalar_field(g), coeff_fun)),
+            g -> laplacian(g) * scaling(set!(coeff_fun, scalar_field(g))),
             # ...including the one whose localization copies a padded coefficient:
             # that copy belongs to prepare, and must not reappear per apply.
-            g -> laplacian(g) * diffusion(g, set!(scalar_field(g), diff_coeff_fun)),
+            g -> laplacian(g) * diffusion(g, set!(diff_coeff_fun, scalar_field(g))),
         )
             small, large = steady((16, 16), mk, dist_mul!), steady((32, 32), mk, dist_mul!)
             @test large < 2 * small
@@ -1165,9 +1165,9 @@ end
             # ...and the slab diffusion leaf, whose adjoint is the interior stencil
             # plus the ghost-plane gather over a padded κ (issue #77) rather than the
             # self-adjoint shortcut.
-            g -> diffusion(g, set!(scalar_field(g), diff_coeff_fun)) + laplacian(g),
+            g -> diffusion(g, set!(diff_coeff_fun, scalar_field(g))) + laplacian(g),
             # ...in both Added slots: as `node.b` the leaf accumulates with β = true.
-            g -> laplacian(g) + diffusion(g, set!(scalar_field(g), diff_coeff_fun)),
+            g -> laplacian(g) + diffusion(g, set!(diff_coeff_fun, scalar_field(g))),
         )
             small = steady((16, 16), mk, dist_adjoint!)
             large = steady((32, 32), mk, dist_adjoint!)
@@ -1200,7 +1200,7 @@ end
 
             # slice 2b: a real coefficient field on an undistributed CartesianGrid
             # is sliceable onto the slabs, so it joins the whitelist.
-            κ = set!(scalar_field(g), x -> 1 + x[1])
+            κ = set!(x -> 1 + x[1], scalar_field(g))
             @test distributable(scaling(κ))
             @test distributable(scaling(κ) + laplacian(g))
             @test distributable(laplacian(g) * scaling(κ))
@@ -1210,7 +1210,7 @@ end
             # terms. Its face averaging reads κ across the cut, which
             # `_slab_field`'s padded window supplies at localization time — no exchange, so
             # nothing further to require of it here.
-            κp = set!(scalar_field(g), x -> 1 + x[1] + x[2])
+            κp = set!(x -> 1 + x[1] + x[2], scalar_field(g))
             for avg in (ArithmeticMean(), HarmonicMean())
                 @test distributable(diffusion(g, κp; averaging=avg))
             end
@@ -1220,7 +1220,7 @@ end
         end
 
         @testset "rejected: field-valued parameters" begin
-            v = set!(vector_field(g), x -> SVector(1.0, 0.0))
+            v = set!(x -> SVector(1.0, 0.0), vector_field(g))
             @test !distributable(advection(g, v))
             # the message must name the reason, not just the type
             err = try
@@ -1271,11 +1271,11 @@ end
             @test all(p -> !same(p, g), partition_grid(g, 2))
             @test same(Adapt.adapt(Array, g), g)
 
-            κg = set!(scalar_field(g), x -> 1 + x[1])
+            κg = set!(x -> 1 + x[1], scalar_field(g))
             @test check1(laplacian(g) + scaling(κg), g) isa MatrixFreeOperators.Added
             @test check1(laplacian(g) + scaling(κg), g2) isa MatrixFreeOperators.Added
 
-            κc = set!(scalar_field(gc), x -> 1 + x[1])
+            κc = set!(x -> 1 + x[1], scalar_field(gc))
             @test distributable(laplacian(g) + scaling(κc))        # each leaf is fine alone
             err = try
                 check1(laplacian(g) + scaling(κc), g)
@@ -1308,7 +1308,7 @@ end
             # face coefficient, no error), or a `BoundsError` on a smaller one.
             for szκ in ((12, 10), (4, 6))
                 gκ = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), szκ)
-                κκ = set!(scalar_field(gκ), x -> 1 + x[1])
+                κκ = set!(x -> 1 + x[1], scalar_field(gκ))
                 Dbad = MatrixFreeOperators.Diffusion(g, κκ, ArithmeticMean())
                 @test distributable(Dbad)          # the leaf alone cannot tell
                 for L in (Dbad, laplacian(g) + Dbad)
@@ -1359,7 +1359,7 @@ end
 
         # The error points at the offending node, not merely at the tree root.
         @testset "message names the culprit inside a tree" begin
-            v = set!(vector_field(g), x -> SVector(1.0, 0.0))
+            v = set!(x -> SVector(1.0, 0.0), vector_field(g))
             err = try
                 check(laplacian(g) + 2.0 * advection(g, v))
             catch e
diff --git a/test/prepare_linalg.jl b/test/prepare_linalg.jl
index cd7fb22..931c37e 100644
--- a/test/prepare_linalg.jl
+++ b/test/prepare_linalg.jl
@@ -69,11 +69,11 @@
         function poisson_error(n)
             g = CartesianGrid(((0.0, 1.0),), (n,))
             A = prepare(laplacian(g))
-            f = set!(scalar_field(g), x -> π^2 * sin(π * x[1]))
+            f = set!(x -> π^2 * sin(π * x[1]), scalar_field(g))
             b = -flatten(f)
             u, stats = Krylov.minres(A, b)
             @test stats.solved
-            ustar = flatten(set!(scalar_field(g), x -> sin(π * x[1])))
+            ustar = flatten(set!(x -> sin(π * x[1]), scalar_field(g)))
             return maximum(abs, u .- ustar)
         end
         e32 = poisson_error(32)
@@ -101,7 +101,7 @@
             ((0.0, 1.0), (0.0, 1.0)), (8, 6);
             bc=((Dirichlet(), Dirichlet()), (Periodic(), Periodic())),
         )
-        κ = set!(scalar_field(g), x -> 1 + x[1])
+        κ = set!(x -> 1 + x[1], scalar_field(g))
         K = divergence(g) * scaling(κ) * MatrixFreeOperators.gradient(g)
         P = prepare(K, scalar_field(g))
         @test P.op isa MatrixFreeOperators.PreparedComposed
@@ -126,7 +126,7 @@
             ((0.0, 1.0), (0.0, 1.0)), (5, 4);
             bc=((Dirichlet(), Neumann()), (Periodic(), Periodic())),
         )
-        v = set!(vector_field(g), x -> SVector(1 + x[1], x[2]))
+        v = set!(x -> SVector(1 + x[1], x[2]), vector_field(g))
         Adv = advection(g, v)
         P = prepare(Adv, scalar_field(g))
         Pt = prepare(adjoint(Adv), scalar_field(g))
@@ -163,7 +163,7 @@
         @test P.op isa MatrixFreeOperators.PreparedComposed
         @test materialize(P) ≈ materialize(prepare(inner))'
 
-        κ = set!(scalar_field(g), x -> 1 + x[1] * x[2])
+        κ = set!(x -> 1 + x[1] * x[2], scalar_field(g))
         K = divergence(g) * scaling(κ)                           # vector → scalar
         Kt = prepare(MatrixFreeOperators.AdjointOp(K), scalar_field(g))
         B = materialize(prepare(K, vector_field(g)))
@@ -190,14 +190,14 @@
             σf(x) = 1 + x[1] * x[2]
             ustar(x) = sin(π * x[1]) * sin(π * x[2])
             f(x) = 2 * π^2 * ustar(x) + σf(x) * ustar(x)
-            σ = set!(scalar_field(g), σf)
+            σ = set!(σf, scalar_field(g))
             K = scaling(σ) - laplacian(g)
             @test isselfadjoint(K)
             P = prepare(K, scalar_field(g))
-            b = flatten(set!(scalar_field(g), f))
+            b = flatten(set!(f, scalar_field(g)))
             u, stats = Krylov.cg(P, b)
             @test stats.solved
-            return maximum(abs, u .- flatten(set!(scalar_field(g), ustar)))
+            return maximum(abs, u .- flatten(set!(ustar, scalar_field(g))))
         end
         e16 = helmholtz_error(16)
         e32 = helmholtz_error(32)
@@ -214,13 +214,13 @@
         @test all(iszero, collect(interior(apply(L, zero_in))))      # islinear ⇒ L(0) = 0
 
         b = boundary_rhs(L, g)
-        f = set!(scalar_field(g), x -> π^2 * sin(π * x[1]))
+        f = set!(x -> π^2 * sin(π * x[1]), scalar_field(g))
         rhs = -flatten(f) .- flatten(b)                               # Δu = -f  ⇒  A·u = -f - b
         P = prepare(L)
         u, stats = Krylov.minres(P, rhs)
         @test stats.solved
         ustar = flatten(
-            set!(scalar_field(g), x -> sin(π * x[1]) + (1 - x[1]) * a + x[1] * c)
+            set!(x -> sin(π * x[1]) + (1 - x[1]) * a + x[1] * c, scalar_field(g))
         )
         @test maximum(abs, u .- ustar) < 0.01
     end
@@ -240,9 +240,9 @@
             ((0.0, 1.0), (0.0, 1.0)), (6, 5);
             bc=((Dirichlet(), Neumann()), (Periodic(), Periodic())),
         )
-        κ = set!(scalar_field(g), x -> 1 + x[1] * x[2])
-        v = set!(vector_field(g), x -> SVector(1 + x[1], x[2]))
-        u0 = set!(scalar_field(g), x -> sin(π * x[1]))
+        κ = set!(x -> 1 + x[1] * x[2], scalar_field(g))
+        v = set!(x -> SVector(1 + x[1], x[2]), vector_field(g))
+        u0 = set!(x -> sin(π * x[1]), scalar_field(g))
         Lap = laplacian(g)
         Adv = advection(g, v)
         S = scaling(κ)
@@ -323,8 +323,8 @@
         # pins the coefficient's type so the rewrite is actually exercised, then holds
         # the rebuilt tree to the traits of the tree handed in — including the ones
         # that are false (isdiagonal of Diffusion, isselfadjoint of Advection).
-        κf = set!(scalar_field(bf), x -> 1 + x[1] * x[2])
-        vf = set!(vector_field(bf), x -> SVector(1 + x[1], x[2]))
+        κf = set!(x -> 1 + x[1] * x[2], scalar_field(bf))
+        vf = set!(x -> SVector(1 + x[1], x[2]), vector_field(bf))
         Sf = scaling(κf)
         Dif = diffusion(bf, κf)
         Af = advection(bf, vf)
@@ -379,13 +379,13 @@
         @test MFO.adjoint_operator(P.op) === D1
         # a diagonal leaf wrapped by hand — no prepare path builds one, but the
         # declared transpose must still be the leaf, not a lazy wrapper of the twin
-        S = scaling(set!(scalar_field(g), x -> 1 + x[1]))
+        S = scaling(set!(x -> 1 + x[1], scalar_field(g)))
         Sᵀ = MFO.PreparedAdjoint(S, scalar_field(g))
         @test isdiagonal(Sᵀ)
         @test MFO.adjoint_operator(Sᵀ) === S
         @test !(MFO.adjoint_operator(Sᵀ) isa AdjointOp)
 
-        ȳ = set!(scalar_field(g), x -> sin(2π * x[2]) * (1 - x[1]) + x[1]^2)
+        ȳ = set!(x -> sin(2π * x[2]) * (1 - x[1]) + x[1]^2, scalar_field(g))
         x̄ = scalar_field(g)
         apply_adjoint!(x̄, P.op, ȳ, g)
         yref = apply(D1, ȳ)
@@ -399,7 +399,7 @@
         )
         bf = BlockForest(gp; blocksize=(4, 4), maxlevel=2)
         Df = derivative(bf, 1; order=1)
-        Sf = scaling(set!(scalar_field(bf), x -> 1 + x[1] * x[2]))
+        Sf = scaling(set!(x -> 1 + x[1] * x[2], scalar_field(bf)))
         # PreparedAdjoint(Df) alone, and nested as a factor of a PreparedComposed:
         # (Sfᵀ ∘ Dfᵀ)ᵀ = Df ∘ Sf — the transpose walk recurses into the twin.
         for (label, Lt, Lref) in (
@@ -428,7 +428,7 @@
     @testset "boundary_rhs through combinators" begin
         g = CartesianGrid(((0.0, 1.0),), (8,); bc=((Dirichlet(2.0), Neumann(1.0)),))
         L = laplacian(g)
-        S = scaling(set!(scalar_field(g), x -> 1 + x[1]))
+        S = scaling(set!(x -> 1 + x[1], scalar_field(g)))
         bL = collect(interior(boundary_rhs(L, g)))
         @test collect(interior(boundary_rhs(3 * L, g))) ≈ 3 .* bL
         @test collect(interior(boundary_rhs(L + L, g))) ≈ 2 .* bL
diff --git a/test/reactant_parity.jl b/test/reactant_parity.jl
index ff2f768..6c5249a 100644
--- a/test/reactant_parity.jl
+++ b/test/reactant_parity.jl
@@ -10,7 +10,7 @@ Reactant.set_default_backend("cpu")
             bc=((Periodic(), Periodic()), (Dirichlet(), Neumann())),
         )
         P = prepare(laplacian(g), scalar_field(g))
-        x = flatten(set!(scalar_field(g), x -> sin(x[1]) * x[2]^2))
+        x = flatten(set!(x -> sin(x[1]) * x[2]^2, scalar_field(g)))
         y = zero(x)
         mul!(y, P, x)
 
@@ -25,7 +25,7 @@ Reactant.set_default_backend("cpu")
     @testset "prepared mul! parity (3D)" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0), (0.0, 1.0)), (8, 9, 10))
         P = prepare(laplacian(g), scalar_field(g))
-        x = flatten(set!(scalar_field(g), x -> x[1]^2 + sinpi(x[2]) * x[3]))
+        x = flatten(set!(x -> x[1]^2 + sinpi(x[2]) * x[3], scalar_field(g)))
         y = zero(x)
         mul!(y, P, x)
 
@@ -40,8 +40,8 @@ Reactant.set_default_backend("cpu")
     @testset "apply! α/β accumulation parity" begin
         g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16))
         L = laplacian(g)
-        u = set!(scalar_field(g), x -> sinpi(x[1]) * x[2])
-        w = set!(scalar_field(g), x -> x[1] + x[2])
+        u = set!(x -> sinpi(x[1]) * x[2], scalar_field(g))
+        w = set!(x -> x[1] + x[2], scalar_field(g))
         y = apply!(deepcopy(w), L, deepcopy(u), g, 2.0, 0.5)
 
         ur = Reactant.to_rarray(deepcopy(u))
diff --git a/test/scaling_identity_advection.jl b/test/scaling_identity_advection.jl
index 2adb032..30a6ebb 100644
--- a/test/scaling_identity_advection.jl
+++ b/test/scaling_identity_advection.jl
@@ -1,7 +1,7 @@
 @testset "ScalingOp, IdentityOp, Advection" begin
     @testset "scaling by a Number" begin
         g = CartesianGrid(((0.0, 1.0),), (6,))
-        u = set!(scalar_field(g), x -> x[1])
+        u = set!(x -> x[1], scalar_field(g))
         S = scaling(2.5)
         @test collect(interior(S * u)) ≈ 2.5 .* collect(interior(u))
         @test islinear(S) && isconstant(S) && isdiagonal(S) && isselfadjoint(S)
@@ -11,14 +11,14 @@
 
     @testset "scaling by a coefficient field" begin
         g = CartesianGrid(((0.0, 1.0),), (6,))
-        κ = set!(scalar_field(g), x -> 1 + x[1]^2)
-        u = set!(scalar_field(g), x -> sin(x[1]))
+        κ = set!(x -> 1 + x[1]^2, scalar_field(g))
+        u = set!(x -> sin(x[1]), scalar_field(g))
         S = scaling(κ)
         @test collect(interior(S * u)) ≈ collect(interior(κ)) .* collect(interior(u))
         @test isselfadjoint(S) && adjoint(S) === S
         @test MatrixFreeOperators.operator_grid(S) === g
 
-        v = set!(vector_field(g), x -> SVector(x[1]))
+        v = set!(x -> SVector(x[1]), vector_field(g))
         Sv = S * v
         @test getindex.(collect(interior(Sv)), 1) ≈
             collect(interior(κ)) .* getindex.(collect(interior(v)), 1)
@@ -28,12 +28,12 @@
 
     @testset "identity_op" begin
         g = CartesianGrid(((0.0, 1.0),), (6,))
-        u = set!(scalar_field(g), x -> x[1]^3)
+        u = set!(x -> x[1]^3, scalar_field(g))
         I = identity_op()
         @test collect(interior(I * u)) == collect(interior(u))
         @test islinear(I) && isselfadjoint(I) && isdiagonal(I)
         @test adjoint(I) === I
-        z = set!(scalar_field(g), x -> 1.0)
+        z = set!(x -> 1.0, scalar_field(g))
         MatrixFreeOperators.apply!(z, I, u, g, 2.0, -1.0)
         @test collect(interior(z)) ≈ 2 .* collect(interior(u)) .- 1
     end
@@ -44,15 +44,15 @@
                 ((0.0, 2π), (0.0, 2π)), (n, n);
                 bc=((Periodic(), Periodic()), (Periodic(), Periodic())),
             )
-            v = set!(vector_field(g), x -> SVector(sin(x[2]), cos(x[1])))
-            u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
+            v = set!(x -> SVector(sin(x[2]), cos(x[1])), vector_field(g))
+            u = set!(x -> sin(x[1]) * sin(x[2]), scalar_field(g))
             A = advection(g, v)
             @test islinear(A) && isconstant(A)
             y = A * u
             ref = set!(
-                scalar_field(g),
                 x ->
                     sin(x[2]) * cos(x[1]) * sin(x[2]) + cos(x[1]) * sin(x[1]) * cos(x[2]),
+                scalar_field(g),
             )
             return maximum(abs, collect(interior(y)) .- collect(interior(ref)))
         end
@@ -74,9 +74,9 @@
         @test !islinear(A)
         @test_throws ArgumentError adjoint(A)
 
-        u = set!(vector_field(g), x -> SVector(sin(x[1])))
+        u = set!(x -> SVector(sin(x[1])), vector_field(g))
         y = A * u
-        ref = set!(vector_field(g), x -> SVector(sin(x[1]) * cos(x[1])))
+        ref = set!(x -> SVector(sin(x[1]) * cos(x[1])), vector_field(g))
         @test maximum(norm.(collect(interior(y)) .- collect(interior(ref)))) < 0.01
 
         @test_throws ArgumentError apply(A, scalar_field(g))