diff --git a/mosaic/mosaic_shared/problems/navier_stokes_grid/schemas.py b/mosaic/mosaic_shared/problems/navier_stokes_grid/schemas.py
index 827e64ad..7bd8d430 100644
--- a/mosaic/mosaic_shared/problems/navier_stokes_grid/schemas.py
+++ b/mosaic/mosaic_shared/problems/navier_stokes_grid/schemas.py
@@ -126,6 +126,20 @@ class InputSchema(BaseModel):
             "Requires x_lo to be 'dirichlet' with no value."
         ),
     )
+    forcing: GridVectorField | None = Field(
+        default=None,
+        description=(
+            "Body force per unit mass (an acceleration) at the cell centers, "
+            "constant in time, same shape as v0. 2-D periodic boxes only."
+        ),
+    )
+    linear_drag: Array[(1,), Float32] | None = Field(
+        default=None,
+        description=(
+            "Coefficient α of a linear drag -αu on the right-hand side of the "
+            "momentum equation. 2-D periodic boxes only."
+        ),
+    )
 
     @property
     def extent(self) -> tuple[float, ...]:
@@ -179,6 +193,18 @@ def _check_bcs(self) -> "InputSchema":
                     f"inflow_profile has shape {tuple(self.inflow_profile.shape)}, "
                     f"expected ({ny},) to match v0"
                 )
