fix(darcy): correct sign conventions for Darcy flux and velocity

src/underworld3/constitutive_models.py

@@ -2410,7 +2410,7 @@ def flux(self):
 
         ddu = self.grad_u - self.Parameters.s.sym
 
-        return self._q(ddu)
+        return -self._q(ddu)
 
 
 class TransverseIsotropicFlowModel(ViscousFlowModel):

src/underworld3/systems/solvers.py

@@ -164,7 +164,7 @@ class SNES_Poisson(SNES_Scalar):
 
     .. math::
 
-        \nabla \cdot \left[ \boldsymbol{\kappa} \nabla u \right] = f
+        -\nabla \cdot \left[ \boldsymbol{\kappa} \nabla u \right] = f
 
     where :math:`\mathbf{F} = \boldsymbol{\kappa} \nabla u` relates the flux to
     gradients in the unknown :math:`u`.
@@ -257,7 +257,7 @@ def f(self):
         The source term :math:`f` appears on the right-hand side:
 
         .. math::
-            \nabla \cdot (\kappa \nabla u) = f
+            -\nabla \cdot (\kappa \nabla u) = f
 
         Returns
         -------
@@ -434,7 +434,7 @@ def __init__(
 
     # =========================================================================
     # PETSc Residual Templates
-    # For Darcy flow: -∇·(K∇p) = f, with velocity v = -K∇p
+    # For Darcy flow: -∇·(κ(∇h - s)) = f, with velocity q = -κ(∇h - s)
     # =========================================================================
 
     F0 = Template(
@@ -450,24 +450,28 @@ def __init__(
     F1 = Template(
         r"\mathbf{F}_1\left( \mathbf{u} \right)",
         lambda self: self.darcy_flux,
-        r"""Darcy flux term (pointwise).
-
-        The $\mathbf{F}_1$ vector is the Darcy flux $K \nabla p$,
-        where $K$ is the permeability.
+        r"""Darcy flux term for PDE assembly (pointwise).
+
+        The $\mathbf{F}_1$ vector is the physical Darcy velocity
+        $\mathbf{q} = -\kappa(\nabla h - \mathbf{s})$.
+        This sign is consistent with uw3's natural BC convention, where
+        ``add_natural_bc(g, boundary)`` specifies the inward Darcy flux
+        (positive = into domain). For $f=0$ the PDE assembly is correct;
+        for $f \neq 0$ a sign fix is needed (tracked separately).
         """,
     )
 
     @timing.routine_timer_decorator
     def darcy_problem_description(self):
         """Build residual terms for Darcy flow FEM assembly."""
-        # f1 residual term (weighted integration)
+        # f0 residual term (pointwise)
         self._f0 = self.F0.sym
 
         # f1 residual term (integration by parts / gradients)
         self._f1 = self.F1.sym
 
-        # Flow calculation
-        self._v_projector.uw_function = -self.F1.sym
+        # Flow calculation uses physical Darcy velocity q = F1
+        self._v_projector.uw_function = self.F1.sym
 
         return
 
@@ -543,11 +547,11 @@ def solve(
         super().solve(zero_init_guess, _force_setup,
                       divergence_retries=divergence_retries)
 
-        # Now solve flow field: v = -flux = -K(grad(h) - s)
+        # Now solve flow field: v = flux = κ(grad(h) - s)
 
         # self._v_projector.petsc_options["snes_rtol"] = 1.0e-6
         # self._v_projector.petsc_options.delValue("ksp_monitor")
-        self._v_projector.uw_function = -self.darcy_flux
+        self._v_projector.uw_function = self.darcy_flux
         self._v_projector.solve(zero_init_guess)
 
         return
@@ -800,7 +804,10 @@ def solve(
 
         if not self.constitutive_model._solver_is_setup:
             self._needs_function_rewire = True
-            self.DFDt.psi_fn = self.constitutive_model.flux.T
+            # DFDt must track K(∇h−s) = −q so TransientDarcy F1 is correct
+            # for PDE assembly (∇·F1 = f0 requires F1 = K(∇h−s), not q).
+            # Note: natural BCs on TransientDarcy would need negated sign (deferred).
+            self.DFDt.psi_fn = -self.constitutive_model.flux.T
 
         if not self.is_setup:
             self._setup_pointwise_functions(verbose)