From aec3ff75fb6888ae000652951b6f685a2d7080c9 Mon Sep 17 00:00:00 2001 From: Juan Carlos Graciosa Date: Tue, 23 Sep 2025 21:24:50 +1000 Subject: [PATCH] Fix for bug #25 along with changes for natural BC in Scalar solver to work. Also added pytest for Poisson (scalar) natural BC --- .../cython/petsc_generic_snes_solvers.pyx | 43 ++++++----- tests/test_1000_poissonNaturalBC.py | 77 +++++++++++++++++++ 2 files changed, 99 insertions(+), 21 deletions(-) create mode 100644 tests/test_1000_poissonNaturalBC.py diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index a2070f00..b42b1668 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -700,7 +700,7 @@ class SNES_Scalar(SolverBaseClass): print(" - field: {}".format(bc.f_id)) print(" - components: {}".format(bc.components)) print(" - boundary: {}".format(bc.boundary)) - print(" - fn: {} ".format(bc.fn)) + print(" - fn: {} ".format(bc.fn_f)) boundary = bc.boundary value = mesh.boundaries[bc.boundary].value @@ -846,38 +846,39 @@ class SNES_Scalar(SolverBaseClass): # continue if bc.fn_f is not None: - + bd_F0 = sympy.Array(bc.fn_f) - self._bd_f0 = sympy.ImmutableDenseMatrix(bd_F0) - - G0 = sympy.derive_by_array(self._bd_f0, U) - G1 = sympy.derive_by_array(self._bd_f0, self.Unknowns.L) + bc.fns["u_f0"] = sympy.ImmutableDenseMatrix(bd_F0) + + G0 = sympy.derive_by_array(bd_F0, U) + G1 = sympy.derive_by_array(bd_F0, L) - self._bd_G0 = sympy.ImmutableMatrix(G0) # sympy.ImmutableMatrix(sympy.permutedims(G0, permutation).reshape(dim,dim)) - self._bd_G1 = sympy.ImmutableMatrix(G1.reshape(dim)) # sympy.ImmutableMatrix(sympy.permutedims(G1, permutation).reshape(dim,dim*dim)) + bc.fns["uu_G0"] = sympy.ImmutableMatrix(G0.reshape(1, 1)) + bc.fns["uu_G1"] = sympy.ImmutableMatrix(G1.reshape(dim, 1)) - fns_bd_residual += [self._bd_f0] - fns_bd_jacobian += [self._bd_G0, self._bd_G1] + fns_bd_residual += [bc.fns["u_f0"]] + fns_bd_jacobian += [bc.fns["uu_G0"], bc.fns["uu_G1"]] + # Similar to SNES_Vector, will leave these out for now, perhaps a different user-interface altogether is required for flux-like bcs - if bc.fn_F is not None: + # if bc.fn_F is not None: - bd_F1 = sympy.Array(bc.fn_F).reshape(dim) - self._bd_f1 = sympy.ImmutableDenseMatrix(bd_F1) + # bd_F1 = sympy.Array(bc.fn_F).reshape(dim) + # self._bd_f1 = sympy.ImmutableDenseMatrix(bd_F1) - G2 = sympy.derive_by_array(self._bd_f1, U) - G3 = sympy.derive_by_array(self._bd_f1, self.Unknowns.L) + # G2 = sympy.derive_by_array(self._bd_f1, U) + # G3 = sympy.derive_by_array(self._bd_f1, self.Unknowns.L) - self._bd_uu_G2 = sympy.ImmutableMatrix(G2.reshape(dim)) # sympy.ImmutableMatrix(sympy.permutedims(G2, permutation).reshape(dim*dim,dim)) - self._bd_uu_G3 = sympy.ImmutableMatrix(G3.reshape(dim,dim)) # sympy.ImmutableMatrix(sympy.permutedims(G3, permutation).reshape(dim*dim,dim*dim)) + # self._bd_uu_G2 = sympy.ImmutableMatrix(G2.reshape(dim)) # sympy.ImmutableMatrix(sympy.permutedims(G2, permutation).reshape(dim*dim,dim)) + # self._bd_uu_G3 = sympy.ImmutableMatrix(G3.reshape(dim,dim)) # sympy.ImmutableMatrix(sympy.permutedims(G3, permutation).reshape(dim*dim,dim*dim)) - fns_bd_residual += [self._bd_f1] - fns_bd_jacobian += [self._bd_G2, self._bd_G3] + # fns_bd_residual += [self._bd_f1] + # fns_bd_jacobian += [self._bd_G2, self._bd_G3] self._fns_bd_residual = fns_bd_residual self._fns_bd_jacobian = fns_bd_jacobian - + # generate JIT code. # first, we must specify the primary fields. # these are fields for which the corresponding sympy functions @@ -1284,7 +1285,7 @@ class SNES_Vector(SolverBaseClass): print(" - field: {}".format(bc.f_id)) print(" - component: {}".format(bc.components)) print(" - boundary: {}".format(bc.boundary)) - print(" - fn: {} ".format(bc.fn)) + print(" - fn: {} ".format(bc.fn_f)) boundary = bc.boundary value = mesh.boundaries[bc.boundary].value diff --git a/tests/test_1000_poissonNaturalBC.py b/tests/test_1000_poissonNaturalBC.py new file mode 100644 index 00000000..a3bb7e9c --- /dev/null +++ b/tests/test_1000_poissonNaturalBC.py @@ -0,0 +1,77 @@ +import pytest +import underworld3 as uw +import sympy +import numpy as np + +''' +Unit test for Natural BCs in a Poisson (scalar) problem. +''' + +res = 16 + +width = 1 +height = 1 + +minX, maxX = 0, width +minY, maxY = 0, height + +mesh_simp_reg = uw.meshing.UnstructuredSimplexBox( minCoords = (minX, minY), + maxCoords = (maxX, maxY), + cellSize = 1 / res, + qdegree = 3, + regular = True) +mesh_simp_irreg = uw.meshing.UnstructuredSimplexBox(minCoords = (minX, minY), + maxCoords =(maxX, maxY), + cellSize = 1 / res, + regular = False) +mesh_quad = uw.meshing.StructuredQuadBox(minCoords = (minX, minY), + maxCoords = (maxX, maxY), + elementRes = (res, res)) + +@pytest.mark.parametrize("mesh", + [mesh_simp_reg, + mesh_simp_irreg, + mesh_quad] +) + +def test_poisson_natural_bc(mesh): + + T_soln = uw.discretisation.MeshVariable("T", mesh, 1, degree = 2) + + poisson = uw.systems.Poisson(mesh = mesh, + u_Field = T_soln, + degree = 2, + verbose = False, + ) + + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.petsc_options.delValue("ksp_monitor") + + # set the source based on analytical solution + x, y = mesh.N.x, mesh.N.y + ana_soln = (x**2) * y + + #poisson.tolerance = 1e-6 # increase tolerance to decrease atol in allclose + poisson.f = -2 * y + + poisson.add_natural_bc(0.0, "Left") + poisson.add_natural_bc([2*y], "Right") + poisson.add_natural_bc([x**2], "Bottom") + poisson.add_natural_bc([x**2], "Top") + + poisson.add_dirichlet_bc(0.0, "Left") + poisson.add_dirichlet_bc([y], "Right") + poisson.add_dirichlet_bc(0.0, "Bottom") + poisson.add_dirichlet_bc([x**2], "Top") + + poisson.solve() + + with mesh.access(): + num = T_soln.data[:].squeeze() + ana = uw.function.evaluate(ana_soln, T_soln.coords) + + assert np.allclose(ana, num, atol = 1e-4), "Numerical and analytical solutions differ!" + + del poisson + del mesh \ No newline at end of file