Skip to content

Fix for bug #25 #26

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
jcgraciosa merged 1 commit into from
Oct 8, 2025
Merged

Fix for bug #25 #26

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
43 changes: 22 additions & 21 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 @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
77 changes: 77 additions & 0 deletions tests/test_1000_poissonNaturalBC.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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