From 11a699e400a400a1445bc1fc27ffad4d62e01271 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com> Date: Fri, 10 May 2024 18:03:59 +0100 Subject: [PATCH] Overhaul Monge-Ampere tests (#80) Partially addresses #19. Closes #74. Partially addresses #41. Partially addresses #78. Closes #87. This PR overhauls the Monge-Ampere testing framework to use `unittest` and cover much more of the functionality. It also cuts down the test suite run time significantly by running CI tests at lower resolution and with larger tolerances. --- Makefile | 3 +- demos/monge_ampere1.py | 18 +- demos/monge_ampere_3d.py | 8 +- demos/monge_ampere_helmholtz.py | 12 +- movement/laplacian_smoothing.py | 6 +- movement/monge_ampere.py | 377 +++++++++++++++----------------- movement/mover.py | 2 +- movement/solver_parameters.py | 11 +- test/test_monge_ampere.py | 302 +++++++++++++++---------- 9 files changed, 407 insertions(+), 332 deletions(-) diff --git a/Makefile b/Makefile index 606b385..7e269cb 100644 --- a/Makefile +++ b/Makefile @@ -23,7 +23,8 @@ test: lint coverage: @echo "Generating coverage report..." @python3 -m coverage erase - @python3 -m coverage run --source=movement -m pytest -v test + @python3 -m coverage run --source=movement -m pytest -v test \ + --durations=20 @python3 -m coverage html @echo "Done." diff --git a/demos/monge_ampere1.py b/demos/monge_ampere1.py index 662b4a3..b615e1f 100644 --- a/demos/monge_ampere1.py +++ b/demos/monge_ampere1.py @@ -39,13 +39,17 @@ # # We begin the example by importing from the namespaces of Firedrake and Movement. +# To start with a simple example, consider a uniform mesh of the unit square. Feel free +# to ignore the `"MOVEMENT_REGRESSION_TEST"`, as it is only used when this demo is run +# in the test suite (to reduce its runtime). :: +import os + from firedrake import * from movement import * -# To start with a simple example, consider a uniform mesh of the unit square. - -n = 20 +test = os.environ.get("MOVEMENT_REGRESSION_TEST") +n = 10 if test else 20 mesh = UnitSquareMesh(n, n) # We can plot the initial mesh using Matplotlib as follows. @@ -87,9 +91,13 @@ def ring_monitor(mesh): # With an initial mesh and a monitor function, we are able to construct a -# :class:`~.MongeAmpereMover` instance and adapt the mesh. +# :class:`~.MongeAmpereMover` instance and adapt the mesh. By default, the Monge-Ampère +# equation is solved to a relative tolerance of :math:`10^{-8}`. However, for the +# purposes of continuous integration testing, a tolerance of :math:`10^{-3}` is used +# instead to further reduce the runtime. :: -mover = MongeAmpereMover(mesh, ring_monitor, method="quasi_newton") +rtol = 1.0e-03 if test else 1.0e-08 +mover = MongeAmpereMover(mesh, ring_monitor, method="quasi_newton", rtol=rtol) mover.move() # The adapted mesh can be accessed via the `mesh` attribute of the mover. Plotting it, diff --git a/demos/monge_ampere_3d.py b/demos/monge_ampere_3d.py index b62fee2..4ea4c5a 100644 --- a/demos/monge_ampere_3d.py +++ b/demos/monge_ampere_3d.py @@ -72,9 +72,13 @@ def monitor(mesh): # we use the `"relaxation"` method (see :cite:`McRae:2018`), # which gives faster convergence for this case. -n = 20 +import os + +test = os.environ.get("MOVEMENT_REGRESSION_TEST") +rtol = 1.0e-03 if test else 1.0e-08 +n = 10 if test else 20 mesh = UnitCubeMesh(n, n, n) -mover = movement.MongeAmpereMover(mesh, monitor, method="relaxation") +mover = movement.MongeAmpereMover(mesh, monitor, method="relaxation", rtol=rtol) mover.move() # The results will be written to the `monitor.pvd` file which represents diff --git a/demos/monge_ampere_helmholtz.py b/demos/monge_ampere_helmholtz.py index bd4b8bf..3e1f063 100644 --- a/demos/monge_ampere_helmholtz.py +++ b/demos/monge_ampere_helmholtz.py @@ -41,11 +41,16 @@ # Because our chosen solution does not satisfy homogeneous Neumann boundary conditions, # we instead apply Dirichlet boundary conditions based on the chosen analytical solution. +import os + from firedrake import * from movement import MongeAmpereMover -mesh = UnitSquareMesh(20, 20) # initial mesh +test = os.environ.get("MOVEMENT_REGRESSION_TEST") +n = 10 if test else 20 + +mesh = UnitSquareMesh(n, n) # initial mesh def u_exact(mesh): @@ -132,7 +137,8 @@ def monitor(mesh): # :figwidth: 60% # :align: center -mover = MongeAmpereMover(mesh, monitor, method="quasi_newton") +rtol = 1.0e-03 if test else 1.0e-08 +mover = MongeAmpereMover(mesh, monitor, method="quasi_newton", rtol=rtol) mover.move() # For every iteration the MongeAmpereMover prints the minimum to maximum ratio of @@ -219,7 +225,7 @@ def monitor2(mesh): return m -mover = MongeAmpereMover(mesh, monitor2, method="quasi_newton") +mover = MongeAmpereMover(mesh, monitor2, method="quasi_newton", rtol=rtol) mover.move() u_h = solve_helmholtz(mover.mesh) diff --git a/movement/laplacian_smoothing.py b/movement/laplacian_smoothing.py index 39fe578..1ee752e 100644 --- a/movement/laplacian_smoothing.py +++ b/movement/laplacian_smoothing.py @@ -58,7 +58,7 @@ def _setup_solver(self, boundary_conditions): ) self._solver = firedrake.LinearVariationalSolver( problem, - solver_parameters=solver_parameters.cg, + solver_parameters=solver_parameters.cg_ilu, ) self._solver.solve() @@ -90,5 +90,5 @@ def move(self, time, update_boundary_velocity=None, boundary_conditions=None): # Update mesh coordinates self.displacement[:] = self.v.dat.data_with_halos * self.dt - self._x.dat.data_with_halos[:] += self.displacement - self.mesh.coordinates.assign(self._x) + self.x.dat.data_with_halos[:] += self.displacement + self.mesh.coordinates.assign(self.x) diff --git a/movement/monge_ampere.py b/movement/monge_ampere.py index 8f2fe38..715cee2 100644 --- a/movement/monge_ampere.py +++ b/movement/monge_ampere.py @@ -1,3 +1,6 @@ +import abc +from warnings import warn + import firedrake import numpy as np import ufl @@ -11,7 +14,6 @@ "MongeAmpereMover_Relaxation", "MongeAmpereMover_QuasiNewton", "MongeAmpereMover", - "monge_ampere", ] @@ -22,38 +24,45 @@ def tangential(v, n): def MongeAmpereMover(mesh, monitor_function, method="relaxation", **kwargs): r""" - Movement of a `mesh` is determined by a `monitor_function` - :math:`m` and the Monge-Ampère type equation - - .. math:: - m(x)\det(I + H(\phi)) = \theta, - - for a scalar potential :math:`\phi`, where :math:`I` is the - identity matrix, :math:`\theta` is a normalisation coefficient - and :math:`H(\phi)` denotes the Hessian of :math:`\phi` with - respect to the coordinates :math:`\xi` of the computational mesh. - - The physical mesh coordinates :math:`x` are updated according to - - .. math:: - x = \xi + \nabla\phi. - - :arg mesh: the physical mesh - :arg monitor_function: a Python function which takes a mesh as input - :kwarg method: choose from 'relaxation' and 'quasi_newton' - :kwarg phi_init: initial guess for the scalar potential - :kwarg sigma_init: initial guess for the Hessian - :return: converged scalar potential and Hessian + Movement of a `mesh` is determined by a `monitor_function` + :math:`m` and the Monge-Ampère type equation + + .. math:: + m(x)\det(I + H(\phi)) = \theta, + + for a scalar potential :math:`\phi`, where :math:`I` is the + identity matrix, :math:`\theta` is a normalisation coefficient + and :math:`H(\phi)` denotes the Hessian of :math:`\phi` with + respect to the coordinates :math:`\xi` of the computational mesh. + + The physical mesh coordinates :math:`x` are updated according to + + .. math:: + x = \xi + \nabla\phi. + + :arg mesh: the physical mesh + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :arg monitor_function: a Python function which takes a mesh as input + :type monitor_function: :class:`~.Callable` + :kwarg method: choose from 'relaxation' and 'quasi_newton' + :type method: :class:`str` + :kwarg phi_init: initial guess for the scalar potential + :type phi_init: :class:`firedrake.function.Function` + :kwarg sigma_init: initial guess for the Hessian + :type sigma_init: :class:`firedrake.function.Function` + :return: the Monge-Ampere Mover object + :rtype: :class:`MongeAmpereMover_Relaxation` or + :class:`MongeAmpereMover_QuasiNewton` """ if method == "relaxation": return MongeAmpereMover_Relaxation(mesh, monitor_function, **kwargs) elif method == "quasi_newton": return MongeAmpereMover_QuasiNewton(mesh, monitor_function, **kwargs) else: - raise ValueError(f"Method {method} not recognised.") + raise ValueError(f"Method '{method}' not recognised.") -class MongeAmpereMover_Base(PrimeMover): +class MongeAmpereMover_Base(PrimeMover, metaclass=abc.ABCMeta): """ Base class for mesh movers based on the solution of Monge-Ampere type equations. @@ -62,19 +71,26 @@ class MongeAmpereMover_Base(PrimeMover): def __init__(self, mesh, monitor_function, **kwargs): """ :arg mesh: the physical mesh + :type mesh: :class:`firedrake.mesh.MeshGeometry` :arg monitor_function: a Python function which takes a mesh as input + :type monitor_function: :class:`~.Callable` :kwarg phi_init: initial guess for the scalar potential + :type phi_init: :class:`firedrake.function.Function` :kwarg sigma_init: initial guess for the Hessian + :type sigma_init: :class:`firedrake.function.Function` :kwarg maxiter: maximum number of iterations for the relaxation + :type maxiter: :class:`int` :kwarg rtol: relative tolerance for the residual + :type rtol: :class:`float` :kwarg dtol: divergence tolerance for the residual + :type dtol: :class:`float` :kwarg fix_boundary_nodes: should all boundary nodes remain fixed? + :type fix_boundary_nodes: :class:`bool` """ if monitor_function is None: - raise ValueError("Please supply a monitor function") + raise ValueError("Please supply a monitor function.") # Collect parameters before calling super - self.pseudo_dt = firedrake.Constant(kwargs.pop("pseudo_timestep", 0.1)) self.maxiter = kwargs.pop("maxiter", 1000) self.rtol = kwargs.pop("rtol", 1.0e-08) self.dtol = kwargs.pop("dtol", 2.0) @@ -88,41 +104,42 @@ def __init__(self, mesh, monitor_function, **kwargs): self.P1_ten = firedrake.TensorFunctionSpace(self.mesh, "CG", 1) # Create objects used during the mesh movement + self._create_functions() self.theta = firedrake.Constant(0.0) - self.monitor = firedrake.Function(self.P1, name="Monitor function") self.monitor.interpolate(self.monitor_function(self.mesh)) - self.volume = firedrake.Function(self.P0, name="Mesh volume") self.volume.interpolate(ufl.CellVolume(self.mesh)) self.original_volume = firedrake.Function(self.volume) self.total_volume = firedrake.assemble(firedrake.Constant(1.0) * self.dx) self.L_P0 = firedrake.TestFunction(self.P0) * self.monitor * self.dx + + @abc.abstractmethod + def _create_functions(self): + self.monitor = firedrake.Function(self.P1, name="Monitor function") + self.volume = firedrake.Function(self.P0, name="Mesh volume") self._grad_phi = firedrake.Function(self.P1_vec) self.grad_phi = firedrake.Function(self.mesh.coordinates) - @PETSc.Log.EventDecorator("MongeAmpereBase.apply_initial_guess") - def apply_initial_guess(self, phi_init=None, sigma_init=None, **kwargs): + @PETSc.Log.EventDecorator() + def apply_initial_guess(self, phi_init, sigma_init): """ - Initialise the approximations to the scalar potential - and its hessian with an initial guess. + Initialise the approximations to the scalar potential and its Hessian with an + initial guess. - By default, both are initialised to zero, which corresponds - to the case where the computational and physical meshes - coincide. + By default, both are initialised to zero, which corresponds to the case where + the computational and physical meshes coincide. - :kwarg phi_init: initial guess for the scalar potential - :kwarg sigma_init: initial guess for the Hessian + :arg phi_init: initial guess for the scalar potential + :type phi_init: :class:`firedrake.function.Function` + :arg sigma_init: initial guess for the Hessian + :type sigma_init: :class:`firedrake.function.Function` """ if phi_init is not None and sigma_init is not None: - assert hasattr(self, "phi") - assert hasattr(self, "sigma") - assert hasattr(self, "phi_old") - assert hasattr(self, "sigma_old") - self.phi.assign(phi_init) - self.sigma.assign(sigma_init) - self.phi_old.assign(phi_init) - self.sigma_old.assign(sigma_init) + self.phi.project(phi_init) + self.sigma.project(sigma_init) + self.phi_old.project(phi_init) + self.sigma_old.project(sigma_init) elif phi_init is not None or sigma_init is not None: - raise ValueError("Need to initialise both phi *and* sigma") + raise ValueError("Need to initialise both phi *and* sigma.") @property @PETSc.Log.EventDecorator() @@ -147,32 +164,34 @@ def _diagnostics(self): residual_l2_rel = residual_l2 / norm_l2 return minmax, residual_l2_rel, cv - @property - @PETSc.Log.EventDecorator("MongeAmpereBase.update_coordinates") - def x(self): - """ - Update the coordinate :class:`Function` using - the recovered gradient. + @PETSc.Log.EventDecorator() + def _update_coordinates(self): + r""" + Update the physical coordinates :math:`\mathbf{x}` using the recovered gradient: + .. math:: + \mathbf{x} = \boldsymbol{\xi} + \nabla_{\boldsymbol{\xi}}\phi. """ try: self.grad_phi.assign(self._grad_phi) except Exception: self.grad_phi.interpolate(self._grad_phi) - self._x.assign(self.xi + self.grad_phi) # x = ξ + grad(φ) - return self._x + self.x.assign(self.xi + self.grad_phi) + self.mesh.coordinates.assign(self.x) @property - @PETSc.Log.EventDecorator("MongeAmpereBase.create_l2_projector") + @PETSc.Log.EventDecorator() def l2_projector(self): """ - Create a linear solver for obtaining the gradient - of the potential using an L2 projection. + Create a linear solver for obtaining the gradient of the potential using an + :math:`L^2` projection. Boundary conditions are imposed as a post-processing step. + + :return: the linear solver + :rtype: :class:`~.LinearVariationalSolver` """ if hasattr(self, "_l2_projector"): return self._l2_projector - assert hasattr(self, "phi_old") u_cts = firedrake.TrialFunction(self.P1_vec) v_cts = firedrake.TestFunction(self.P1_vec) @@ -183,26 +202,26 @@ def l2_projector(self): # Enforce no movement normal to boundary n = ufl.FacetNormal(self.mesh) bcs = [] - for i in self.mesh.exterior_facets.unique_markers: + for tag in self.mesh.exterior_facets.unique_markers: + # TODO: Write tests for the boundary conditions block below (#79) if self.fix_boundary_nodes: - bcs.append(firedrake.DirichletBC(self.P1_vec, 0, i)) + bcs.append(firedrake.DirichletBC(self.P1_vec, 0, tag)) continue # Check for axis-aligned boundaries - _n = [firedrake.assemble(abs(n[j]) * self.ds(i)) for j in range(self.dim)] + _n = [firedrake.assemble(abs(n[j]) * self.ds(tag)) for j in range(self.dim)] iszero = [np.allclose(ni, 0.0) for ni in _n] nzero = sum(iszero) - if nzero == self.dim: - raise ValueError(f"Invalid normal vector {_n}") - elif nzero == self.dim - 1: + assert nzero < self.dim + if nzero == self.dim - 1: idx = iszero.index(False) - bcs.append(firedrake.DirichletBC(self.P1_vec.sub(idx), 0, i)) + bcs.append(firedrake.DirichletBC(self.P1_vec.sub(idx), 0, tag)) continue # Enforce no mesh movement normal to boundaries a_bc = ufl.dot(v_cts, n) * ufl.dot(u_cts, n) * self.ds L_bc = ufl.dot(v_cts, n) * firedrake.Constant(0.0) * self.ds - bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, i)) + bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, tag)) # Allow tangential movement, but only up until the end of boundary segments a_bc = ufl.dot(tangential(v_cts, n), tangential(u_cts, n)) * self.ds @@ -214,71 +233,61 @@ def l2_projector(self): if len(edges) == 0: bbc = None # Periodic case else: - from warnings import warn - warn( "Have you checked that all straight line segments are uniquely tagged?" ) corners = [(i, j) for i in edges for j in edges.difference([i])] bbc = firedrake.DirichletBC(self.P1_vec, 0, corners) - bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, i, bcs=bbc)) + bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, tag, bcs=bbc)) # Create solver problem = firedrake.LinearVariationalProblem(a, L, self._grad_phi, bcs=bcs) - sp = { - "ksp_type": "cg", - "pc_type": "bjacobi", - "sub_pc_type": "ilu", - } self._l2_projector = firedrake.LinearVariationalSolver( - problem, solver_parameters=sp + problem, solver_parameters=solver_parameters.cg_ilu ) return self._l2_projector class MongeAmpereMover_Relaxation(MongeAmpereMover_Base): r""" - The elliptic Monge-Ampere equation is solved in a parabolised - form using a pseudo-time relaxation, - - .. math:: - -\frac\partial{\partial\tau}\Delta\phi = m(x)\det(I + H(\phi)) - \theta, - - where :math:`\tau` is the pseudo-time variable. Forward Euler is - used for the pseudo-time integration (see McRae et al. 2018 for - details). - - References - ========== - A. T. T. McRae, C. J. Cotter, C. J. Budd, Optimal-transport-based - mesh adaptivity on the plane and sphere using finite elements, SIAM - Journal on Scientific Computing 40 (2) (2018) 1121–1148. - doi:10.1137/16M1109515. + The elliptic Monge-Ampere equation is solved in a parabolised form using a + pseudo-time relaxation, + + .. math:: + -\frac\partial{\partial\tau}\Delta\phi = m(x)\det(I + H(\phi)) - \theta, + + where :math:`\tau` is the pseudo-time variable. Forward Euler is used for the + pseudo-time integration (see :cite:`MCB:18` for details). """ - @PETSc.Log.EventDecorator("MongeAmpereMover.__init__") - def __init__(self, mesh, monitor_function, **kwargs): + @PETSc.Log.EventDecorator() + def __init__( + self, mesh, monitor_function, phi_init=None, sigma_init=None, **kwargs + ): """ :arg mesh: the physical mesh + :type mesh: :class:`firedrake.mesh.MeshGeometry` :arg monitor_function: a Python function which takes a mesh as input + :type monitor_function: :class:`~.Callable` :kwarg phi_init: initial guess for the scalar potential + :type phi_init: :class:`firedrake.function.Function` :kwarg sigma_init: initial guess for the Hessian + :type sigma_init: :class:`firedrake.function.Function` :kwarg pseudo_timestep: pseudo-timestep to use for the relaxation + :type pseudo_timestep: :class:`float` :kwarg maxiter: maximum number of iterations for the relaxation + :type maxiter: :class:`int` :kwarg rtol: relative tolerance for the residual + :type rtol: :class:`float` :kwarg dtol: divergence tolerance for the residual + :type dtol: :class:`float` """ self.pseudo_dt = firedrake.Constant(kwargs.pop("pseudo_timestep", 0.1)) super().__init__(mesh, monitor_function=monitor_function, **kwargs) - # Create functions to hold solution data - self.phi = firedrake.Function(self.P1) - self.sigma = firedrake.Function(self.P1_ten) - self.phi_old = firedrake.Function(self.P1) - self.sigma_old = firedrake.Function(self.P1_ten) - # Initialise phi and sigma - self.apply_initial_guess(**kwargs) + if phi_init or sigma_init: + self.apply_initial_guess(phi_init, sigma_init) # Setup residuals I = ufl.Identity(self.dim) @@ -288,16 +297,24 @@ def __init__(self, mesh, monitor_function, **kwargs): self._residual_l2_form = psi * self.residual * self.dx self._norm_l2_form = psi * self.theta * self.dx + def _create_functions(self): + super()._create_functions() + self.phi = firedrake.Function(self.P1) + self.sigma = firedrake.Function(self.P1_ten) + self.phi_old = firedrake.Function(self.P1) + self.sigma_old = firedrake.Function(self.P1_ten) + @property - @PETSc.Log.EventDecorator("MongeAmpereMover.create_pseudotimestepper") + @PETSc.Log.EventDecorator() def pseudotimestepper(self): """ Setup the pseudo-timestepper for the relaxation method. + + :return: the pseudo-timestepper + :rtype: :class:`~.LinearVariationalSolver` """ if hasattr(self, "_pseudotimestepper"): return self._pseudotimestepper - assert hasattr(self, "phi") - assert hasattr(self, "phi_old") phi = firedrake.TrialFunction(self.P1) psi = firedrake.TestFunction(self.P1) a = ufl.inner(ufl.grad(psi), ufl.grad(phi)) * self.dx @@ -306,29 +323,26 @@ def pseudotimestepper(self): + self.pseudo_dt * psi * self.residual * self.dx ) problem = firedrake.LinearVariationalProblem(a, L, self.phi) - sp = { - "ksp_type": "cg", - "pc_type": "gamg", - } nullspace = firedrake.VectorSpaceBasis(constant=True) self._pseudotimestepper = firedrake.LinearVariationalSolver( problem, - solver_parameters=sp, + solver_parameters=solver_parameters.cg_gamg, nullspace=nullspace, transpose_nullspace=nullspace, ) return self._pseudotimestepper @property - @PETSc.Log.EventDecorator("MongeAmpereMover.create_equidistributor") + @PETSc.Log.EventDecorator() def equidistributor(self): """ Setup the equidistributor for the relaxation method. + + :return: the equidistributor + :rtype: :class:`~.LinearVariationalSolver` """ if hasattr(self, "_equidistributor"): return self._equidistributor - assert hasattr(self, "phi") - assert hasattr(self, "sigma") n = ufl.FacetNormal(self.mesh) sigma = firedrake.TrialFunction(self.P1_ten) tau = firedrake.TestFunction(self.P1_ten) @@ -338,31 +352,22 @@ def equidistributor(self): + ufl.dot(ufl.dot(tangential(ufl.grad(self.phi), n), tau), n) * self.ds ) problem = firedrake.LinearVariationalProblem(a, L, self.sigma) - sp = { - "ksp_type": "cg", - "pc_type": "bjacobi", - "sub_pc_type": "ilu", - } self._equidistributor = firedrake.LinearVariationalSolver( - problem, solver_parameters=sp + problem, solver_parameters=solver_parameters.cg_ilu ) return self._equidistributor - @PETSc.Log.EventDecorator("MongeAmpereMover.move") + @PETSc.Log.EventDecorator() def move(self): """ Run the relaxation method to convergence and update the mesh. + + :return: the iteration count + :rtype: :class:`int` """ - assert hasattr(self, "phi") - assert hasattr(self, "sigma") - assert hasattr(self, "phi_old") - assert hasattr(self, "sigma_old") for i in range(self.maxiter): - # L2 project self.l2_projector.solve() - - # Update mesh coordinates - self.mesh.coordinates.assign(self.x) + self._update_coordinates() # Update monitor function self.monitor.interpolate(self.monitor_function(self.mesh)) @@ -383,14 +388,17 @@ def move(self): f" Residual {residual:10.4e}" f" Variation (σ/μ) {cv:10.4e}" ) + plural = "s" if i != 0 else "" if residual < self.rtol: - PETSc.Sys.Print(f"Converged in {i+1} iterations.") + PETSc.Sys.Print(f"Converged in {i+1} iteration{plural}.") break if residual > self.dtol * initial_norm: - raise firedrake.ConvergenceError(f"Diverged after {i+1} iterations.") + raise firedrake.ConvergenceError( + f"Diverged after {i+1} iteration{plural}." + ) if i == self.maxiter - 1: raise firedrake.ConvergenceError( - f"Failed to converge in {i+1} iterations." + f"Failed to converge in {i+1} iteration{plural}." ) # Apply pseudotimestepper and equidistributor @@ -398,43 +406,41 @@ def move(self): self.equidistributor.solve() self.phi_old.assign(self.phi) self.sigma_old.assign(self.sigma) - self.mesh.coordinates.assign(self.x) + self._update_coordinates() return i class MongeAmpereMover_QuasiNewton(MongeAmpereMover_Base): r""" - The elliptic Monge-Ampere equation is solved using a quasi-Newton - method (see McRae et al. 2018 for details). - - References - ========== - A. T. T. McRae, C. J. Cotter, C. J. Budd, Optimal-transport-based - mesh adaptivity on the plane and sphere using finite elements, SIAM - Journal on Scientific Computing 40 (2) (2018) 1121–1148. - doi:10.1137/16M1109515. + The elliptic Monge-Ampere equation is solved using a quasi-Newton method (see + :cite:`MCB:18` for details). """ - @PETSc.Log.EventDecorator("MongeAmpereMover.__init__") - def __init__(self, mesh, monitor_function, **kwargs): + @PETSc.Log.EventDecorator() + def __init__( + self, mesh, monitor_function, phi_init=None, sigma_init=None, **kwargs + ): """ :arg mesh: the physical mesh + :type mesh: :class:`firedrake.mesh.MeshGeometry` :arg monitor_function: a Python function which takes a mesh as input - :kwarg maxiter: maximum number of iterations for the relaxation + :type monitor_function: :class:`~.Callable` + :kwarg phi_init: initial guess for the scalar potential + :type phi_init: :class:`firedrake.function.Function` + :kwarg sigma_init: initial guess for the Hessian + :type sigma_init: :class:`firedrake.function.Function` + :kwarg maxiter: maximum number of iterations for the Quasi-Newton solver + :type maxiter: :class:`int` :kwarg rtol: relative tolerance for the residual + :type rtol: :class:`float` :kwarg dtol: divergence tolerance for the residual + :type dtol: :class:`float` """ super().__init__(mesh, monitor_function=monitor_function, **kwargs) - # Create functions to hold solution data - self.V = self.P1 * self.P1_ten - self.phisigma = firedrake.Function(self.V) - self.phi, self.sigma = self.phisigma.subfunctions - self.phisigma_old = firedrake.Function(self.V) - self.phi_old, self.sigma_old = self.phisigma_old.subfunctions - # Initialise phi and sigma - self.apply_initial_guess(**kwargs) + if phi_init or sigma_init: + self.apply_initial_guess(phi_init, sigma_init) # Setup residuals I = ufl.Identity(self.dim) @@ -444,15 +450,25 @@ def __init__(self, mesh, monitor_function, **kwargs): self._residual_l2_form = psi * self.residual * self.dx self._norm_l2_form = psi * self.theta * self.dx + def _create_functions(self): + super()._create_functions() + self.V = self.P1 * self.P1_ten + self.phisigma = firedrake.Function(self.V) + self.phi, self.sigma = self.phisigma.subfunctions + self.phisigma_old = firedrake.Function(self.V) + self.phi_old, self.sigma_old = self.phisigma_old.subfunctions + @property - @PETSc.Log.EventDecorator("MongeAmpereMover.create_equidistributor") + @PETSc.Log.EventDecorator() def equidistributor(self): """ Setup the equidistributor for the quasi-newton method. + + :return: the equidistributor + :rtype: :class:`~.NonlinearVariationalSolver` """ if hasattr(self, "_equidistributor"): return self._equidistributor - assert hasattr(self, "phisigma") n = ufl.FacetNormal(self.mesh) I = ufl.Identity(self.dim) phi, sigma = firedrake.split(self.phisigma) @@ -470,11 +486,10 @@ def update_monitor(cursol): """ Callback for updating the monitor function. """ - assert hasattr(self, "phisigma_old") with self.phisigma_old.dat.vec as v: cursol.copy(v) self.l2_projector.solve() - self.mesh.coordinates.assign(self.x) + self._update_coordinates() self.monitor.interpolate(self.monitor_function(self.mesh)) self.mesh.coordinates.assign(self.xi) self.theta.assign( @@ -519,7 +534,7 @@ def monitor(snes, i, rnorm): """ cursol = snes.getSolution() update_monitor(cursol) - self.mesh.coordinates.assign(self.x) + self._update_coordinates() firedrake.assemble(self.L_P0, tensor=self.volume) self.volume.interpolate(self.volume / self.original_volume) self.mesh.coordinates.assign(self.xi) @@ -535,52 +550,24 @@ def monitor(snes, i, rnorm): self.snes.setMonitor(monitor) return self._equidistributor - @PETSc.Log.EventDecorator("MongeAmpereMover.move") + @PETSc.Log.EventDecorator() def move(self): """ Run the quasi-Newton method to convergence and update the mesh. + + :return: the iteration count + :rtype: :class:`int` """ try: self.equidistributor.solve() i = self.snes.getIterationNumber() - PETSc.Sys.Print(f"Converged in {i} iterations.") + plural = "s" if i != 1 else "" + PETSc.Sys.Print(f"Converged in {i} iteration{plural}.") except firedrake.ConvergenceError: i = self.snes.getIterationNumber() - raise firedrake.ConvergenceError(f"Failed to converge in {i} iterations.") - self.mesh.coordinates.assign(self.x) + plural = "s" if i != 1 else "" + raise firedrake.ConvergenceError( + f"Failed to converge in {i} iteration{plural}." + ) + self._update_coordinates() return i - - -def monge_ampere(mesh, monitor_function, method="relaxation", **kwargs): - r""" - Movement of a `mesh` is determined by a `monitor_function` - :math:`m` and the Monge-Ampère type equation - - .. math:: - m(x)\det(I + H(\phi)) = \theta, - - for a scalar potential :math:`\phi`, where :math:`I` is the - identity matrix, :math:`\theta` is a normalisation coefficient - and :math:`H(\phi)` denotes the Hessian of :math:`\phi` with - respect to the coordinates :math:`\xi` of the computational mesh. - - The physical mesh coordinates :math:`x` are updated according to - - .. math:: - x = \xi + \nabla\phi. - - :arg mesh: the physical mesh - :arg monitor_function: a Python function which takes a mesh as input - :kwarg method: choose from 'relaxation' and 'quasi_newton' - :kwarg phi_init: initial guess for the scalar potential - :kwarg sigma_init: initial guess for the Hessian - :return: converged scalar potential and Hessian - """ - if method == "relaxation": - mover = MongeAmpereMover_Relaxation(mesh, monitor_function, **kwargs) - elif method == "quasi_newton": - mover = MongeAmpereMover_QuasiNewton(mesh, monitor_function, **kwargs) - else: - raise ValueError(f"Monge-Ampere solver {method} not recognised.") - mover.adapt() - return mover.phi, mover.sigma diff --git a/movement/mover.py b/movement/mover.py index 16c443c..cbd1ab9 100644 --- a/movement/mover.py +++ b/movement/mover.py @@ -27,7 +27,7 @@ def __init__(self, mesh, monitor_function=None, **kwargs): # Mesh coordinate functions self.coord_space = self.mesh.coordinates.function_space() - self._x = firedrake.Function(self.mesh.coordinates, name="Physical coordinates") + self.x = firedrake.Function(self.mesh.coordinates, name="Physical coordinates") self.xi = firedrake.Function( self.mesh.coordinates, name="Computational coordinates" ) diff --git a/movement/solver_parameters.py b/movement/solver_parameters.py index 4648ff0..b99379e 100644 --- a/movement/solver_parameters.py +++ b/movement/solver_parameters.py @@ -2,7 +2,7 @@ "ksp_type": "gmres", "pc_type": "fieldsplit", "pc_fieldsplit_type": "multiplicative", - "pc_fieldsplit_off_diag_use_amat": True, + "pc_fieldsplit_off_diag_use_amat": None, "fieldsplit_0_pc_type": "gamg", "fieldsplit_0_ksp_type": "preonly", "fieldsplit_0_mg_levels_ksp_max_it": 5, @@ -24,7 +24,7 @@ "ksp_type": "gmres", "pc_type": "fieldsplit", "pc_fieldsplit_type": "multiplicative", - "pc_fieldsplit_off_diag_use_amat": True, + "pc_fieldsplit_off_diag_use_amat": None, "fieldsplit_0_pc_type": "gamg", "fieldsplit_0_ksp_type": "preonly", "fieldsplit_0_mg_levels_ksp_max_it": 5, @@ -59,12 +59,17 @@ "pc_type": "jacobi", } -cg = { +cg_ilu = { "ksp_type": "cg", "pc_type": "bjacobi", "pc_sub_type": "ilu", } +cg_gamg = { + "ksp_type": "cg", + "pc_type": "gamg", +} + lu = { "mat_type": "aij", "snes_type": "ksponly", diff --git a/test/test_monge_ampere.py b/test/test_monge_ampere.py index ecdfa79..6f3506f 100644 --- a/test/test_monge_ampere.py +++ b/test/test_monge_ampere.py @@ -1,133 +1,197 @@ +import unittest +from unittest.mock import MagicMock + import numpy as np -import pytest from monitors import * +from parameterized import parameterized from movement import * -@pytest.fixture(params=["relaxation", "quasi_newton"]) -def method(request): - return request.param - - -@pytest.fixture(params=[True, False]) -def fix_boundary(request): - return request.param - - -@pytest.fixture(params=[2, 3]) -def dimension(request): - return request.param - - -def uniform_mesh(dim, n): - if dim == 2: - return UnitSquareMesh(n, n) - else: - return UnitCubeMesh(n, n, n) - - -def test_uniform_monitor(dimension, method, exports=False): +class TestMongeAmpere(unittest.TestCase): """ - Test that the mesh mover converges in one - iteration for a constant monitor function. + Unit tests for Monge-Ampere methods. """ - n = 10 - mesh = uniform_mesh(dimension, n) - coords = mesh.coordinates.dat.data.copy() - mover = MongeAmpereMover(mesh, const_monitor, method=method) - num_iterations = mover.move() + def mesh(self, dim, n=10): + self.assertTrue(dim in (2, 3)) + return UnitSquareMesh(n, n) if dim == 2 else UnitCubeMesh(n, n, n) - assert np.allclose(coords, mover.mesh.coordinates.dat.data) - assert num_iterations == 0 + @property + def dummy_mesh(self): + return MagicMock(UnitSquareMesh) + @property + def dummy_monitor(self): + return lambda *args: MagicMock(Function) -def test_continue(dimension, method, exports=False): - """ - Test that providing a good initial guess - benefits the solver. - """ - n = 20 - mesh = uniform_mesh(dimension, n) - rtol = 1.0e-03 - - # Solve the problem to a weak tolerance - mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=0.1) - num_it_init = mover.move() - if exports: - File("outputs/init.pvd").write(mover.phi, mover.sigma) - - # Continue with a tighter tolerance - mover.rtol = rtol - num_it_continue = mover.move() - if exports: - File("outputs/continue.pvd").write(mover.phi, mover.sigma) - - # Solve the problem again to a tight tolerance - mesh = uniform_mesh(dimension, n) - mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=rtol) - num_it_naive = mover.move() - if exports: - File("outputs/naive.pvd").write(mover.phi, mover.sigma) - - assert num_it_continue <= num_it_naive - assert num_it_init + num_it_continue <= num_it_naive - # FIXME: Looks like the mesh is tangled or close to tangling - # for the relaxation method, which is concerning. - - -def test_change_monitor(dimension, method, exports=False): - """ - Test that the mover can handle changes to - the monitor function, such as would happen - during timestepping. - """ - n = 20 - mesh = uniform_mesh(dimension, n) - dim = mesh.geometric_dimension() - coords = mesh.coordinates.dat.data.copy() - tol = 1.0e-03 - - # Adapt to a ring monitor - mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=1e3 * tol**dim) - mover.move() - if exports: - File("outputs/ring.pvd").write(mover.phi, mover.sigma) - assert not np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol) - - # Adapt to a constant monitor - mover.monitor_function = const_monitor - mover.move() - if exports: - File("outputs/const.pvd").write(mover.phi, mover.sigma) - assert np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol) - - -@pytest.mark.slow -def test_bcs(dimension, method, fix_boundary): - """ - Test that domain boundaries are fixed by - the Monge-Ampere movers. - """ - n = 20 - mesh = uniform_mesh(dimension, n) - one = Constant(1.0) - bnd = assemble(one * ds(domain=mesh)) - bnodes = DirichletBC(mesh.coordinates.function_space(), 0, "on_boundary").nodes - bnd_coords = mesh.coordinates.dat.data.copy()[bnodes] - - # Adapt to a ring monitor - mover = MongeAmpereMover( - mesh, ring_monitor, method=method, fix_boundary_nodes=fix_boundary - ) - mover.move() + def test_method_valueerror(self): + with self.assertRaises(ValueError) as cm: + MongeAmpereMover(self.dummy_mesh, self.dummy_monitor, method="method") + self.assertEqual(str(cm.exception), "Method 'method' not recognised.") - # Check boundary lengths are preserved - bnd_new = assemble(one * ds(domain=mover.mesh)) - assert np.isclose(bnd, bnd_new) + def test_no_monitor_valueerror(self): + with self.assertRaises(ValueError) as cm: + MongeAmpereMover(self.dummy_mesh, None) + self.assertEqual(str(cm.exception), "Please supply a monitor function.") - # Check boundaries are indeed fixed - if fix_boundary: - bnd_coords_new = mover.mesh.coordinates.dat.data[bnodes] - assert np.allclose(bnd_coords, bnd_coords_new) + @parameterized.expand( + [(2, "relaxation"), (2, "quasi_newton"), (3, "relaxation"), (3, "quasi_newton")] + ) + def test_uniform_monitor(self, dim, method): + """ + Test that the mesh mover converges in one iteration for a constant monitor + function. + """ + mesh = self.mesh(dim) + coords = mesh.coordinates.dat.data.copy() + P1 = FunctionSpace(mesh, "CG", 1) + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + + mover = MongeAmpereMover( + mesh, + const_monitor, + method=method, + phi_init=Function(P1), + sigma_init=Function(P1_ten), + rtol=1e-3, + ) + num_iterations = mover.move() + + self.assertTrue(np.allclose(coords, mover.mesh.coordinates.dat.data)) + self.assertEqual(num_iterations, 0) + + @parameterized.expand( + [(2, "relaxation"), (2, "quasi_newton"), (3, "relaxation"), (3, "quasi_newton")] + ) + def test_continue(self, dim, method): + """ + Test that providing a good initial guess benefits the solver. + """ + mesh = self.mesh(dim) + rtol = 1.0e-03 + + # Solve the problem to a weak tolerance + mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=0.1) + num_it_init = mover.move() + + # Continue with a tighter tolerance + mover.rtol = rtol + num_it_continue = mover.move() + + # Solve the problem again to a tight tolerance + mesh = self.mesh(dim) + mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=rtol) + num_it_naive = mover.move() + + self.assertLessEqual(num_it_continue, num_it_naive) + self.assertLessEqual(num_it_init + num_it_continue, num_it_naive) + # FIXME: Looks like the mesh is tangled or close to tangling + # for the relaxation method, which is concerning. + + @parameterized.expand( + [(2, "relaxation"), (2, "quasi_newton"), (3, "relaxation"), (3, "quasi_newton")] + ) + def test_change_monitor(self, dim, method): + """ + Test that the mover can handle changes to the monitor function, such as would + happen during timestepping. + """ + mesh = self.mesh(dim) + gdim = mesh.geometric_dimension() + coords = mesh.coordinates.dat.data.copy() + atol = 1.0e-03 + rtol = 1.0e02 * atol**gdim + + # Adapt to a ring monitor + mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=rtol) + mover.move() + moved_coords = mover.mesh.coordinates.dat.data + self.assertFalse(np.allclose(coords, moved_coords, atol=atol)) + + # Adapt to a constant monitor + mover.monitor_function = const_monitor + mover.move() + moved_coords = mover.mesh.coordinates.dat.data + self.assertTrue(np.allclose(coords, moved_coords, atol=atol)) + + @parameterized.expand( + [ + (2, "relaxation", True), + (2, "relaxation", False), + (2, "quasi_newton", True), + (2, "quasi_newton", False), + (3, "relaxation", True), + (3, "relaxation", False), + (3, "quasi_newton", True), + (3, "quasi_newton", False), + ] + ) + def test_bcs(self, dim, method, fix_boundary): + """ + Test that domain boundaries are fixed by the Monge-Ampere movers. + """ + mesh = self.mesh(dim) + bnd = assemble(Constant(1.0) * ds(domain=mesh)) + bnodes = DirichletBC(mesh.coordinates.function_space(), 0, "on_boundary").nodes + bnd_coords = mesh.coordinates.dat.data.copy()[bnodes] + + # Adapt to a ring monitor + mover = MongeAmpereMover( + mesh, + ring_monitor, + method=method, + fix_boundary_nodes=fix_boundary, + rtol=1e-3, + ) + mover.move() + + # Check boundary lengths are preserved + bnd_new = assemble(Constant(1.0) * ds(domain=mover.mesh)) + self.assertAlmostEqual(bnd, bnd_new) + + # Check boundaries are indeed fixed + if fix_boundary: + bnd_coords_new = mover.mesh.coordinates.dat.data[bnodes] + self.assertTrue(np.allclose(bnd_coords, bnd_coords_new)) + + @parameterized.expand([("relaxation"), ("quasi_newton")]) + def test_maxiter_convergenceerror(self, method): + """ + Test that the mesh mover raises a :class:`~.ConvergenceError` if it reaches the + maximum number of iterations. + """ + mesh = self.mesh(2, n=4) + mover = MongeAmpereMover(mesh, ring_monitor, method=method, maxiter=1) + with self.assertRaises(ConvergenceError) as cm: + mover.move() + self.assertEqual(str(cm.exception), "Failed to converge in 1 iteration.") + + def test_divergence_convergenceerror(self): + """ + Test that the mesh mover raises a :class:`~.ConvergenceError` if it diverges. + """ + mesh = self.mesh(2, n=4) + mover = MongeAmpereMover_Relaxation(mesh, ring_monitor, dtol=1.0e-08) + with self.assertRaises(ConvergenceError) as cm: + mover.move() + self.assertEqual(str(cm.exception), "Diverged after 1 iteration.") + + def test_initial_guess_valueerror(self): + mesh = self.mesh(2, n=2) + phi_init = Function(FunctionSpace(mesh, "CG", 1)) + with self.assertRaises(ValueError) as cm: + MongeAmpereMover_Relaxation(mesh, ring_monitor, phi_init=phi_init) + self.assertEqual(str(cm.exception), "Need to initialise both phi *and* sigma.") + + def test_coordinate_update(self): + mesh = self.mesh(2, n=2) + dg_coords = Function(VectorFunctionSpace(mesh, "DG", 1)) + dg_coords.project(mesh.coordinates) + mover = MongeAmpereMover(Mesh(dg_coords), const_monitor) + mover._grad_phi.assign(2.0) + mover._update_coordinates() + self.assertNotEqual(mover.grad_phi.dof_dset.size, mover._grad_phi.dof_dset.size) + self.assertAlmostEqual(errornorm(mover.grad_phi, mover._grad_phi), 0) + mover.xi += 2 + self.assertAlmostEqual(errornorm(mover.x, mover.xi), 0)