Merge branch 'main' into 101_tangling_check

demos/laplacian_smoothing.py

@@ -141,8 +141,6 @@ def update_boundary_velocity(t):
         update_boundary_velocity=update_boundary_velocity,
         boundary_conditions=boundary_conditions,
     )
-    displacement = np.linalg.norm(mover.displacement)
-    print(f"time = {time:.1f} s, displacement = {displacement:.2f} m")
 
     # Plot the current mesh, adding a time label
     ax = axes[idx // 4, idx % 4]
@@ -152,6 +150,22 @@ def update_boundary_velocity(t):
 axes[0, 1].axis(False)
 plt.savefig("laplacian_smoothing-adapted_meshes.jpg")
 
+# This should give command line output similar to the following:
+#
+# .. code-block:: none
+#
+#    0.00 s   Volume ratio  1.00   Variation (σ/μ) 5.16e-16   Displacement 0.00 m
+#    0.10 s   Volume ratio  1.03   Variation (σ/μ) 7.35e-03   Displacement 0.04 m
+#    0.20 s   Volume ratio  1.08   Variation (σ/μ) 1.90e-02   Displacement 0.06 m
+#    0.30 s   Volume ratio  1.13   Variation (σ/μ) 3.04e-02   Displacement 0.06 m
+#    0.40 s   Volume ratio  1.17   Variation (σ/μ) 3.73e-02   Displacement 0.04 m
+#    0.50 s   Volume ratio  1.17   Variation (σ/μ) 3.73e-02   Displacement 0.00 m
+#    0.60 s   Volume ratio  1.13   Variation (σ/μ) 3.04e-02   Displacement 0.04 m
+#    0.70 s   Volume ratio  1.08   Variation (σ/μ) 1.90e-02   Displacement 0.06 m
+#    0.80 s   Volume ratio  1.03   Variation (σ/μ) 7.35e-03   Displacement 0.06 m
+#    0.90 s   Volume ratio  1.00   Variation (σ/μ) 1.34e-10   Displacement 0.04 m
+#    1.00 s   Volume ratio  1.00   Variation (σ/μ) 1.34e-10   Displacement 0.00 m
+#
 # .. figure:: laplacian_smoothing-adapted_meshes.jpg
 #    :figwidth: 100%
 #    :align: center

demos/lineal_spring.py

@@ -158,8 +158,6 @@ def update_boundary_displacement(t):
         update_boundary_displacement=update_boundary_displacement,
         boundary_conditions=boundary_conditions,
     )