+        if self.forcing is not None or self.linear_drag is not None:
+            if not bc.is_fully_periodic or self.obstacle is not None:
+                raise ValueError(
+                    "forcing and linear_drag require a periodic box without an obstacle"
+                )
+            if self.v0.shape[-1] != 2:
+                raise ValueError("forcing and linear_drag are only defined in 2-D")
+        if self.forcing is not None and self.forcing.shape != self.v0.shape:
+            raise ValueError(
+                f"forcing has shape {tuple(self.forcing.shape)}, expected "
+                f"{tuple(self.v0.shape)} to match v0"
+            )
         return self
 
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/exponax/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/exponax/tesseract_api.py
index c0504774..11a517ab 100644
--- a/mosaic/tesseracts/navier-stokes-grid/exponax/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/exponax/tesseract_api.py
@@ -49,6 +49,8 @@ class InputSchema(
     @model_validator(mode="after")
     def _check_periodic_bcs(self) -> "InputSchema":
         check_periodic(self, "exponax")
+        if self.forcing is not None and self.kolmogorov_forcing:
+            raise ValueError("pass either forcing or kolmogorov_forcing, not both")
         return self
 
 
@@ -113,13 +115,24 @@ def exponax_fwd(
     injection_mode: int,
     injection_scale: float,
     boundary_conditions: dict | None = None,
+    forcing: jnp.ndarray | None = None,
+    linear_drag: float | None = None,
     **_kwargs: Any,
 ) -> jnp.ndarray:
     """Run 2D or 3D incompressible Navier-Stokes using exponax.
 
     2D uses the pseudo-spectral streamfunction-vorticity stepper with an
     internal velocity<->vorticity conversion so the interface matches the
-    other CFD tesseracts.
+    other CFD tesseracts. Exponax's 2D steppers evolve q = -ω: they take the
+    stream function from Δψ = q and the velocity from (∂ᵧψ, -∂ₓψ), whose
+    vorticity is -q.
+
+    A ``forcing`` field enters as the injection term of Exponax's Kolmogorov
+    stepper, the spectral curl of the force on the grid as given. Its own
+    ``injection_mode``/``injection_scale`` cannot express the canonical force:
+    Exponax puts sample j at y = j·h, the canonical layout at (j + ½)h, so its
+    sin(κy) is half a cell off. ``linear_drag`` α maps to Exponax's drag λ,
+    which adds λω and so damps for λ < 0.
 
     3D uses the Leray-projected velocity stepper directly.
 
@@ -135,6 +148,8 @@ def exponax_fwd(
         injection_mode: Forcing wavenumber.
         injection_scale: Forcing amplitude.
         boundary_conditions: Boundary condition config (must be fully periodic).
+        forcing: Body force at the grid points, shape of v0 (2D only).
+        linear_drag: Coefficient α of a linear drag -αu (2D only).
 
     Returns:
         Final velocity field, same shape as v0.
@@ -152,15 +167,18 @@ def exponax_fwd(
     }
 
     if ndim == 2:
+        if linear_drag is not None:
+            stepper_kwargs["drag"] = drag - linear_drag
+
         # v0: (N, N, 1, 2) -> squeeze z -> (N, N, 2)
         v = v0[:, :, 0, :]
 
-        # Mean velocity is conserved by incompressible periodic NS (Galilean invariance).
-        # The vorticity formulation strips it (∂const/∂x = 0), so we save and restore it.
+        # The vorticity formulation strips the mean velocity (∂const/∂x = 0),
+        # so we save it and evolve it separately.
         v_mean = jnp.mean(v, axis=(0, 1), keepdims=True)  # (1, 1, 2)
 
-        # velocity -> vorticity: (1, N, N)  (operate on fluctuation only)
-        vort = _velocity_to_vorticity(v - v_mean, domain_extent)
+        # velocity -> Exponax's q = -ω: (1, N, N)  (operate on fluctuation only)
+        q = -_velocity_to_vorticity(v - v_mean, domain_extent)
 
         if kolmogorov_forcing:
             stepper = ex.stepper.KolmogorovFlowVorticity(
@@ -169,6 +187,18 @@ def exponax_fwd(
                 injection_scale=injection_scale,
                 **stepper_kwargs,
             )
+        elif forcing is not None:
+            stepper = ex.stepper.KolmogorovFlowVorticity(
+                num_spatial_dims=2, injection_scale=0.0, **stepper_kwargs
+            )
+            # The force drives q with -curl f. Exponax builds its injection as
+            # a real array; the force's coefficients need their phase.
+            curl = _velocity_to_vorticity(forcing[:, :, 0, :], domain_extent)
+            stepper = eqx.tree_at(
+                lambda st: st._integrator._nonlinear_fun.injection,
+                stepper,
+                jnp.fft.rfft2(-curl).astype(jnp.complex64),
+            )
         else:
             stepper = ex.stepper.NavierStokesVorticity(
                 num_spatial_dims=2, **stepper_kwargs
@@ -177,10 +207,25 @@ def exponax_fwd(
         def step_fn(carry: Any, _: Any) -> tuple[Any, None]:
             return stepper(carry).astype(carry.dtype), None
 
-        result_vort, _ = jax.lax.scan(step_fn, vort, None, length=steps)
+        result_q, _ = jax.lax.scan(step_fn, q, None, length=steps)
 
-        # vorticity -> velocity: (N, N, 2) -> (N, N, 1, 2); restore conserved mean
-        result = _vorticity_to_velocity(result_vort, domain_extent) + v_mean
+        if forcing is not None or linear_drag is not None:
+            # dū/dt = λū + f̄ under Exponax's drag λ and the mean force f̄.
+            rate = stepper_kwargs["drag"]
+            t_end = dt * steps
+            f_mean = (
+                0.0
+                if forcing is None
+                else jnp.mean(forcing[:, :, 0, :], axis=(0, 1), keepdims=True)
+            )
+            growth = jnp.exp(rate * t_end)
+            # (e^{λT} - 1)/λ, which tends to T as λ → 0 (rate may be traced).
+            safe = jnp.where(rate == 0, 1.0, rate)
+            gain = jnp.where(rate == 0, t_end, (growth - 1.0) / safe)
+            v_mean = v_mean * growth + f_mean * gain
+
+        # q -> velocity: (N, N, 2) -> (N, N, 1, 2); restore the mean
+        result = _vorticity_to_velocity(-result_q, domain_extent) + v_mean
         return result[:, :, None, :]
 
     else:  # 3D
@@ -219,13 +264,13 @@ def apply_jit(inputs: dict) -> dict:
     return {"result": exponax_fwd(**inputs)}
 
 
-_SCALAR_KEYS = ("dt", "viscosity", "drag", "injection_scale")
+_SCALAR_KEYS = ("dt", "viscosity", "drag", "injection_scale", "linear_drag")
 
 
 def _unpack_scalars(d: dict) -> dict:
     """Extract Python floats from 1-element arrays for JIT-static scalar params."""
     for key in _SCALAR_KEYS:
-        if key in d:
+        if d.get(key) is not None:
             d[key] = float(d[key][0])
     return d
 
@@ -244,7 +289,8 @@ def _diff_apply(inputs: dict) -> dict:
     """Forward pass differentiated by jax_vjp, with the scalar params traced."""
     fwd_kwargs = dict(inputs)
     for k in _SCALAR_KEYS:
-        fwd_kwargs[k] = jnp.reshape(fwd_kwargs[k], ())
+        if fwd_kwargs.get(k) is not None:
+            fwd_kwargs[k] = jnp.reshape(fwd_kwargs[k], ())
     return {"result": exponax_fwd(**fwd_kwargs)}
 
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/ins-jl/solver.jl b/mosaic/tesseracts/navier-stokes-grid/ins-jl/solver.jl
index 14c13c13..d7feab66 100644
--- a/mosaic/tesseracts/navier-stokes-grid/ins-jl/solver.jl
+++ b/mosaic/tesseracts/navier-stokes-grid/ins-jl/solver.jl
@@ -209,6 +209,43 @@ function ns_forward_3d(v0::AbstractArray, rhs, setup, psolver,
 end
 
 
+# ---------------------------------------------------------------------------
+# Body force and linear drag (2-D)
+# ---------------------------------------------------------------------------
+
+"""
+Collocated body force (n,n,2) → projected, ghosted staggered force, or `nothing`.
+
+The force is averaged onto the faces like v0 and projected once, so that
+adding it to the projected right-hand side keeps the stages divergence-free.
+"""
+function staggered_force(force, setup, psolver, n::Int)
+    isnothing(force) && return nothing
+    f = add_ghosts_2d(coloc_to_stag_2d(to_device(Float32.(force)), n), n)
+    return project(f, setup; psolver)
+end
+
+"""
+`rhs` plus a constant body force `f` and a linear drag `-drag * u`.
+
+The drag acts on the stage velocity, which is divergence-free, so the sum
+stays projected. Returns `rhs` itself when there is neither.
+"""
+function forced_rhs(rhs, f, drag::Real)
+    isnothing(f) && drag == 0 && return rhs
+    return function (u, p, t)
+        du = rhs(u, p, t)
+        if !isnothing(f)
+            du = du .+ f
+        end
+        if drag != 0
+            du = du .- eltype(u)(drag) .* u
+        end
+        return du
+    end
+end
+
+
 # ---------------------------------------------------------------------------
 # Public API (called from Python via juliacall)
 # ---------------------------------------------------------------------------
@@ -226,12 +263,24 @@ used for the VJP primal. `RKMethods.RK44` applies the same per-stage BC +
 projection as the projected RHS in `ns_forward_*`, so the two paths agree to
 round-off.
 """
-function ns_apply(v0_np, nu::Float64, dt::Float64, steps::Int, n::Int, L::Float64)
+function ns_apply(v0_np, nu::Float64, dt::Float64, steps::Int, n::Int, L::Float64;
+                  force = nothing, drag::Float64 = 0.0)
     T = Float32
     v0   = to_device(T.(v0_np))
     ndim = size(v0, ndims(v0))  # last dim: 2 or 3
     setup, psolver = get_setup_and_psolver(n, L, ndim)
 
+    if !isnothing(force) || drag != 0
+        # A body force or drag runs the same RK4 loop as the VJP, so that both
+        # integrate the same right-hand side. 2-D only.
+        rhs = forced_rhs(
+            create_right_hand_side(setup, psolver),
+            staggered_force(force, setup, psolver, n),
+            drag,
+        )
+        return Array(T.(ns_forward_2d(v0, rhs, setup, psolver, T(nu), T(dt), steps, n)))
+    end
+
     u = initial_state(v0, setup, psolver, n)
 
     state, _ = solve_unsteady(;
@@ -265,13 +314,18 @@ exact), which removed the need for the finite-difference fallback (and its
 two extra forward rollouts per VJP call) that earlier versions required.
 """
 function ns_vjp(v0_np, cotangent_np, nu::Float64, dt::Float64,
-                steps::Int, n::Int, L::Float64)
+                steps::Int, n::Int, L::Float64;
+                force = nothing, drag::Float64 = 0.0)
     v0   = to_device(Float32.(v0_np))
     cot  = to_device(Float32.(cotangent_np))
     ndim = size(v0, ndims(v0))
 
     setup, psolver = get_setup_and_psolver(n, L, ndim)
-    rhs  = create_right_hand_side(setup, psolver)
+    rhs  = forced_rhs(
+        create_right_hand_side(setup, psolver),
+        staggered_force(force, setup, psolver, n),
+        drag,
+    )
 
     if ndim == 2
         fwd = (v, dt_, nu_) -> ns_forward_2d(v, rhs, setup, psolver, nu_, dt_, steps, n)
diff --git a/mosaic/tesseracts/navier-stokes-grid/ins-jl/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/ins-jl/tesseract_api.py
index f4727a85..cf42c981 100644
--- a/mosaic/tesseracts/navier-stokes-grid/ins-jl/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/ins-jl/tesseract_api.py
@@ -57,12 +57,26 @@ def _to_numpy(jl_arr: Any) -> np.ndarray:
     return np.asarray(jl_arr).copy()
 
 
+def _body_force(inputs: InputSchema) -> dict[str, Any]:
+    """Keyword arguments for the body force and drag, empty when there are none.
+
+    Both are 2-D only (the canonical schema rejects them in 3-D).
+    """
+    kwargs: dict[str, Any] = {}
+    if inputs.forcing is not None:
+        kwargs["force"] = _to_julia(np.asarray(inputs.forcing)[:, :, 0, :])
+    if inputs.linear_drag is not None:
+        kwargs["drag"] = float(inputs.linear_drag[0])
+    return kwargs
+
+
 def _run(
     v0: np.ndarray,
     viscosity: float,
     dt: float,
     steps: int,
     domain_extent: float,
+    **body_force: Any,
 ) -> np.ndarray:
     """Forward pass.
 
@@ -84,6 +98,7 @@ def _run(
                 int(steps),
                 int(nx),
                 float(domain_extent),
+                **body_force,
             )
         )  # (nx, ny, 2)
         return result[:, :, np.newaxis, :]  # (nx, ny, 1, 2)
