From 12075efe024a17d107774a5538c3869857289ef7 Mon Sep 17 00:00:00 2001
From: Louis Moresi <louis.moresi@anu.edu.au>
Date: Sun, 1 Mar 2026 08:29:05 +1100
Subject: [PATCH] Fix Poisson solver natural (Neumann) boundary conditions
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Two bugs prevented natural BCs from working in the scalar Poisson solver
(SNES_Scalar), while the Stokes solver (SNES_Vector) worked correctly:

1. _setup_discretisation registered natural BCs with per-boundary labels
   (e.g., "Bottom") instead of the consolidated "UW_Boundaries" label
   that _setup_solver expects.

2. _setup_solver had a placeholder comment but no code to register the
   compiled boundary residual/jacobian functions with PETSc via
   UW_PetscDSSetBdResidual/BdJacobian/BdJacobianPreconditioner.

The fix mirrors the working SNES_Vector pattern. Verified against an
analytical solution (T = x²y) on simplex and quad meshes.

Closes #25

Underworld development team with AI support from Claude Code
---
 .../cython/petsc_generic_snes_solvers.pyx     | 45 ++++++++++++++++---
 tests/test_1000_poissonNaturalBC.py           |  1 -
 2 files changed, 40 insertions(+), 6 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index d3f92714d..0c052a22d 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -1173,10 +1173,6 @@ class SNES_Scalar(SolverBaseClass):
 
 
         self.essential_bcs = []
-        # TODO(BUG): add_natural_bc() causes PETSc error 73 ("Object in wrong state")
-        # when used with this solver. The Stokes solver's natural BCs work correctly,
-        # suggesting a setup/ordering issue specific to scalar Poisson.
-        # See planning file: underworld.md (Bugs section, 2026-01-19)
         self.natural_bcs = []
         self.bcs = self.essential_bcs
         self.boundary_conditions = False
@@ -1321,7 +1317,7 @@ class SNES_Scalar(SolverBaseClass):
             bc = PetscDSAddBoundary_UW(cdm.dm,
                                 bc_type,
                                 str(boundary+f"{bc.components}").encode('utf8'),
-                                str(boundary).encode('utf8'),
+                                "UW_Boundaries".encode('utf8'),   # consolidated boundary label
                                 bc.f_id,  # field ID in the DM
                                 num_constrained_components,
                                 <const PetscInt *> &comps_view[0],
@@ -1540,7 +1536,46 @@ class SNES_Scalar(SolverBaseClass):
 
         ## Now add the boundary residual / jacobian terms
 
+        cdef DMLabel c_label
+
+        for bc in self.natural_bcs:
+
+            boundary = bc.boundary
+            boundary_id = bc.PETScID
+
+            value = self.mesh.boundaries[bc.boundary].value
+            bc_label = self.dm.getLabel("UW_Boundaries")
+
+            label_val = value
+
+            i_bd_res = self.ext_dict.bd_res
+            i_bd_jac = self.ext_dict.bd_jac
+
+            c_label = bc_label
 
+            if bc.fn_f is not None:
+
+                UW_PetscDSSetBdResidual(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0,
+                                ext.fns_bd_residual[i_bd_res[bc.fns["u_f0"]]],
+                                NULL,
+                                )
+
+                UW_PetscDSSetBdJacobian(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0, 0,
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G0"]]],
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G1"]]],
+                                NULL,
+                                NULL,
+                                )
+
+                UW_PetscDSSetBdJacobianPreconditioner(ds.ds, c_label.dmlabel, label_val, boundary_id,
+                                0, 0, 0,
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G0"]]],
+                                ext.fns_bd_jacobian[i_bd_jac[bc.fns["uu_G1"]]],
+                                NULL,
+                                NULL,
+                                )
 
         # Rebuild this lot
 
diff --git a/tests/test_1000_poissonNaturalBC.py b/tests/test_1000_poissonNaturalBC.py
index e14944af9..f3b26cc02 100644
--- a/tests/test_1000_poissonNaturalBC.py
+++ b/tests/test_1000_poissonNaturalBC.py
@@ -53,7 +53,6 @@
 )
 
 
-@pytest.mark.skip(reason="Poisson natural BC setup failing - needs solver investigation (PETSc error 73: Object in wrong state)")
 @pytest.mark.parametrize("mesh", [mesh_simp_reg, mesh_simp_irreg, mesh_quad])
 def test_poisson_natural_bc(mesh):
     """Test Poisson solver with natural (flux) boundary conditions."""