-    displacement = np.linalg.norm(mover.displacement)
-    print(f"time = {time:.1f} s, displacement = {displacement:.2f} m")
 
     # Plot the current mesh, adding a time label
     ax = axes[idx // 4, idx % 4]
@@ -169,6 +167,22 @@ def update_boundary_displacement(t):
 axes[0, 1].axis(False)
 plt.savefig("lineal_spring-adapted_meshes.jpg")
 
+# This should give command line output similar to the following:
+#
+# .. code-block:: none
+#
+#    0.00   Volume ratio  1.00   Variation (σ/μ) 5.16e-16   Displacement 0.00 m
+#    0.10   Volume ratio  1.02   Variation (σ/μ) 4.23e-03   Displacement 0.05 m
+#    0.20   Volume ratio  1.05   Variation (σ/μ) 1.11e-02   Displacement 0.08 m
+#    0.30   Volume ratio  1.09   Variation (σ/μ) 1.82e-02   Displacement 0.08 m
+#    0.40   Volume ratio  1.11   Variation (σ/μ) 2.27e-02   Displacement 0.05 m
+#    0.50   Volume ratio  1.11   Variation (σ/μ) 2.27e-02   Displacement 0.00 m
+#    0.60   Volume ratio  1.09   Variation (σ/μ) 1.81e-02   Displacement 0.05 m
+#    0.70   Volume ratio  1.05   Variation (σ/μ) 1.08e-02   Displacement 0.08 m
+#    0.80   Volume ratio  1.02   Variation (σ/μ) 3.82e-03   Displacement 0.08 m
+#    0.90   Volume ratio  1.01   Variation (σ/μ) 1.32e-03   Displacement 0.05 m
+#    1.00   Volume ratio  1.01   Variation (σ/μ) 1.32e-03   Displacement 0.00 m
+#
 # .. figure:: lineal_spring-adapted_meshes.jpg
 #    :figwidth: 100%
 #    :align: center

demos/monge_ampere1.py

@@ -102,6 +102,65 @@ def ring_monitor(mesh):
 mover = MongeAmpereMover(mesh, ring_monitor, method="quasi_newton", rtol=rtol)
 mover.move()
 
+# This should give command line output similar to the following:
+#
+# .. code-block:: none
+#
+#       0   Volume ratio 11.49   Variation (σ/μ) 9.71e-01   Residual 9.39e-01
+#       1   Volume ratio  8.32   Variation (σ/μ) 6.84e-01   Residual 5.35e-01
+#       2   Volume ratio  5.74   Variation (σ/μ) 5.55e-01   Residual 3.83e-01
+#       3   Volume ratio  6.86   Variation (σ/μ) 4.92e-01   Residual 3.06e-01
+#       4   Volume ratio  5.91   Variation (σ/μ) 4.53e-01   Residual 2.69e-01
+#       5   Volume ratio  8.38   Variation (σ/μ) 4.20e-01   Residual 2.22e-01
+#       6   Volume ratio  7.34   Variation (σ/μ) 4.12e-01   Residual 2.14e-01
+#       7   Volume ratio  7.68   Variation (σ/μ) 4.02e-01   Residual 2.03e-01
+#       8   Volume ratio  7.93   Variation (σ/μ) 3.84e-01   Residual 1.84e-01
+#       9   Volume ratio  7.81   Variation (σ/μ) 3.83e-01   Residual 1.86e-01
+#      10   Volume ratio  7.60   Variation (σ/μ) 3.93e-01   Residual 1.97e-01
+#      11   Volume ratio  7.99   Variation (σ/μ) 4.14e-01   Residual 2.13e-01
+#      12   Volume ratio  8.22   Variation (σ/μ) 4.21e-01   Residual 2.20e-01
+#      13   Volume ratio 10.79   Variation (σ/μ) 4.54e-01   Residual 2.13e-01
+#      14   Volume ratio  9.66   Variation (σ/μ) 4.15e-01   Residual 1.33e-01
+#      15   Volume ratio 10.52   Variation (σ/μ) 3.75e-01   Residual 9.77e-02
+#      16   Volume ratio 10.00   Variation (σ/μ) 3.90e-01   Residual 8.64e-02
+#      17   Volume ratio  9.00   Variation (σ/μ) 3.61e-01   Residual 6.33e-02
+#      18   Volume ratio  9.53   Variation (σ/μ) 3.73e-01   Residual 4.41e-02
+#      19   Volume ratio  8.86   Variation (σ/μ) 3.60e-01   Residual 3.71e-02
+#      20   Volume ratio  9.38   Variation (σ/μ) 3.65e-01   Residual 2.71e-02
+#      21   Volume ratio  8.95   Variation (σ/μ) 3.57e-01   Residual 2.23e-02
+#      22   Volume ratio  9.15   Variation (σ/μ) 3.57e-01   Residual 1.32e-02
+#      23   Volume ratio  8.90   Variation (σ/μ) 3.52e-01   Residual 8.93e-03
+#      24   Volume ratio  8.87   Variation (σ/μ) 3.50e-01   Residual 3.93e-03
+#      25   Volume ratio  8.80   Variation (σ/μ) 3.48e-01   Residual 2.61e-03
+#      26   Volume ratio  8.85   Variation (σ/μ) 3.49e-01   Residual 1.51e-03
+#      27   Volume ratio  8.83   Variation (σ/μ) 3.48e-01   Residual 1.15e-03
+#      28   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 7.98e-04
+#      29   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 6.27e-04
+#      30   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 4.46e-04
+#      31   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 3.46e-04
+#      32   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 2.39e-04
+#      33   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 1.77e-04
+#      34   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 1.14e-04
+#      35   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 7.82e-05
+#      36   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 4.69e-05
+#      37   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 2.96e-05
+#      38   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 1.77e-05
+#      39   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 1.11e-05
+#      40   Volume ratio  8.85   Variation (σ/μ) 3.48e-01   Residual 7.43e-06
+#      41   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 5.07e-06
+#      42   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 3.86e-06
+#      43   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 2.85e-06
+#      44   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 2.30e-06
+#      45   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 1.72e-06
+#      46   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 1.38e-06
+#      47   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 9.75e-07
+#      48   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 7.42e-07
+#      49   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 4.50e-07
+#      50   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 3.00e-07
+#      51   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 1.42e-07
+#      52   Volume ratio  8.84   Variation (σ/μ) 3.48e-01   Residual 7.93e-08
+#    Solver converged in 52 iterations.
+#
 # The adapted mesh can be accessed via the `mesh` attribute of the mover. Plotting it,
 # we see that the adapted mesh has its resolution focused around a ring, as expected.
 

demos/monge_ampere_helmholtz.py

@@ -142,30 +142,30 @@ def monitor(mesh):
 mover.move()
 
 # For every iteration the MongeAmpereMover prints the minimum to maximum ratio of
-# the cell areas in the mesh, the residual in the Monge Ampere equation, and the
+# the cell areas in the mesh, the residual in the Monge-Ampère equation, and the
 # coefficient of variation of the cell areas:
 #
 # .. code-block:: none
 #
-#    0   Min/Max 2.0268e-01   Residual 4.7659e-01   Variation (σ/μ) 9.9384e-01
-#    1   Min/Max 3.7852e-01   Residual 2.4133e-01   Variation (σ/μ) 9.9659e-01
-#    2   Min/Max 5.9791e-01   Residual 1.2442e-01   Variation (σ/μ) 9.9774e-01
-#    3   Min/Max 7.1000e-01   Residual 6.5811e-02   Variation (σ/μ) 9.9804e-01
-#    4   Min/Max 7.7704e-01   Residual 3.4929e-02   Variation (σ/μ) 9.9818e-01
-#    5   Min/Max 8.3434e-01   Residual 1.7261e-02   Variation (σ/μ) 9.9829e-01
-#    6   Min/Max 8.5805e-01   Residual 7.7528e-03   Variation (σ/μ) 9.9833e-01
-#    7   Min/Max 8.6653e-01   Residual 3.1551e-03   Variation (σ/μ) 9.9835e-01
-#    8   Min/Max 8.6796e-01   Residual 1.1644e-03   Variation (σ/μ) 9.9835e-01
-#    9   Min/Max 8.6792e-01   Residual 3.8816e-04   Variation (σ/μ) 9.9835e-01
-#   10   Min/Max 8.6784e-01   Residual 1.1574e-04   Variation (σ/μ) 9.9835e-01
-#   11   Min/Max 8.6778e-01   Residual 1.5645e-05   Variation (σ/μ) 9.9835e-01
-#   12   Min/Max 8.6776e-01   Residual 7.5654e-06   Variation (σ/μ) 9.9835e-01
-#   13   Min/Max 8.6776e-01   Residual 3.5803e-06   Variation (σ/μ) 9.9835e-01
-#   14   Min/Max 8.6775e-01   Residual 1.5113e-06   Variation (σ/μ) 9.9835e-01
-#   15   Min/Max 8.6775e-01   Residual 5.7080e-07   Variation (σ/μ) 9.9835e-01
-#   16   Min/Max 8.6775e-01   Residual 1.9357e-07   Variation (σ/μ) 9.9835e-01
-#   17   Min/Max 8.6775e-01   Residual 5.8585e-08   Variation (σ/μ) 9.9835e-01
-#   Converged in 17 iterations.
+#       0   Volume ratio  4.93   Variation (σ/μ) 4.73e-01   Residual 4.77e-01
+#       1   Volume ratio  2.64   Variation (σ/μ) 2.44e-01   Residual 2.41e-01
+#       2   Volume ratio  1.67   Variation (σ/μ) 1.31e-01   Residual 1.24e-01
+#       3   Volume ratio  1.41   Variation (σ/μ) 7.57e-02   Residual 6.58e-02
+#       4   Volume ratio  1.29   Variation (σ/μ) 4.77e-02   Residual 3.49e-02
+#       5   Volume ratio  1.20   Variation (σ/μ) 3.37e-02   Residual 1.73e-02
+#       6   Volume ratio  1.17   Variation (σ/μ) 2.81e-02   Residual 7.75e-03
+#       7   Volume ratio  1.15   Variation (σ/μ) 2.64e-02   Residual 3.16e-03
+#       8   Volume ratio  1.15   Variation (σ/μ) 2.59e-02   Residual 1.16e-03
+#       9   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 3.88e-04
+#      10   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 1.16e-04
+#      11   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 1.56e-05
+#      12   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 7.57e-06
+#      13   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 3.58e-06
+#      14   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 1.51e-06
+#      15   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 5.71e-07
+#      16   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 1.94e-07
+#      17   Volume ratio  1.15   Variation (σ/μ) 2.57e-02   Residual 5.86e-08
+#    Solver converged in 17 iterations.
 #
 # Plotting the resulting mesh:
 

movement/laplacian_smoothing.py

@@ -45,15 +45,18 @@ def __init__(self, mesh, timestep, **kwargs):
         dim = self.mesh.topological_dimension()
         self.displacement = np.zeros((self.mesh.num_vertices(), dim))
 
+    def _create_functions(self):
+        super()._create_functions()
+        self.rhs = firedrake.Function(self.coord_space, name="Zero RHS")
+
     @PETSc.Log.EventDecorator()
     def _setup_solver(self, boundary_conditions):
         if not hasattr(self, "_solver"):
             test = firedrake.TestFunction(self.coord_space)
             trial = firedrake.TrialFunction(self.coord_space)
-            f = firedrake.Function(self.coord_space, name="Zero RHS")
 
             a = ufl.inner(ufl.grad(trial), ufl.grad(test)) * self.dx
-            L = ufl.inner(f, test) * self.dx
+            L = ufl.inner(self.rhs, test) * self.dx
             problem = firedrake.LinearVariationalProblem(
                 a, L, self.v, bcs=boundary_conditions
             )
@@ -96,6 +99,12 @@ 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"
+            f"   Volume ratio {self.volume_ratio:5.2f}"
+            f"   Variation (σ/μ) {self.coefficient_of_variation:8.2e}"
+            f"   Displacement {np.linalg.norm(self.displacement):.2f} m"
+        )