@@ -108,6 +123,7 @@ def _vjp(
     dt: float,
     steps: int,
     domain_extent: float,
+    **body_force: Any,
 ) -> dict:
     """VJP: cotangent w.r.t. result -> gradients w.r.t. v0, viscosity, dt, domain_extent."""
     nx = v0.shape[0]
@@ -125,6 +141,7 @@ def _vjp(
             int(steps),
             int(nx),
             float(domain_extent),
+            **body_force,
         )
         grad_v0_2d = _to_numpy(result[0])
         grad_nu = float(result[1])
@@ -171,6 +188,7 @@ def apply(inputs: InputSchema) -> OutputSchema:
         float(inputs.dt[0]),
         inputs.steps,
         inputs.domain_extent,
+        **_body_force(inputs),
     )
     return {"result": result.astype(np.float32)}
 
@@ -193,6 +211,7 @@ def vector_jacobian_product(
         float(inputs.dt[0]),
         inputs.steps,
         inputs.domain_extent,
+        **_body_force(inputs),
     )
     return {key: grads[key] for key in vjp_inputs if key in grads}
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/jax-cfd/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/jax-cfd/tesseract_api.py
index f15e637a..253838e1 100644
--- a/mosaic/tesseracts/navier-stokes-grid/jax-cfd/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/jax-cfd/tesseract_api.py
@@ -51,17 +51,25 @@ def cfd_fwd(
     steps: int,
     inner_steps: int,
     domain_extent: float,
+    forcing: jnp.ndarray | None = None,
+    linear_drag: float | None = None,
     **_kwargs: Any,
 ) -> jax.Array:
     """Compute the final velocity field of the incompressible Navier-Stokes equations.
 
     Supports 2D and 3D on a fully periodic domain, inferred from
     v0.shape[-1]. The result has the same shape as v0.
+
+    ``forcing`` and ``linear_drag`` enter the explicit right-hand side that
+    every RK4 stage evaluates, the force averaged onto the faces like v0 and
+    the drag ``-linear_drag * u`` on the current face velocity.
     """
     ndim = v0.shape[-1]  # 2 for 2D, 3 for 3D
     if ndim == 2:
         # squeeze the dummy nz=1 axis: (nx, ny, 1, 2) → (nx, ny, 2)
         v0 = v0[:, :, 0, :]
+        if forcing is not None:
+            forcing = forcing[:, :, 0, :]
     domain_sizes = (domain_extent,) * ndim
 
     bc = cfd.boundaries.periodic_boundary_conditions(ndim)
@@ -92,6 +100,19 @@ def make_grid_vars(v_2d: jnp.ndarray) -> Any:
     def interp_to_coloc(v_stag: Any, own_axis: int) -> Any:
         return 0.5 * (v_stag + jnp.roll(v_stag, 1, axis=own_axis))
 
+    def body_force(v: Any) -> tuple:
+        # jax-cfd's own sum_forcings calls jax.tree_map, which JAX 0.10
+        # removed, so the force and the drag are summed here.
+        terms = []
+        for i, u in enumerate(v):
+            term = 0.0
+            if forcing is not None:
+                term = interp_to_face(forcing[..., i], own_axis=i)
+            if linear_drag is not None:
+                term = term - linear_drag * u.data
+            terms.append(cfd.grids.GridArray(term, u.offset, grid))
+        return tuple(terms)
+
     # --- NS step function ------------------------------------------------
     # The library's default (van Leer + forward Euler) uses a Lax-Wendroff
     # flux whose time correction assumes pure advection. In incompressible
@@ -106,6 +127,7 @@ def interp_to_coloc(v_stag: Any, own_axis: int) -> Any:
             grid=grid,
             convect=cfd.advection.convect_linear,
             time_stepper=cfd.time_stepping.classic_rk4,
+            forcing=(None if forcing is None and linear_drag is None else body_force),
         ),
         steps=inner_steps,
     )
@@ -135,8 +157,8 @@ def apply_jit(inputs: dict) -> dict:
 
 def _unpack_scalars(d: dict) -> dict:
     """Extract Python floats from 1-element arrays for JIT-static scalar params."""
-    for key in ("density", "viscosity", "dt"):
-        if key in d:
+    for key in ("density", "viscosity", "dt", "linear_drag"):
+        if d.get(key) is not None:
             d[key] = float(d[key][0])
     return d
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/openfoam/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/openfoam/tesseract_api.py
index 6d017940..15ab9ff5 100644
--- a/mosaic/tesseracts/navier-stokes-grid/openfoam/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/openfoam/tesseract_api.py
@@ -989,6 +989,13 @@ def apply(inputs: InputSchema) -> OutputSchema:
 
     With an obstacle, also return the final pressure and the force history.
     """
+    if inputs.forcing is not None or inputs.linear_drag is not None:
+        # Supporting a source term needs the foamRun incompressibleFluid
+        # solver with a coded fvModel in place of icoFoam.
+        raise NotImplementedError(
+            "openfoam runs icoFoam, which has no momentum source term, so it "
+            "cannot apply forcing or linear_drag"
+        )
     if inputs.obstacle is not None:
         return _run_channel(inputs)
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/phiflow/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/phiflow/tesseract_api.py
index bd529483..0ac418ad 100644
--- a/mosaic/tesseracts/navier-stokes-grid/phiflow/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/phiflow/tesseract_api.py
@@ -95,6 +95,8 @@ def periodic_fwd(
     dt: float,
     steps: int,
     domain_extent: float,
+    forcing: jnp.ndarray | None = None,
+    linear_drag: float | None = None,
 ) -> jnp.ndarray:
     """Run a 2D or 3D incompressible Navier-Stokes simulation in a periodic box.
 
