Skip to content

Fix Poisson solver natural (Neumann) boundary conditions #60

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lmoresi merged 1 commit into from
Feb 28, 2026
Merged

Fix Poisson solver natural (Neumann) boundary conditions #60

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 40 additions & 5 deletions src/underworld3/cython/petsc_generic_snes_solvers.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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],
Expand Down Expand Up @@ -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

Expand Down
1 change: 0 additions & 1 deletion tests/test_1000_poissonNaturalBC.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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."""
Expand Down