Skip to content
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

Support float and double precision, as well as test complex numbers #19

Merged
merged 1 commit into from
Sep 18, 2024
Merged

Support float and double precision, as well as test complex numbers #19

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
4 changes: 2 additions & 2 deletions src/scifem/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@ def create_real_functionspace(mesh: dolfinx.mesh.Mesh, value_shape: tuple[int, .

"""


ufl_e = basix.ufl.element("P", mesh.basix_cell(), 0, dtype=float, discontinuous=True,
dtype = mesh.geometry.x.dtype
ufl_e = basix.ufl.element("P", mesh.basix_cell(), 0, dtype=dtype, discontinuous=True,
shape=value_shape)

if (dtype:=mesh.geometry.x.dtype) == np.float64:
Expand Down
82 changes: 58 additions & 24 deletions tests/test_real_functionspace.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,70 +11,73 @@
@pytest.mark.parametrize("L", [0.1, 0.2, 0.3])
@pytest.mark.parametrize("H", [1.3, 0.8, 0.2])
@pytest.mark.parametrize("cell_type", [dolfinx.mesh.CellType.triangle, dolfinx.mesh.CellType.quadrilateral])
def test_real_function_space_mass(L, H, cell_type):
@pytest.mark.parametrize("dtype", [np.float64, np.float32])
def test_real_function_space_mass(L, H, cell_type, dtype):
"""
Check that real space mass matrix is the same as assembling the volume of the mesh
"""

mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[0.,0.],[L,H]], [7,9],cell_type)
mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[0.,0.],[L,H]], [7,9],cell_type, dtype=dtype)

V = scifem.create_real_functionspace(mesh)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(u, v) * ufl.dx

A = dolfinx.fem.assemble_matrix(dolfinx.fem.form(a), bcs=[])
A = dolfinx.fem.assemble_matrix(dolfinx.fem.form(a, dtype=dtype), bcs=[])
A.scatter_reverse()

tol = 100 * np.finfo(dtype).eps

assert len(A.data) == 1
if MPI.COMM_WORLD.rank == 0:
assert np.isclose(A.data[0], L*H)

assert np.isclose(A.data[0], L*H, atol=tol)

@pytest.mark.parametrize("dtype", [np.float64, np.float32])
@pytest.mark.parametrize("cell_type", [dolfinx.mesh.CellType.tetrahedron, dolfinx.mesh.CellType.hexahedron])
def test_real_function_space_vector(cell_type):
def test_real_function_space_vector(cell_type, dtype):
"""
Test that assembling against a real space test function is equivalent to assembling a vector
"""


mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 2,3,5, cell_type)
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 2,3,5, cell_type, dtype=dtype)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))
v = ufl.TrialFunction(V)

R = scifem.create_real_functionspace(mesh)
u = ufl.TestFunction(R)
a_R = ufl.inner(u, v) * ufl.dx
form_rhs = dolfinx.fem.form(a_R)
form_rhs = dolfinx.fem.form(a_R, dtype=dtype)

A_R = dolfinx.fem.assemble_matrix(form_rhs, bcs=[])
A_R.scatter_reverse()

L = ufl.inner(ufl.constantvalue.IntValue(1), v) * ufl.dx
form_lhs = dolfinx.fem.form(L)
form_lhs = dolfinx.fem.form(L, dtype=dtype)
b = dolfinx.fem.assemble_vector(form_lhs)
b.scatter_reverse(dolfinx.la.InsertMode.add)
b.scatter_forward()

row_map = A_R.index_map(0)
num_local_rows = row_map.size_local
num_dofs = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
tol = 100 * np.finfo(dtype).eps
if MPI.COMM_WORLD.rank == 0:
assert num_local_rows == 1
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
assert A_R.indptr[1] - A_R.indptr[0] == num_dofs_global
np.testing.assert_allclose(A_R.indices, np.arange(num_dofs_global))
np.testing.assert_allclose(b.array[:num_dofs], A_R.data[:num_dofs])
np.testing.assert_allclose(b.array[:num_dofs], A_R.data[:num_dofs], atol=tol)
else:
assert num_local_rows == 0


@pytest.mark.parametrize("dtype", [PETSc.RealType])
@pytest.mark.parametrize("tensor", [0, 1, 2])
@pytest.mark.parametrize("degree", range(1, 5))
def test_singular_poisson(tensor, degree):
def test_singular_poisson(tensor, degree, dtype):
M = 25
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle)
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, M, M, dolfinx.mesh.CellType.triangle, dtype=dtype)

if tensor == 0:
value_shape = ()
Expand All @@ -93,21 +96,21 @@ def test_singular_poisson(tensor, degree):
x = ufl.SpatialCoordinate(mesh)
pol = x[0]**degree - 2*x[1]**degree
# Compute average value of polynomial to make mean 0
C = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(pol*ufl.dx)), op=MPI.SUM)
u_scalar = pol - dolfinx.fem.Constant(mesh, C)
C = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(pol*ufl.dx, dtype=dtype)), op=MPI.SUM)
u_scalar = pol - dolfinx.fem.Constant(mesh, dtype(C))
if tensor == 0:
u_ex = u_scalar
zero = dolfinx.fem.Constant(mesh, 0.0)
zero = dolfinx.fem.Constant(mesh, dtype(0.0))
elif tensor == 1:
u_ex = ufl.as_vector([u_scalar, -u_scalar])
zero = dolfinx.fem.Constant(mesh, (0.0,0.0))
zero = dolfinx.fem.Constant(mesh, dtype((0.0,0.0)))
else:
u_ex = ufl.as_tensor([[u_scalar, 2*u_scalar], [3*u_scalar, -u_scalar],
[u_scalar, 2*u_scalar],
])
zero = dolfinx.fem.Constant(mesh, ((0.0,0.0),
zero = dolfinx.fem.Constant(mesh, dtype(((0.0,0.0),
(0.0,0.0),
(0.0,0.0)))
(0.0,0.0))))

dx = ufl.Measure("dx", domain=mesh)
f = -ufl.div(ufl.grad(u_ex))
Expand All @@ -120,8 +123,8 @@ def test_singular_poisson(tensor, degree):
L0 = ufl.inner(f , v) * dx + ufl.inner(g, v) * ufl.ds
L1 = ufl.inner(zero, d) * dx

a = dolfinx.fem.form([[a00, a01], [a10, None]])
L = dolfinx.fem.form([L0, L1])
a = dolfinx.fem.form([[a00, a01], [a10, None]], dtype=dtype)
L = dolfinx.fem.form([L0, L1], dtype=dtype)


A = dolfinx.fem.petsc.assemble_matrix_block(a)
Expand All @@ -148,8 +151,39 @@ def test_singular_poisson(tensor, degree):
uh.x.array[:len(x_local[0])] = x_local[0]
uh.x.scatter_forward()

error = dolfinx.fem.form(ufl.inner(u_ex - uh, u_ex - uh) * dx)
error = dolfinx.fem.form(ufl.inner(u_ex - uh, u_ex - uh) * dx, dtype=dtype)

e_local = dolfinx.fem.assemble_scalar(error)
tol = 2e3 * np.finfo(dtype).eps
e_global = np.sqrt(mesh.comm.allreduce(e_local, op=MPI.SUM))
assert np.isclose(e_global, 0)
assert np.isclose(e_global, 0, atol=tol)


@pytest.mark.parametrize("ftype, stype", [(np.float32, np.complex64), (np.float64, np.complex128)])
def test_complex_real_space(ftype, stype):
mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 13, dtype=ftype)

val = (2 + 3j, -4 + 5j)
value_shape = (2,)
R = scifem.create_real_functionspace(mesh, value_shape=value_shape)
u = dolfinx.fem.Function(R, dtype=stype)
u.x.array[0] = val[0]
u.x.array[1] = val[1]
u.x.scatter_forward()

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, value_shape))
v = ufl.TestFunction(V)
L = ufl.inner(u, v) * ufl.dx

b = dolfinx.fem.assemble_vector(dolfinx.fem.form(L, dtype=stype))
b.scatter_reverse(dolfinx.la.InsertMode.add)
b.scatter_forward()
const = dolfinx.fem.Constant(mesh, stype(val))
L_const = ufl.inner(const, v) * ufl.dx
b_const = dolfinx.fem.assemble_vector(dolfinx.fem.form(L_const, dtype=stype))
b_const.scatter_reverse(dolfinx.la.InsertMode.add)
b_const.scatter_forward()


tol = 100 * np.finfo(stype).eps
np.testing.assert_allclose(b.array, b_const.array, atol=tol)