From 86c66b8cd966f6e69a28e69778da8f923d0c48a5 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com> Date: Sun, 7 Jul 2024 14:09:27 -0400 Subject: [PATCH] Unit tests for Laplacian smoothing (#111) Closes #19. Closes #78. This will finish covering Movement with unit tests once #95 and #109 are merged. --- movement/laplacian_smoothing.py | 35 +++++-- test/test_forced_movement.py | 151 +++++++++++++++++++++++++++++++ test/test_laplacian_smoothing.py | 20 ++++ test/test_spring.py | 84 ----------------- 4 files changed, 196 insertions(+), 94 deletions(-) create mode 100644 test/test_forced_movement.py create mode 100644 test/test_laplacian_smoothing.py diff --git a/movement/laplacian_smoothing.py b/movement/laplacian_smoothing.py index d216a92..38c9e99 100644 --- a/movement/laplacian_smoothing.py +++ b/movement/laplacian_smoothing.py @@ -4,7 +4,7 @@ import ufl from firedrake.petsc import PETSc -import movement.solver_parameters as solver_parameters +import movement.solver_parameters as sp from movement.mover import PrimeMover __all__ = ["LaplacianSmoother"] @@ -50,7 +50,15 @@ def _create_functions(self): self.rhs = firedrake.Function(self.coord_space, name="Zero RHS") @PETSc.Log.EventDecorator() - def _setup_solver(self, boundary_conditions): + def _solve(self, boundary_conditions, solver_parameters=None): + """ + Solve the Laplace system. + + :kwarg boundary_conditions: Dirichlet boundary conditions to be enforced + :type boundary_conditions: :class:`~.DirichletBC` or :class:`list` thereof + :kwarg solver_parameters: solver parameters for solving the Laplace equation + :type solver_parameters: :class:`dict` + """ if not hasattr(self, "_solver"): test = firedrake.TestFunction(self.coord_space) trial = firedrake.TrialFunction(self.coord_space) @@ -62,23 +70,31 @@ def _setup_solver(self, boundary_conditions): ) self._solver = firedrake.LinearVariationalSolver( problem, - solver_parameters=solver_parameters.cg_ilu, + solver_parameters=solver_parameters or sp.cg_ilu, ) self._solver.solve() @PETSc.Log.EventDecorator() - def move(self, time, update_boundary_velocity=None, boundary_conditions=None): + def move( + self, + time, + update_boundary_velocity=None, + boundary_conditions=None, + solver_parameters=None, + ): """ Assemble and solve the Laplacian system and update the coordinates. :arg time: the current time :type time: :class:`float` - :kwarg update_boundary_velocity: function that updates the boundary conditions at - the current time + :kwarg update_boundary_velocity: function that updates the boundary conditions + at the current time :type update_boundary_velocity: :class:`~.Callable` with a single argument of :class:`float` type :kwarg boundary_conditions: Dirichlet boundary conditions to be enforced :type boundary_conditions: :class:`~.DirichletBC` or :class:`list` thereof + :kwarg solver_parameters: solver parameters for solving the Laplace equation + :type solver_parameters: :class:`dict` """ if update_boundary_velocity is not None: update_boundary_velocity(time) @@ -86,12 +102,11 @@ def move(self, time, update_boundary_velocity=None, boundary_conditions=None): boundary_conditions = firedrake.DirichletBC( self.coord_space, 0, "on_boundary" ) - self._setup_solver(boundary_conditions) # Solve on computational mesh self.mesh.coordinates.assign(self.xi) try: - self._solver.solve() + self._solve(boundary_conditions, solver_parameters=solver_parameters) except fexc.ConvergenceError as conv_err: self._convergence_error(exception=conv_err) @@ -99,8 +114,6 @@ def move(self, time, update_boundary_velocity=None, boundary_conditions=None): self.displacement[:] = self.v.dat.data_with_halos * self.dt self.x.dat.data_with_halos[:] += self.displacement self.mesh.coordinates.assign(self.x) - if hasattr(self, "tangling_checker"): - self.tangling_checker.check() self.volume.interpolate(ufl.CellVolume(self.mesh)) PETSc.Sys.Print( f"{time:.2f} s" @@ -108,3 +121,5 @@ def move(self, time, update_boundary_velocity=None, boundary_conditions=None): f" Variation (σ/μ) {self.coefficient_of_variation:8.2e}" f" Displacement {np.linalg.norm(self.displacement):.2f} m" ) + if hasattr(self, "tangling_checker"): + self.tangling_checker.check() diff --git a/test/test_forced_movement.py b/test/test_forced_movement.py new file mode 100644 index 0000000..6926eba --- /dev/null +++ b/test/test_forced_movement.py @@ -0,0 +1,151 @@ +import unittest + +import numpy as np +from firedrake import * + +from movement import LaplacianSmoother, SpringMover + + +class BaseClasses: + """ + Base classes for testing mesh movement under forcings. + """ + + class TestBoundaryForcing(unittest.TestCase): + """ + Unit tests for mesh movement under boundary forcings. + """ + + def move( + self, + mesh, + fixed_boundary_tags=None, + moving_boundary_tags=None, + vector=[1, 0], + **kwargs, + ): + mover = self.mover(mesh) + bcs = [] + if fixed_boundary_tags: + bcs.append(DirichletBC(mover.coord_space, 0, fixed_boundary_tags)) + if moving_boundary_tags: + f = Function(mover.coord_space).interpolate(as_vector(vector)) + bcs.append(DirichletBC(mover.coord_space, f, moving_boundary_tags)) + mover.move(0.0, boundary_conditions=bcs, **kwargs) + return mover.mesh + + @staticmethod + def shifted_mesh(nx, ny, shift_x=0, shift_y=0): + mesh = UnitSquareMesh(nx, ny) + x, y = SpatialCoordinate(mesh) + shifted_coords = Function(mesh.coordinates) + shifted_coords.interpolate(as_vector([x + shift_x, y + shift_y])) + return Mesh(shifted_coords) + + @staticmethod + def stretched_mesh(nx, ny, stretch_x=1, stretch_y=1): + return RectangleMesh(nx, ny, stretch_x, stretch_y) + + def test_fixed_triangle(self): + mesh = UnitTriangleMesh() + coords = mesh.coordinates.dat.data + self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data)) + + def test_fixed_square(self): + mesh = UnitSquareMesh(1, 1) + coords = mesh.coordinates.dat.data + self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data)) + + def test_force_right_free(self): + mesh = UnitSquareMesh(10, 10) + coord_array = mesh.coordinates.dat.data + new_coord_array = self.move( + mesh, moving_boundary_tags=1 + ).coordinates.dat.data + self.assertFalse(np.allclose(coord_array, new_coord_array)) + shifted_coord_array = self.shifted_mesh( + 10, 10, shift_x=1 + ).coordinates.dat.data + self.assertTrue(np.allclose(shifted_coord_array, new_coord_array)) + + def test_force_right_fixed(self): + mesh = UnitSquareMesh(10, 10) + coord_array = mesh.coordinates.dat.data + new_mesh = self.move( + mesh, moving_boundary_tags=1, fixed_boundary_tags=2, vector=[0.1, 0] + ) + new_coord_array = new_mesh.coordinates.dat.data + self.assertFalse(np.allclose(coord_array, new_coord_array)) + # # TODO: Implement no-slip BCs for segments 3 and 4 (#99) + # stretched_mesh = self.stretched_mesh(10, 10, stretch_x=2) + # stretched_coord_array = stretched_mesh.coordinates.dat.data + # self.assertTrue(np.allclose(stretched_coord_array, new_coord_array)) + + def test_force_right_left_free(self): + mesh = UnitSquareMesh(10, 10) + coord_array = mesh.coordinates.dat.data + mesh = self.move(mesh, moving_boundary_tags=1, vector=[1, 0]) + self.assertFalse(np.allclose(coord_array, mesh.coordinates.dat.data)) + mesh = self.move(mesh, moving_boundary_tags=1, vector=[-1, 0]) + self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data)) + + +class TestSpringMover(BaseClasses.TestBoundaryForcing): + """ + Unit tests for the lineal spring method under boundary forcings. + """ + + @staticmethod + def mover(mesh): + return SpringMover(mesh, 1.0) + + def test_update_boundary_displacement(self): + mesh = UnitSquareMesh(10, 10) + coord_array = mesh.coordinates.dat.data + mover = self.mover(mesh) + forcing = Function(mover.coord_space) + bc = DirichletBC(mover.coord_space, forcing, 1) + + def update_bc(time): + forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)])) + + mover.move(0.0, boundary_conditions=bc, update_boundary_displacement=update_bc) + self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data)) + mover.move(1.0, boundary_conditions=bc, update_boundary_displacement=update_bc) + self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data)) + + +class TestLaplacianSmoother(BaseClasses.TestBoundaryForcing): + """ + Unit tests for Laplacian smoothing under boundary forcings. + """ + + @staticmethod + def mover(mesh): + return LaplacianSmoother(mesh, 1.0) + + def test_convergence_error(self): + with self.assertRaises(ConvergenceError) as cm: + self.move( + UnitSquareMesh(10, 10), + moving_boundary_tags=1, + fixed_boundary_tags=2, + vector=[1.0, 0], + solver_parameters={"ksp_type": "cg", "ksp_max_it": 0}, + ) + self.assertEqual(str(cm.exception), "Solver failed to converge.") + + def test_update_boundary_velocity(self): + mesh = UnitSquareMesh(10, 10) + coord_array = mesh.coordinates.dat.data + mover = LaplacianSmoother(mesh, 1.0) + forcing = Function(mover.coord_space) + bc = DirichletBC(mover.coord_space, forcing, 1) + + def update_bc(time): + forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)])) + + mover.move(0.0, boundary_conditions=bc, update_boundary_velocity=update_bc) + self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data)) + mover.move(1.0, boundary_conditions=bc, update_boundary_velocity=update_bc) + self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data)) diff --git a/test/test_laplacian_smoothing.py b/test/test_laplacian_smoothing.py new file mode 100644 index 0000000..6726b3d --- /dev/null +++ b/test/test_laplacian_smoothing.py @@ -0,0 +1,20 @@ +import unittest +from unittest.mock import MagicMock + +from firedrake import * + +from movement import LaplacianSmoother + + +class TestExceptions(unittest.TestCase): + """ + Unit tests for exceptions raised during Laplacian smoothing. + """ + + def test_tangling_valueerror(self): + mover = LaplacianSmoother(UnitSquareMesh(3, 3), 1.0) + mover.x.dat.data[3] += 0.2 + mover._solve = MagicMock() + with self.assertRaises(ValueError) as cm: + mover.move(0) + self.assertEqual(str(cm.exception), "Mesh has 1 tangled element.") diff --git a/test/test_spring.py b/test/test_spring.py index 8109cbf..4a55998 100644 --- a/test/test_spring.py +++ b/test/test_spring.py @@ -146,87 +146,3 @@ def test_angles(self): """ angles = self.rad2deg(self.mover.angles.dat.data) self.assertTrue(np.allclose(angles, [0, 135, 90])) - - -class TestMovement(unittest.TestCase): - """ - Unit tests for movement under spring-based methods. - """ - - @staticmethod - def move(mesh, fixed_boundary_tags=None, moving_boundary_tags=None, vector=[1, 0]): - mover = SpringMover(mesh, 1.0) - bcs = [] - if fixed_boundary_tags: - bcs.append(DirichletBC(mover.coord_space, 0, fixed_boundary_tags)) - if moving_boundary_tags: - f = Function(mover.coord_space).interpolate(as_vector(vector)) - bcs.append(DirichletBC(mover.coord_space, f, moving_boundary_tags)) - mover.move(0.0, boundary_conditions=bcs) - return mover.mesh - - @staticmethod - def shifted_mesh(nx, ny, shift_x=0, shift_y=0): - mesh = UnitSquareMesh(nx, ny) - x, y = SpatialCoordinate(mesh) - shifted_coords = Function(mesh.coordinates) - shifted_coords.interpolate(as_vector([x + shift_x, y + shift_y])) - return Mesh(shifted_coords) - - @staticmethod - def stretched_mesh(nx, ny, stretch_x=1, stretch_y=1): - return RectangleMesh(nx, ny, stretch_x, stretch_y) - - def test_fixed_triangle(self): - mesh = UnitTriangleMesh() - coords = mesh.coordinates.dat.data - self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data)) - - def test_fixed_square(self): - mesh = UnitSquareMesh(1, 1) - coords = mesh.coordinates.dat.data - self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data)) - - def test_force_right_free(self): - mesh = UnitSquareMesh(10, 10) - coord_array = mesh.coordinates.dat.data - new_coord_array = self.move(mesh, moving_boundary_tags=1).coordinates.dat.data - self.assertFalse(np.allclose(coord_array, new_coord_array)) - shifted_coord_array = self.shifted_mesh(10, 10, shift_x=1).coordinates.dat.data - self.assertTrue(np.allclose(shifted_coord_array, new_coord_array)) - - def test_force_right_fixed(self): - mesh = UnitSquareMesh(10, 10) - coord_array = mesh.coordinates.dat.data - new_mesh = self.move( - mesh, moving_boundary_tags=1, fixed_boundary_tags=2, vector=[0.1, 0] - ) - new_coord_array = new_mesh.coordinates.dat.data - self.assertFalse(np.allclose(coord_array, new_coord_array)) - # # TODO: Implement no-slip BCs for segments 3 and 4 (#99) - # stretched_mesh = self.stretched_mesh(10, 10, stretch_x=2) - # stretched_coord_array = stretched_mesh.coordinates.dat.data - # self.assertTrue(np.allclose(stretched_coord_array, new_coord_array)) - - def test_force_right_left_free(self): - mesh = UnitSquareMesh(10, 10) - coord_array = mesh.coordinates.dat.data - mesh = self.move(mesh, moving_boundary_tags=1, vector=[1, 0]) - self.assertFalse(np.allclose(coord_array, mesh.coordinates.dat.data)) - mesh = self.move(mesh, moving_boundary_tags=1, vector=[-1, 0]) - self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data)) - - def test_update_boundary_displacement(self): - mesh = UnitSquareMesh(10, 10) - coord_array = mesh.coordinates.dat.data - mover = SpringMover(mesh, 1.0) - forcing = Function(mover.coord_space) - bc = DirichletBC(mover.coord_space, forcing, 1) - - def update_bc(time): - forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)])) - - mover.move(0.0, boundary_conditions=bc, update_boundary_displacement=update_bc) - self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data)) - mover.move(1.0, boundary_conditions=bc, update_boundary_displacement=update_bc) - self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))