@@ -109,6 +111,9 @@ def periodic_fwd(
         dt: Timestep size.
         steps: Number of simulation steps.
         domain_extent: Side length of the periodic square/cubic domain.
+        forcing: Body force at the cell centers, shape of v0. It is resampled
+            to the faces like v0 and added to the momentum rate.
+        linear_drag: Coefficient α of a drag -αu in the momentum rate.
 
     Returns:
         Final velocity at the cell centers, same shape as v0.
@@ -117,6 +122,8 @@ def periodic_fwd(
     if ndim == 2:
         # squeeze the dummy nz=1 axis: (nx, ny, 1, 2) → (nx, ny, 2)
         v0 = v0[:, :, 0, :]
+        if forcing is not None:
+            forcing = forcing[:, :, 0, :]
 
     if ndim == 2:
         nx, ny = v0.shape[:2]
@@ -150,10 +157,37 @@ def staggered_to_faces(vel: StaggeredGrid) -> jnp.ndarray:
     # at step counts ≥ 20 (default solver produces NaN in 3D; 2D unaffected).
     _cg_solve = math.Solve("CG", _CG_TOL, _CG_TOL)
 
+    def to_faces(v: jnp.ndarray) -> jnp.ndarray:
+        """Collocated (*spatial, ndim) → face-centred (ndim, *spatial).
+
+        PhiFlow's native CenteredGrid→StaggeredGrid resampling.
+        """
+        v_arr = jnp.moveaxis(v, -1, 0)  # (*spatial, ndim) → (ndim, *spatial)
+        channel = math.stack(
+            [math.tensor(v_arr[i], math.spatial(spatial_str)) for i in range(ndim)],
+            dim=math.channel("vector"),
+        )
+        staggered = StaggeredGrid(
+            CenteredGrid(channel, ext, bounds=bounds, **grid_kwargs),
+            ext,
+            bounds=bounds,
+            **grid_kwargs,
+        )
+        return staggered_to_faces(staggered)
+
+    force_faces = None if forcing is None else faces_to_staggered(to_faces(forcing))
+
+    # The force and the drag are part of the momentum rate, so every RK4 stage
+    # sees them, the drag at that stage's velocity.
     def momentum(vel: StaggeredGrid) -> StaggeredGrid:
-        return advect.finite_difference(vel, vel, order=2) + diffuse.finite_difference(
+        rate = advect.finite_difference(vel, vel, order=2) + diffuse.finite_difference(
             vel, viscosity, order=2
         )
+        if force_faces is not None:
+            rate = rate + force_faces
+        if linear_drag is not None:
+            rate = rate - linear_drag * vel
+        return rate
 
     # phiflow's documented pairing for accurate periodic flow: finite
     # differences under incompressible_rk4 instead of operator splitting
@@ -176,21 +210,9 @@ def step(
         )
         return (staggered_to_faces(vel), p.values.native(spatial_str)), None
 
-    # Convert IC from collocated to staggered face values once using PhiFlow's
-    # native CenteredGrid→StaggeredGrid resampling.  The scan carry then stays
+    # Convert the IC to staggered face values once. The scan carry then stays
     # in face-value space so no smoothing accumulates across steps.
-    v_arr = jnp.moveaxis(v0, -1, 0)  # (*spatial, ndim) → (ndim, *spatial)
-    channel = math.stack(
-        [math.tensor(v_arr[i], math.spatial(spatial_str)) for i in range(ndim)],
-        dim=math.channel("vector"),
-    )
-    vel_init = StaggeredGrid(
-        CenteredGrid(channel, ext, bounds=bounds, **grid_kwargs),
-        ext,
-        bounds=bounds,
-        **grid_kwargs,
-    )
-    face_init = staggered_to_faces(vel_init)
+    face_init = to_faces(v0)
 
     # fluid.make_incompressible's CG while_loop is wrapped with custom_vjp,
     # so the backward pass does an adjoint solve rather than differentiating
@@ -399,17 +421,19 @@ def apply_jit(inputs: dict) -> dict:
         inputs["dt"],
         inputs["steps"],
         inputs["domain_extent"],
+        inputs["forcing"],
+        inputs["linear_drag"],
     )
     return {"result": result}
 
 
-_SCALAR_KEYS = ("viscosity", "dt")
+_SCALAR_KEYS = ("viscosity", "dt", "linear_drag")
 
 
 def _unpack_scalars(d: dict) -> dict:
     """Extract Python floats from 1-element arrays for JIT-static scalar params."""
     for key in _SCALAR_KEYS:
-        if key in d:
+        if d.get(key) is not None:
             d[key] = float(d[key][0])
     return d
 
diff --git a/mosaic/tesseracts/navier-stokes-grid/pict/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/pict/tesseract_api.py
index 17e73c28..f596ea12 100644
--- a/mosaic/tesseracts/navier-stokes-grid/pict/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/pict/tesseract_api.py
@@ -754,12 +754,32 @@ def _run_periodic(
     ndim: int,
     dtype: torch.dtype,
     differentiable: bool,
+    force: torch.Tensor | None = None,
+    drag: float | None = None,
 ) -> tuple[torch.Tensor, PISOtorch.Domain]:
-    """Run PICT on the periodic box and return the final velocity in PICT layout."""
+    """Run PICT on the periodic box and return the final velocity in PICT layout.
+
+    ``force`` (PICT layout) and ``drag`` enter as PICT's velocity source,
+    ``force - drag * u`` from the velocity at the start of each step, so the
+    drag is explicit. Both are in PICT's unit-cell units.
+    """
     domain = _make_periodic_domain(viscosity_val, N, ndim, dtype)
     domain.PrepareSolve()  # allocate internal structures first
-    domain.getBlock(0).setVelocity(v0_tensor)
+    block = domain.getBlock(0)
+    block.setVelocity(v0_tensor)
     domain.UpdateDomainData()
+    prep_fn = None
+    if force is not None or drag is not None:
+
+        def _pre_step(domain: PISOtorch.Domain, **_kw: Any) -> None:
+            # Set every step, from the current velocity, so the source stays
+            # in the graph that PISOtorch_diff records.
+            source = torch.zeros_like(block.velocity) if force is None else force
+            if drag is not None:
+                source = source - drag * block.velocity
+            block.setVelocitySource(source.contiguous())
+
+        prep_fn = {"PRE": _pre_step}
     # PICT's float32 default tolerance for the advection and pressure solves is
     # 1e-5. Its residual error in the forward and adjoint solves puts a floor of
     # about 1e-3 under the VJP's agreement with finite differences on the
@@ -769,11 +789,12 @@ def _run_periodic(
         dt_val,
         steps,
         differentiable,
+        prep_fn=prep_fn,
         non_orthogonal=False,
         advection_tol=_PERIODIC_SOLVE_TOL,
         pressure_tol=_PERIODIC_SOLVE_TOL,
     )
-    return domain.getBlock(0).velocity, domain
+    return block.velocity, domain
 
 
 def _run_channel(
@@ -903,11 +924,35 @@ def apply(inputs: InputSchema) -> OutputSchema:
     nu_pict_t = torch.tensor([nu_pict], dtype=dtype, device=torch.device("cpu"))
     v0_t = _v0_to_pict(v0_np, _DEVICE, dtype, requires_grad=False)
     result_t, _ = _run_periodic(
-        v0_t, nu_pict_t, dt_pict, inputs.steps, N, ndim, dtype, differentiable=False
+        v0_t,
+        nu_pict_t,
+        dt_pict,
+        inputs.steps,
+        N,
+        ndim,
+        dtype,
+        differentiable=False,
+        **_periodic_source(inputs, dtype),
     )
     return {"result": _pict_to_v0(result_t, N, ndim).detach().cpu().numpy()}
 
 
+def _periodic_source(inputs: InputSchema, dtype: torch.dtype) -> dict[str, Any]:
+    """Body force and drag in PICT's unit-cell units.
+
+    With x_pict = x · N/L and the velocity unchanged, time scales like length,
+    so an acceleration and a drag rate both divide by N/L.
+    """
+    phys_scale, _, _ = _periodic_scaling(inputs)
+    force = None
+    if inputs.forcing is not None:
+        force = _v0_to_pict(np.asarray(inputs.forcing) / phys_scale, _DEVICE, dtype)
+    drag = None
+    if inputs.linear_drag is not None:
+        drag = float(inputs.linear_drag[0]) / phys_scale
+    return {"force": force, "drag": drag}
+
+
 def vector_jacobian_product(
     inputs: InputSchema,
     vjp_inputs: set[str],
@@ -959,7 +1004,15 @@ def vector_jacobian_product(
         nu_t = torch.tensor([nu_pict], dtype=dtype, requires_grad=want["viscosity"])
         profile_t = None
         result_t, domain = _run_periodic(
-            v0_t, nu_t, dt_pict, inputs.steps, N, ndim, dtype, differentiable=True
+            v0_t,
+            nu_t,
+            dt_pict,
+            inputs.steps,
+            N,
+            ndim,
+            dtype,
+            differentiable=True,
+            **_periodic_source(inputs, dtype),
         )
         if "result" in vjp_outputs and "result" in cotangent_vector:
             outputs.append(result_t)
diff --git a/mosaic/tesseracts/navier-stokes-grid/warp-ns/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/warp-ns/tesseract_api.py
index a7c5f5b7..8a935942 100644
--- a/mosaic/tesseracts/navier-stokes-grid/warp-ns/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/warp-ns/tesseract_api.py
@@ -48,12 +48,17 @@ def tentative_vel_2d_kernel(
     inv_2h: float,
     inv_h2: float,
     nu: wp.array(dtype=wp.float32),
+    fx: wp.array2d(dtype=wp.float32),
+    fy: wp.array2d(dtype=wp.float32),
+    drag: wp.array(dtype=wp.float32),
 ) -> None:
-    """2-D tentative velocity: u* = u + dt·(-u·∇u + ν∇²u).
+    """2-D tentative velocity: u* = u + dt·(-u·∇u + ν∇²u + f - αu).
 
     ``dt`` and ``nu`` are 1-element arrays (rather than plain floats) so that
     Warp's source-to-source autodiff tracks gradients w.r.t. them: Warp only
-    differentiates through array-typed kernel arguments.
+    differentiates through array-typed kernel arguments. The body force
+    ``(fx, fy)`` and the drag coefficient ``drag`` (α) are zero when the inputs
+    carry none.
 
     [2D-only function]
     """
@@ -76,7 +81,7 @@ def tentative_vel_2d_kernel(
         ui * (ux[ip1, j] - ux[im1, j]) * inv_2h
         + vi * (ux[i, jp1] - ux[i, jm1]) * inv_2h
     )
-    ux_star[i, j] = ui + dt_ * (-adv_ux + nu_ * lap_ux)
+    ux_star[i, j] = ui + dt_ * (-adv_ux + nu_ * lap_ux + fx[i, j] - drag[0] * ui)
 
     # uy component
     lap_uy = (uy[im1, j] + uy[ip1, j] + uy[i, jm1] + uy[i, jp1] - 4.0 * vi) * inv_h2
@@ -84,7 +89,7 @@ def tentative_vel_2d_kernel(
         ui * (uy[ip1, j] - uy[im1, j]) * inv_2h
         + vi * (uy[i, jp1] - uy[i, jm1]) * inv_2h
     )
-    uy_star[i, j] = vi + dt_ * (-adv_uy + nu_ * lap_uy)
+    uy_star[i, j] = vi + dt_ * (-adv_uy + nu_ * lap_uy + fy[i, j] - drag[0] * vi)
 
 
 @wp.kernel
@@ -1204,6 +1209,26 @@ def _pass_bufs():
     }
 
 
+def _body_force_2d(
+    forcing_np: np.ndarray | None, linear_drag: float | None, n: int, device: str
+) -> tuple[wp.array, wp.array, wp.array]:
+    """Force components and drag coefficient for ``tentative_vel_2d_kernel``.
+
+    The force is collocated like the velocity, so it is used as given. Absent
+    inputs become zeros, which leave the update unchanged.
+    """
+    if forcing_np is None:
+        forcing_np = np.zeros((n, n, 1, 2), dtype=np.float32)
+    fx = wp.array(forcing_np[:, :, 0, 0], dtype=wp.float32, device=device)
+    fy = wp.array(forcing_np[:, :, 0, 1], dtype=wp.float32, device=device)
+    drag = wp.array(
+        np.array([linear_drag or 0.0], dtype=np.float32),
+        dtype=wp.float32,
+        device=device,
+    )
+    return fx, fy, drag
+
+
 def ns2d_solve_tape(
     v0_np: np.ndarray,
     viscosity: float,
@@ -1212,6 +1237,8 @@ def ns2d_solve_tape(
     domain_extent: float,
     device: str = "cpu",
     track_scalar_grads: bool = False,
+    forcing_np: np.ndarray | None = None,
+    linear_drag: float | None = None,
 ) -> tuple[wp.Tape, wp.array, wp.array, wp.array, wp.array, wp.array, wp.array]:
     """Run periodic 2-D incompressible NS by Chorin projection under a tape.
 
@@ -1282,6 +1309,7 @@ def ns2d_solve_tape(
         requires_grad=track_scalar_grads,
         device=device,
     )
+    fx_wp, fy_wp, drag_wp = _body_force_2d(forcing_np, linear_drag, n, device)
 
     tape = wp.Tape()
     with tape:
@@ -1297,7 +1325,19 @@ def ns2d_solve_tape(
             wp.launch(
                 tentative_vel_2d_kernel,
                 dim=(n, n),
-                inputs=[cur_ux, cur_uy, ux_star, uy_star, dt_wp, inv_2h, inv_h2, nu_wp],
+                inputs=[
+                    cur_ux,
+                    cur_uy,
+                    ux_star,
+                    uy_star,
+                    dt_wp,
+                    inv_2h,
+                    inv_h2,
+                    nu_wp,
+                    fx_wp,
+                    fy_wp,
+                    drag_wp,
+                ],
                 block_dim=_bd_2d,
                 device=device,
             )
@@ -1366,6 +1406,8 @@ def ns2d_solve_forward(
     steps: int,
     domain_extent: float,
     device: str = "cpu",
+    forcing_np: np.ndarray | None = None,
+    linear_drag: float | None = None,
 ) -> np.ndarray:
     """Forward-only periodic 2-D incompressible NS by Chorin projection, no gradients.
 
@@ -1399,6 +1441,7 @@ def ns2d_solve_forward(
         np.array([viscosity], dtype=np.float32), dtype=wp.float32, device=device
     )
     dt_wp = wp.array(np.array([dt], dtype=np.float32), dtype=wp.float32, device=device)
+    fx_wp, fy_wp, drag_wp = _body_force_2d(forcing_np, linear_drag, n, device)
 
     ux_star = wp.zeros((n, n), dtype=wp.float32, device=device)
     uy_star = wp.zeros((n, n), dtype=wp.float32, device=device)
@@ -1419,6 +1462,9 @@ def ns2d_solve_forward(
                 inv_2h,
                 inv_h2,
                 nu_wp,
+                fx_wp,
+                fy_wp,
+                drag_wp,
             ],
             block_dim=_bd_2d,
             device=device,
@@ -2557,11 +2603,30 @@ def apply(inputs: InputSchema) -> OutputSchema:
         )
     else:
         result = ns2d_solve_forward(
-            v0, nu, dt, inputs.steps, inputs.domain_extent, device=device
+            v0,
+            nu,
+            dt,
+            inputs.steps,
+            inputs.domain_extent,
+            device=device,
+            **_body_force_args(inputs),
         )
     return OutputSchema(result=result)
 
 
+def _body_force_args(inputs: InputSchema) -> dict[str, Any]:
+    return {
+        "forcing_np": (
+            None
+            if inputs.forcing is None
+            else np.asarray(inputs.forcing, dtype=np.float32)
+        ),
+        "linear_drag": (
+            None if inputs.linear_drag is None else float(inputs.linear_drag[0])
+        ),
+    }
+
+
 def channel_vjp(
     inputs: InputSchema,
     vjp_inputs: set[str],
@@ -2672,6 +2737,7 @@ def vector_jacobian_product(
             inputs.domain_extent,
             device=device,
             track_scalar_grads=track_scalar_grads,
+            **_body_force_args(inputs),
         )
         cot_result = np.asarray(
             cotangent_vector.get("result", np.zeros_like(v0)), dtype=np.float32
diff --git a/mosaic/tesseracts/navier-stokes-grid/xlb/tesseract_api.py b/mosaic/tesseracts/navier-stokes-grid/xlb/tesseract_api.py
index a5f14e08..4232adf9 100644
--- a/mosaic/tesseracts/navier-stokes-grid/xlb/tesseract_api.py
+++ b/mosaic/tesseracts/navier-stokes-grid/xlb/tesseract_api.py
@@ -196,7 +196,7 @@ def _sub_steps(inputs: dict) -> int:
     """
     if inputs.get("obstacle") is None:
         v0 = np.asarray(inputs["v0"], dtype=np.float64)
-        u_max = float(np.sqrt((v0**2).sum(axis=-1)).max())
+        u_max = max(float(np.sqrt((v0**2).sum(axis=-1)).max()), _forced_speed(inputs))
         ma_target = _MA_TARGET * min(1.0, _MA_TARGET_N / v0.shape[0])
     else:
         u_max = float(np.max(np.abs(np.asarray(inputs["inflow_profile"]))))
@@ -205,6 +205,25 @@ def _sub_steps(inputs: dict) -> int:
     return max(1, math.ceil(ma / ma_target))
 
 
+def _forced_speed(inputs: dict) -> float:
+    """Speed bound for a body force, independent of the flow it drives.
+
+    A steady force F against viscosity and drag drives a laminar flow of at most
+    ‖F‖∞ / (ν κ_min² + α), the amplitude of its largest-scale mode, with
+    κ_min = 2π/L. Sizing the lattice Mach number by it keeps a flow that the
+    force spins up from rest within the sub-step count chosen at the start.
+    """
+    forcing = inputs.get("forcing")
+    if forcing is None:
+        return 0.0
+    kappa_min = 2.0 * math.pi / float(inputs["domain_extent"])
+    damping = inputs["viscosity"] * kappa_min**2 + (inputs.get("linear_drag") or 0.0)
+    return (
+        float(np.sqrt((np.asarray(forcing, np.float64) ** 2).sum(axis=-1)).max())
+        / damping
+    )
+
+
 # mosaic-intervention: consistent-initial-populations
 def _periodic_populations(u: jnp.ndarray, omega: Any, ops: dict) -> jnp.ndarray:
     """Populations consistent with a divergence-free periodic velocity field.
@@ -388,6 +407,8 @@ def xlb_fwd(
     dx: float,
     obstacle: dict | None = None,
     inflow_profile: jnp.ndarray | None = None,
+    forcing: jnp.ndarray | None = None,
+    linear_drag: Any = None,
     _use_f64: bool = False,
     _sub_k: int = 1,
     **_: Any,
@@ -414,6 +435,8 @@ def xlb_fwd(
         obstacle:      Cylinder (``center``, ``radius``, physical units) in the
                        channel, or None for a periodic box.
         inflow_profile: Channel inlet u_x at the cell-center heights, shape (ny,).
+        forcing:       Periodic body force (an acceleration), shape of v0.
+        linear_drag:   Periodic drag coefficient α of a force -αu.
         _use_f64:      If True, run the entire LBM computation in float64 for accurate
                        gradients. Output is cast back to float32 before returning.
                        This prevents float32 cancellation errors in the VJP/JVP paths,
@@ -451,7 +474,41 @@ def xlb_fwd(
         u0 = jnp.moveaxis(v0, -1, 0).astype(fdtype) * scale_eff
 
     out = {}
-    if obstacle is None:
+    if obstacle is None and (forcing is not None or linear_drag is not None):
+        # Body force by the exact difference method (Kupershtokh 2004), which
+        # adds feq(ρ, u + Δu) - feq(ρ, u) after the collision. With the force
+        # a - αu per unit mass, the kick per LBM step is Δu = a - αu* at the
+        # physical velocity u* = u + Δu/2, solved for u* since the drag depends
+        # on it. Lattice units: a dt²/dx and α dt.
+        accel = 0.0
+        if forcing is not None:
+            accel = jnp.moveaxis(forcing[:, :, 0, :], -1, 0).astype(fdtype)
+            accel = accel * scale_eff * dt_eff
+        drag = 0.0 if linear_drag is None else linear_drag * dt_eff
+
+        def kick(u_phys: jnp.ndarray) -> jnp.ndarray:
+            return accel - drag * u_phys
+
+        # mosaic-intervention: consistent-initial-populations
+        # Post-collision populations carry u* + Δu/2: half the kick of their
+        # step is already in them.
+        f0 = _periodic_populations(u0 + 0.5 * kick(u0), omega, ops)
+        # mosaic-intervention-end: consistent-initial-populations
+
+        # mosaic-intervention: manual-lbm-loop
+        def body(f: jnp.ndarray, _: Any) -> tuple[jnp.ndarray, None]:
+            f_s = xlb_stream(f)
+            rho_s, u_s = xlb_macro(f_s)
+            feq = xlb_eq(rho_s, u_s)
+            f_post = xlb_bgk(f_s, feq, rho_s, u_s, omega)
+            u_phys = (u_s + 0.5 * accel) / (1.0 + 0.5 * drag)
+            return f_post + xlb_eq(rho_s, u_s + kick(u_phys)) - feq, None
+
+        f_final, _ = _checkpointed_scan(body, f0, steps_eff)
+        # mosaic-intervention-end: manual-lbm-loop
+        _rho_f, u_post = xlb_macro(f_final)
+        u_out = (u_post - 0.5 * accel) / (1.0 - 0.5 * drag)
+    elif obstacle is None:
         # mosaic-intervention: consistent-initial-populations
         f0 = _periodic_populations(u0, omega, ops)
         # mosaic-intervention-end: consistent-initial-populations
@@ -507,8 +564,9 @@ def apply_jit(inputs: dict) -> dict:
 def _concrete_inputs(inputs: InputSchema) -> dict:
     """Inputs as a dict with Python-float viscosity and dt, and the grid spacing."""
     d = inputs.model_dump()
-    for key in ("viscosity", "dt"):
-        d[key] = float(d[key][0])
+    for key in ("viscosity", "dt", "linear_drag"):
+        if d[key] is not None:
+            d[key] = float(d[key][0])
     d["dx"] = inputs.dx
     return d
 
@@ -575,11 +633,14 @@ def _diff_apply(sub_k: int, dx: float) -> Any:
 
     def fn(inputs: dict) -> dict:
         fwd_kwargs = dict(inputs)
-        for k in ("v0", "inflow_profile"):
+        for k in ("v0", "inflow_profile", "forcing"):
             if fwd_kwargs.get(k) is not None:
                 fwd_kwargs[k] = jnp.asarray(fwd_kwargs[k], dtype=_diff_dtype)
-        for k in ("viscosity", "dt"):
-            fwd_kwargs[k] = jnp.asarray(fwd_kwargs[k], dtype=_diff_dtype).reshape(())
+        for k in ("viscosity", "dt", "linear_drag"):
+            if fwd_kwargs.get(k) is not None:
+                fwd_kwargs[k] = jnp.asarray(fwd_kwargs[k], dtype=_diff_dtype).reshape(
+                    ()
+                )
         fwd_kwargs["dx"] = dx
         return xlb_fwd(
             _use_f64=_USE_F64_DIFF,
diff --git a/tests/test_ns_grid_forcing.py b/tests/test_ns_grid_forcing.py
new file mode 100644
index 00000000..0b7d6102
--- /dev/null
+++ b/tests/test_ns_grid_forcing.py
@@ -0,0 +1,181 @@
+# Copyright 2026 Pasteur Labs. All Rights Reserved.
+# SPDX-License-Identifier: Apache-2.0
+
+"""Body force and linear drag in the periodic ns-grid solvers.
+
+Started from rest, the Kolmogorov force (A sin(κy), 0) keeps the flow a shear
+flow u(y, t). Its advection term vanishes, so the Navier–Stokes equations with
+drag α have the exact solution
+
+    u = (A/λ)(1 − e^{−λt}) sin(κy),  v = 0,  λ = νκ² + α.
+
+One run against it checks the amplitude, sign and wavenumber of the force, its
+placement on staggered grids, the sign of the drag and any unit conversion.
+
+Every test runs against the solver's ``:latest`` image built from this
+checkout and is skipped where that image or Docker is missing. They are never
+built here.
+"""
+
+from __future__ import annotations
+
+import shutil
+
+import numpy as np
+import pytest
+
+from .helpers import image_available
+
+_L = 2 * np.pi
+
+# Solver image and the extra inputs the harness passes it.
+_SOLVERS = {
+    "jax-cfd": (
+        "jax_cfd_navier_stokes_grid:latest",
+        {"density": np.array([1.0], np.float32), "inner_steps": 1},
+    ),
+    "phiflow": ("phiflow_navier_stokes_grid:latest", {}),
+    "warp-ns": ("warp_ns_navier_stokes_grid:latest", {}),
+    "xlb": ("xlb_navier_stokes_grid:latest", {}),
+    "ins-jl": ("ins_navier_stokes_grid:latest", {}),
+    "pict": ("pict_navier_stokes_grid:latest", {}),
+    "exponax": ("exponax_navier_stokes_grid:latest", {}),
+}
+
+_OPENFOAM = "openfoam_navier_stokes_grid:latest"
+
+
+@pytest.fixture(scope="module", params=sorted(_SOLVERS))
+def solver(request):
+    image, extra = _SOLVERS[request.param]
+    if not image_available(image):
+        pytest.skip(f"{image} not built")
+    from tesseract_core import Tesseract
+
+    # CI builds some of these images on runners without GPUs, where the JAX
+    # solvers fall back to the CPU.
+    gpus = ["all"] if shutil.which("nvidia-smi") else None
+    with Tesseract.from_image(image, gpus=gpus) as t:
+        yield request.param, t, extra
+
+
+def _kolmogorov_force(n: int, k: int, amplitude: float = 1.0) -> np.ndarray:
+    """(A sin(κy), 0) at the cell centres, shape (n, n, 1, 2)."""
+    y = (np.arange(n) + 0.5) * _L / n
+    force = np.zeros((n, n, 1, 2), np.float32)
+    force[..., 0, 0] = amplitude * np.sin(2 * np.pi * k * y / _L)[None, :]
+    return force
+
+
+def _inputs(v0, *, nu, dt, steps, extra, **more) -> dict:
+    return {
+        "v0": np.asarray(v0, np.float32),
+        "viscosity": np.array([nu], np.float32),
+        "dt": np.array([dt], np.float32),
+        "steps": steps,
+        "domain_extent": float(_L),
+        **extra,
+        **more,
+    }
+
+
+def _smooth_field(n: int, seed: int = 0, amplitude: float = 0.3) -> np.ndarray:
+    """A divergence-free field from a few low modes of a stream function."""
+    rng = np.random.default_rng(seed)
+    x = (np.arange(n) + 0.5) * _L / n
+    xx, yy = np.meshgrid(x, x, indexing="ij")
+    u = np.zeros((n, n))
+    v = np.zeros((n, n))
+    for kx, ky in ((1, 0), (0, 1), (1, 1), (2, 1), (1, 2)):
+        a, phase = rng.normal(), rng.uniform(0, 2 * np.pi)
+        # ψ = a sin(kx x + ky y + φ); u = ∂ψ/∂y, v = −∂ψ/∂x
+        c = a * np.cos(kx * xx + ky * yy + phase)
+        u += ky * c
+        v -= kx * c
+    field = np.stack([u, v], -1)[:, :, None, :]
+    return (amplitude * field / np.abs(field).max()).astype(np.float32)
+
+
+def test_laminar_kolmogorov_transient(solver):
+    """Within 5% of the exact shear flow at N = 64 after two time units."""
+    _, t, extra = solver
+    n, k, nu, alpha, dt, steps = 64, 4, 0.05, 0.1, 0.01, 200
+    force = _kolmogorov_force(n, k)
+    out = t.apply(
+        _inputs(
+            np.zeros_like(force),
+            nu=nu,
+            dt=dt,
+            steps=steps,
+            extra=extra,
+            forcing=force,
+            linear_drag=np.array([alpha], np.float32),
+        )
+    )
+    rate = nu * k**2 + alpha
+    exact = force * (1 - np.exp(-rate * dt * steps)) / rate
+    result = np.asarray(out["result"], np.float64)
+    error = np.linalg.norm(result - exact) / np.linalg.norm(exact)
+    assert error < 0.05, f"relative error {error:.3g}"
+
+
+def test_zero_force_and_drag_change_nothing(solver):
+    """A zero force and zero drag give the unforced result."""
+    _, t, extra = solver
+    v0 = _smooth_field(16)
+    common = {"nu": 0.01, "dt": 0.05, "steps": 10, "extra": extra}
+    plain = np.asarray(t.apply(_inputs(v0, **common))["result"])
+    zero = np.asarray(
+        t.apply(
+            _inputs(
+                v0,
+                **common,
+                forcing=np.zeros_like(v0),
+                linear_drag=np.zeros(1, np.float32),
+            )
+        )["result"]
+    )
+    np.testing.assert_allclose(zero, plain, rtol=0, atol=1e-6)
+
+
+def test_forced_vjp_matches_central_difference(solver):
+    """d⟨c, u(T)⟩/dv0 along a direction, with force and drag, against FD."""
+    name, t, extra = solver
+    n = 16
+    v0 = _smooth_field(n)
+    force = _kolmogorov_force(n, 2)
+    inputs = _inputs(
+        v0,
+        nu=0.02,
+        dt=0.02,
+        steps=10,
+        extra=extra,
+        forcing=force,
+        linear_drag=np.array([0.1], np.float32),
+    )
+    cotangent = _smooth_field(n, seed=1, amplitude=1.0)
+    direction = _smooth_field(n, seed=2, amplitude=1.0)
+
+    def objective(v):
+        out = t.apply({**inputs, "v0": v.astype(np.float32)})["result"]
+        return float(np.sum(np.asarray(out, np.float64) * cotangent))
+
+    grad = t.vector_jacobian_product(inputs, ["v0"], ["result"], {"result": cotangent})[
+        "v0"
+    ]
+    eps = 1e-2
+    fd = (objective(v0 + eps * direction) - objective(v0 - eps * direction)) / (2 * eps)
+    ad = float(np.sum(np.asarray(grad, np.float64) * direction))
+    assert ad == pytest.approx(fd, rel=1e-2), f"{name}: AD {ad:.6g} vs FD {fd:.6g}"
+
+
+@pytest.mark.skipif(not image_available(_OPENFOAM), reason=f"{_OPENFOAM} not built")
+def test_openfoam_rejects_forcing():
+    from tesseract_core import Tesseract
+
+    v0 = np.zeros((8, 8, 1, 2), np.float32)
+    with (
+        Tesseract.from_image(_OPENFOAM) as t,
+        pytest.raises(Exception, match="icoFoam"),
+    ):
+        t.apply(_inputs(v0, nu=0.01, dt=0.01, steps=1, extra={}, forcing=v0))
diff --git a/tests/test_schemas.py b/tests/test_schemas.py
index a686ffbf..ec9348a9 100644
--- a/tests/test_schemas.py
+++ b/tests/test_schemas.py
@@ -212,3 +212,50 @@ def test_ns_grid_solid_cells_are_the_cells_centered_inside_the_cylinder():
     expected = np.zeros((4, 5), dtype=bool)
     expected[1:3, 1:3] = True
     np.testing.assert_array_equal(solid, expected)
+
+
+def _forcing_case(case):
+    """Inputs that put a body force or drag where it is not defined."""
+    import numpy as np
+
+    if case == "force_in_channel":
+        inputs = _channel_inputs()
+        return {**inputs, "forcing": np.zeros_like(inputs["v0"])}
+    if case == "drag_in_channel":
+        return _channel_inputs(linear_drag=np.array([0.1], dtype=np.float32))
+    if case == "force_in_3d":
+        v0 = np.zeros((8, 8, 8, 3), dtype=np.float32)
+        return {"v0": v0, "forcing": np.zeros_like(v0)}
+    v0 = np.zeros((16, 16, 1, 2), dtype=np.float32)
+    return {"v0": v0, "forcing": np.zeros((8, 8, 1, 2), dtype=np.float32)}
+
+
+@needs_runtime
+@pytest.mark.parametrize(
+    "case, match",
+    [
+        ("force_in_channel", "periodic box without an obstacle"),
+        ("drag_in_channel", "periodic box without an obstacle"),
+        ("force_in_3d", "only defined in 2-D"),
+        ("force_off_grid", r"expected \(16, 16, 1, 2\)"),
+    ],
+)
+def test_ns_grid_forcing_is_a_2d_periodic_input(case, match):
+    from mosaic_shared.problems.navier_stokes_grid import InputSchema
+    from pydantic import ValidationError
+
+    with pytest.raises(ValidationError, match=match):
+        InputSchema(**_forcing_case(case))
+
+
+@needs_runtime
+def test_ns_grid_forcing_and_drag_are_accepted_on_a_periodic_box():
+    import numpy as np
+    from mosaic_shared.problems.navier_stokes_grid import InputSchema
+
+    v0 = np.zeros((16, 16, 1, 2), dtype=np.float32)
+    inputs = InputSchema(
+        v0=v0, forcing=np.ones_like(v0), linear_drag=np.array([0.1], np.float32)
+    )
+    assert inputs.forcing.shape == v0.shape
+    assert InputSchema(v0=v0).forcing is None