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

@@ -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."
 

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,

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

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)

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)