movement/monge_ampere.py

@@ -97,28 +97,24 @@ def __init__(self, mesh, monitor_function, **kwargs):
         self.dtol = kwargs.pop("dtol", 2.0)
         self.fix_boundary_nodes = kwargs.pop("fix_boundary_nodes", False)
         super().__init__(mesh, monitor_function=monitor_function, **kwargs)
+        self.theta = firedrake.Constant(0.0)
 
-        # Create function spaces
-        self.P0 = firedrake.FunctionSpace(self.mesh, "DG", 0)
+    def _create_function_spaces(self):
+        super()._create_function_spaces()
         self.P1 = firedrake.FunctionSpace(self.mesh, "CG", 1)
         self.P1_vec = firedrake.VectorFunctionSpace(self.mesh, "CG", 1)
         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.interpolate(self.monitor_function(self.mesh))
-        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):
+        super()._create_functions()
         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)
+        self.monitor.interpolate(self.monitor_function(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
 
     @PETSc.Log.EventDecorator()
     def apply_initial_guess(self, phi_init, sigma_init):
@@ -143,27 +139,17 @@ def apply_initial_guess(self, phi_init, sigma_init):
             raise ValueError("Need to initialise both phi *and* sigma.")
 
     @property
-    @PETSc.Log.EventDecorator()
-    def _diagnostics(self):
+    def relative_l2_residual(self):
         """
-        Compute the following diagnostics:
-          1) the ratio of the smallest and largest element volumes;
-          2) coefficient of variation (σ/μ) of element volumes;
-          3) relative L2 norm residual.
+        :return: the relative :math:`L^2` norm residual.
+        :rtype: :class:`float`
         """
-        v = self.volume.vector().gather()
-        minmax = v.min() / v.max()
-        mean = v.sum() / v.max()
-        w = v.copy() - mean
-        w *= w
-        std = np.sqrt(w.sum() / w.size)
-        cv = std / mean
         assert hasattr(self, "_residual_l2_form")
         assert hasattr(self, "_norm_l2_form")
-        residual_l2 = firedrake.assemble(self._residual_l2_form).dat.norm
-        norm_l2 = firedrake.assemble(self._norm_l2_form).dat.norm
-        residual_l2_rel = residual_l2 / norm_l2
-        return minmax, residual_l2_rel, cv
+        return (
+            firedrake.assemble(self._residual_l2_form).dat.norm
+            / firedrake.assemble(self._norm_l2_form).dat.norm
+        )
 
     @PETSc.Log.EventDecorator()
     def _update_coordinates(self):
@@ -385,14 +371,14 @@ def move(self):
             self.theta.assign(firedrake.assemble(self.theta_form) / self.total_volume)
 
             # Check convergence criteria
-            minmax, residual, cv = self._diagnostics
+            residual = self.relative_l2_residual
             if i == 0:
                 initial_norm = residual
             PETSc.Sys.Print(
                 f"{i:4d}"
-                f"   Min/Max {minmax:10.4e}"
-                f"   Residual {residual:10.4e}"
-                f"   Variation (σ/μ) {cv:10.4e}"
+                f"   Volume ratio {self.volume_ratio:5.2f}"
+                f"   Variation (σ/μ) {self.coefficient_of_variation:8.2e}"
+                f"   Residual {residual:8.2e}"
             )
             if residual < self.rtol:
                 self._convergence_message(i + 1)
@@ -541,12 +527,11 @@ def monitor(snes, i, rnorm):
             firedrake.assemble(self.L_P0, tensor=self.volume)
             self.volume.interpolate(self.volume / self.original_volume)
             self.mesh.coordinates.assign(self.xi)
-            minmax, residual, cv = self._diagnostics
             PETSc.Sys.Print(
                 f"{i:4d}"
-                f"   Min/Max {minmax:10.4e}"
-                f"   Residual {residual:10.4e}"
-                f"   Variation (σ/μ) {cv:10.4e}"
+                f"   Volume ratio {self.volume_ratio:5.2f}"
+                f"   Variation (σ/μ) {self.coefficient_of_variation:8.2e}"
+                f"   Residual {self.relative_l2_residual:8.2e}"
             )
 
         self.snes = self._equidistributor.snes

movement/mover.py

@@ -1,8 +1,10 @@
+import abc
 from warnings import warn
 
 import firedrake
 import firedrake.exceptions as fexc
 import numpy as np
+import ufl
 from firedrake.cython.dmcommon import create_section
 from firedrake.petsc import PETSc
 
@@ -57,13 +59,23 @@ def __init__(
         self.ds = firedrake.ds(domain=self.mesh, degree=degree)
         self.dS = firedrake.dS(domain=self.mesh, degree=degree)
 
-        # Mesh coordinate functions
+        self._create_function_spaces()
+        self._create_functions()
+
+    @abc.abstractmethod
+    def _create_function_spaces(self):
         self.coord_space = self.mesh.coordinates.function_space()
+        self.P0 = firedrake.FunctionSpace(self.mesh, "DG", 0)
+
+    @abc.abstractmethod
+    def _create_functions(self):
         self.x = firedrake.Function(self.mesh.coordinates, name="Physical coordinates")
         self.xi = firedrake.Function(
             self.mesh.coordinates, name="Computational coordinates"
         )
         self.v = firedrake.Function(self.coord_space, name="Mesh velocity")
+        self.volume = firedrake.Function(self.P0, name="Mesh volume")
+        self.volume.interpolate(ufl.CellVolume(self.mesh))
 
         # Utilities
         if tangling_check is None:
@@ -181,6 +193,25 @@ def coordinate(self, index):
         """
         return self.mesh.coordinates.dat.data_with_halos[self.get_offset(index)]
 
+    @property
+    def volume_ratio(self):
+        """
+        :return: the ratio of the largest and smallest element volumes.
+        :rtype: :class:`float`
+        """
+        volume_array = self.volume.vector().gather()
+        return volume_array.max() / volume_array.min()
+
+    @property
+    def coefficient_of_variation(self):
+        """
+        :return: the coefficient of variation (σ/μ) of element volumes.
+        :rtype: :class:`float`
+        """
+        volume_array = self.volume.vector().gather()
+        mean = volume_array.sum() / volume_array.size
+        return np.sqrt(np.sum((volume_array - mean) ** 2) / volume_array.size) / mean
+
     def move(self):
         """
         Move the mesh according to the method of choice.