From a807b5528804c8cb8349fa9bf890ed3ce3ffc5b2 Mon Sep 17 00:00:00 2001 From: Juan Carlos Graciosa Date: Tue, 16 Jun 2026 10:59:08 +1000 Subject: [PATCH 1/2] fix(darcy): correct sign conventions for Darcy flux and velocity --- src/underworld3/constitutive_models.py | 2 +- src/underworld3/systems/solvers.py | 31 +++++++++++++++----------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index cc75c20b..69c3415e 100644 --- a/src/underworld3/constitutive_models.py +++ b/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): diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 2e88f87e..6072e20f 100644 --- a/src/underworld3/systems/solvers.py +++ b/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( @@ -449,25 +449,28 @@ def __init__( F1 = Template( r"\mathbf{F}_1\left( \mathbf{u} \right)", - lambda self: self.darcy_flux, - r"""Darcy flux term (pointwise). + lambda self: -self.darcy_flux, + r"""Gradient flux term for PDE assembly (pointwise). - The $\mathbf{F}_1$ vector is the Darcy flux $K \nabla p$, - where $K$ is the permeability. + The $\mathbf{F}_1$ vector is $\kappa(\nabla h - \mathbf{s})$, the negative + of the physical Darcy velocity. This sign is required so that + $\nabla \cdot \mathbf{F}_1 = f_0$ recovers the correct strong form + $-\nabla \cdot [\kappa(\nabla h - \mathbf{s})] = f$. """, ) @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) + # F1 = κ(∇h − s) = −q, correct for PDE assembly self._f1 = self.F1.sym - # Flow calculation - self._v_projector.uw_function = -self.F1.sym + # Flow calculation uses physical Darcy velocity q, not F1 + self._v_projector.uw_function = self.darcy_flux return @@ -543,11 +546,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 +803,9 @@ 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 κ(∇h−s) = −q (gradient flux, not physical velocity) + # so that adams_moulton_flux() produces the correct F1 for PDE assembly. + self.DFDt.psi_fn = -self.constitutive_model.flux.T if not self.is_setup: self._setup_pointwise_functions(verbose) From c50aaac4687d03932539157664c08fdac0fe5dda Mon Sep 17 00:00:00 2001 From: Juan Carlos Graciosa Date: Tue, 16 Jun 2026 12:58:11 +1000 Subject: [PATCH 2/2] fix(darcy): restore natural BC compatibility and correct TransientDarcy DFDt --- src/underworld3/systems/solvers.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 6072e20f..da328bd4 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -449,13 +449,15 @@ def __init__( F1 = Template( r"\mathbf{F}_1\left( \mathbf{u} \right)", - lambda self: -self.darcy_flux, - r"""Gradient flux term for PDE assembly (pointwise). - - The $\mathbf{F}_1$ vector is $\kappa(\nabla h - \mathbf{s})$, the negative - of the physical Darcy velocity. This sign is required so that - $\nabla \cdot \mathbf{F}_1 = f_0$ recovers the correct strong form - $-\nabla \cdot [\kappa(\nabla h - \mathbf{s})] = f$. + lambda self: self.darcy_flux, + 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). """, ) @@ -466,11 +468,10 @@ def darcy_problem_description(self): self._f0 = self.F0.sym # f1 residual term (integration by parts / gradients) - # F1 = κ(∇h − s) = −q, correct for PDE assembly self._f1 = self.F1.sym - # Flow calculation uses physical Darcy velocity q, not F1 - self._v_projector.uw_function = self.darcy_flux + # Flow calculation uses physical Darcy velocity q = F1 + self._v_projector.uw_function = self.F1.sym return @@ -803,8 +804,9 @@ def solve( if not self.constitutive_model._solver_is_setup: self._needs_function_rewire = True - # DFDt must track κ(∇h−s) = −q (gradient flux, not physical velocity) - # so that adams_moulton_flux() produces the correct F1 for PDE assembly. + # 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: