diff --git a/demos/monge_ampere_ring.py b/demos/monge_ampere_ring.py index f71fe00..fb68071 100644 --- a/demos/monge_ampere_ring.py +++ b/demos/monge_ampere_ring.py @@ -73,21 +73,16 @@ # # for some values of the parameters :math:`\alpha`, :math:`\beta`, and :math:`\gamma`. # Unity is added at the start to ensure that the monitor function doesn't get too -# close to zero. +# close to zero. Here we can think of :math:`\alpha` as relating to the amplitude of the +# monitor function, :math:`\beta` as relating to the width of the ring, and +# :math:`\gamma` as the radius of the ring. # -# Here we can think of :math:`\alpha` as relating to the amplitude of the monitor -# function, :math:`\beta` as relating to the width of the ring, and :math:`\gamma` as -# the radius of the ring. - - -def ring_monitor(mesh): - alpha = Constant(20.0) - beta = Constant(200.0) - gamma = Constant(0.15) - x, y = SpatialCoordinate(mesh) - r = (x - 0.5) ** 2 + (y - 0.5) ** 2 - return Constant(1.0) + alpha / cosh(beta * (r - gamma)) ** 2 +# For convenience, Movement provides a builder class for ring monitors, +# :class:`~movement.monitor.RingMonitorBuilder`, amongst other commonly used +# monitor functions. :: +mb = RingMonitorBuilder(centre=(0.5, 0.5), radius=0.4, amplitude=20.0, width=200.0) +ring_monitor = mb.get_monitor() # With an initial mesh and a monitor function, we are able to construct a # :class:`~movement.monge_ampere.MongeAmpereMover` instance and adapt the mesh. By @@ -102,60 +97,82 @@ def ring_monitor(mesh): # # .. 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. +# 0 Volume ratio 12.91 Variation (σ/μ) 1.21e+00 Residual 1.16e+00 +# 1 Volume ratio 6.77 Variation (σ/μ) 6.15e-01 Residual 4.85e-01 +# 2 Volume ratio 5.84 Variation (σ/μ) 5.69e-01 Residual 4.26e-01 +# 3 Volume ratio 6.34 Variation (σ/μ) 5.14e-01 Residual 3.74e-01 +# 4 Volume ratio 6.38 Variation (σ/μ) 4.88e-01 Residual 3.22e-01 +# 5 Volume ratio 12.05 Variation (σ/μ) 4.55e-01 Residual 2.57e-01 +# 6 Volume ratio 11.69 Variation (σ/μ) 4.26e-01 Residual 2.30e-01 +# 7 Volume ratio 11.98 Variation (σ/μ) 4.23e-01 Residual 1.98e-01 +# 8 Volume ratio 11.81 Variation (σ/μ) 4.17e-01 Residual 1.95e-01 +# 9 Volume ratio 12.00 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 10 Volume ratio 12.07 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 11 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 12 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 13 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 14 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 15 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 16 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 17 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 18 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 19 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 20 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 21 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 22 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 23 Volume ratio 12.09 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 24 Volume ratio 12.10 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 25 Volume ratio 12.18 Variation (σ/μ) 4.16e-01 Residual 1.94e-01 +# 26 Volume ratio 12.38 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 27 Volume ratio 12.43 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 28 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 29 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 30 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 31 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 32 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 33 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 34 Volume ratio 12.44 Variation (σ/μ) 4.16e-01 Residual 1.93e-01 +# 35 Volume ratio 13.34 Variation (σ/μ) 4.18e-01 Residual 1.94e-01 +# 36 Volume ratio 15.21 Variation (σ/μ) 4.38e-01 Residual 2.09e-01 +# 37 Volume ratio 17.20 Variation (σ/μ) 4.76e-01 Residual 2.22e-01 +# 38 Volume ratio 11.89 Variation (σ/μ) 4.30e-01 Residual 1.45e-01 +# 39 Volume ratio 12.41 Variation (σ/μ) 4.13e-01 Residual 1.01e-01 +# 40 Volume ratio 8.76 Variation (σ/μ) 3.79e-01 Residual 7.36e-02 +# 41 Volume ratio 10.00 Variation (σ/μ) 3.85e-01 Residual 6.36e-02 +# 42 Volume ratio 8.57 Variation (σ/μ) 3.65e-01 Residual 5.17e-02 +# 43 Volume ratio 9.48 Variation (σ/μ) 3.72e-01 Residual 4.17e-02 +# 44 Volume ratio 8.70 Variation (σ/μ) 3.60e-01 Residual 3.51e-02 +# 45 Volume ratio 9.31 Variation (σ/μ) 3.64e-01 Residual 2.69e-02 +# 46 Volume ratio 8.82 Variation (σ/μ) 3.57e-01 Residual 2.23e-02 +# 47 Volume ratio 9.20 Variation (σ/μ) 3.58e-01 Residual 1.44e-02 +# 48 Volume ratio 8.95 Variation (σ/μ) 3.53e-01 Residual 1.03e-02 +# 49 Volume ratio 9.08 Variation (σ/μ) 3.52e-01 Residual 4.91e-03 +# 50 Volume ratio 9.05 Variation (σ/μ) 3.50e-01 Residual 3.98e-03 +# 51 Volume ratio 9.13 Variation (σ/μ) 3.51e-01 Residual 2.81e-03 +# 52 Volume ratio 9.11 Variation (σ/μ) 3.50e-01 Residual 2.28e-03 +# 53 Volume ratio 9.17 Variation (σ/μ) 3.50e-01 Residual 1.51e-03 +# 54 Volume ratio 9.15 Variation (σ/μ) 3.50e-01 Residual 1.11e-03 +# 55 Volume ratio 9.19 Variation (σ/μ) 3.50e-01 Residual 5.28e-04 +# 56 Volume ratio 9.18 Variation (σ/μ) 3.49e-01 Residual 2.89e-04 +# 57 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 1.04e-04 +# 58 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 7.69e-05 +# 59 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 5.31e-05 +# 60 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 3.80e-05 +# 61 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 2.66e-05 +# 62 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 1.87e-05 +# 63 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 1.36e-05 +# 64 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 9.55e-06 +# 65 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 7.22e-06 +# 66 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 5.04e-06 +# 67 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 3.88e-06 +# 68 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 2.50e-06 +# 69 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 1.87e-06 +# 70 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 9.04e-07 +# 71 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 6.32e-07 +# 72 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 2.67e-07 +# 73 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 1.55e-07 +# 74 Volume ratio 9.19 Variation (σ/μ) 3.49e-01 Residual 8.22e-08 +# Solver converged in 74 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. diff --git a/movement/__init__.py b/movement/__init__.py index c00d0e1..c5770f1 100644 --- a/movement/__init__.py +++ b/movement/__init__.py @@ -1,3 +1,4 @@ +from movement.monitor import * # noqa from movement.tangling import * # noqa from movement.laplacian_smoothing import * # noqa diff --git a/movement/monitor.py b/movement/monitor.py new file mode 100644 index 0000000..6ee5e50 --- /dev/null +++ b/movement/monitor.py @@ -0,0 +1,365 @@ +import abc + +import ufl +from animate.recovery import recover_gradient_l2, recover_hessian_clement +from animate.utility import norm +from firedrake import SpatialCoordinate +from firedrake.constant import Constant +from firedrake.function import Function +from firedrake.functionspace import ( + FunctionSpace, + TensorFunctionSpace, + VectorFunctionSpace, +) + +__all__ = [ + "ConstantMonitorBuilder", + "BallMonitorBuilder", + "RingMonitorBuilder", + "GradientMonitorBuilder", + "HessianMonitorBuilder", + "GradientHessianMonitorBuilder", +] + + +class MonitorBuilder(metaclass=abc.ABCMeta): + """ + Abstract base class for monitor function factories. + """ + + def __init__(self, dim): + """ + :arg dim: mesh dimension + :type dim: :class:`int` + """ + self.dim = dim + + @abc.abstractmethod + def _monitor(self, mesh): + """ + Abstract method to create a monitor function. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: monitor function evaluated on given mesh + :rtype: :class:`firedrake.function.Function` + """ + pass # pragma: no cover + + def get_monitor(self): + """ + Returns a callable monitor function whose only argument is the mesh. + + :return: monitor function + :rtype: callable monitor function with a single argument + """ + + def monitor(mesh): + m = self._monitor(mesh) + if not isinstance(m, (Constant, Function)): + m = Function(FunctionSpace(mesh, "CG", 1)).interpolate(m) + return m + + return monitor + + def __call__(self): + """ + Alias for :meth:`get_monitor`. + + :return: monitor function + :rtype: callable monitor function with a single argument + """ + return self.get_monitor() + + +class ConstantMonitorBuilder(MonitorBuilder): + """ + Builder class for constant monitor functions. + """ + + def _monitor(self, mesh): + """ + Creates a constant monitor function. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: constant monitor function + :rtype: :class:`firedrake.constant.Constant` + """ + return Constant(1.0) + + +class BallMonitorBuilder(MonitorBuilder): + r""" + Builder class for monitor functions focused around ball shapes: + + .. math:: + m(\mathbf{x}) = 1 + \frac{\alpha} + {\cosh^2\left(\beta\left((\mathbf{x}-\mathbf{c})\cdot(\mathbf{x}-\mathbf{c}) + -\gamma^2\right)\right)}, + + where :math:`\mathbf{c}` is the centre point, :math:`\alpha` is the amplitude of the + monitor function, :math:`\beta` is the width of the transition region, and + :math:`\gamma` is the radius of the ball. + """ + + def __init__(self, dim, centre, radius, amplitude, width): + r""" + :arg dim: mesh dimension + :type dim: :class:`int` + :arg centre: the centre of the ball + :type centre: :class:`tuple` of :class:`float`\s + :arg radius: the radius of the ball + :type radius: :class:`float` + :arg amplitude: the amplitude of the monitor function + :type amplitude: :class:`float` + :arg width: the width of the transition region + :type width: :class:`float` + """ + assert len(centre) == dim + assert radius > 0 + assert amplitude > 0 + assert width > 0 + super().__init__(dim) + self.centre = ufl.as_vector([Constant(c) for c in centre]) + self.radius = Constant(radius) + self.amplitude = Constant(amplitude) + self.width = Constant(width) + + def _monitor(self, mesh): + """ + Creates a monitor function focused around a ball shape. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: expression of the ball-shaped monitor function + :rtype: :class:`ufl.core.expr.Expr` + """ + diff = SpatialCoordinate(mesh) - self.centre + dist = ufl.dot(diff, diff) + return ( + Constant(1.0) + + self.amplitude / ufl.cosh(self.width * (dist - self.radius**2)) ** 2 + ) + + +class RingMonitorBuilder(BallMonitorBuilder): + r""" + Builder class for monitor functions focused around 2D ring shapes: + + .. math:: + m(x,y) = 1 + \frac{\alpha} + {\cosh^2\left(\beta\left((x-c_0)^2+(y-c_1)^2\right) + -\gamma^2\right)}, + + where :math:`(c_0,c_1)` is the centre point, :math:`\alpha` is the amplitude of the + monitor function, :math:`\beta` is the width of the transition region, and + :math:`\gamma` is the radius of the ball. + """ + + def __init__(self, centre, radius, amplitude, width): + r""" + :arg centre: the centre of the ball + :type centre: :class:`tuple` of :class:`float`\s + :arg radius: the radius of the ball + :type radius: :class:`float` + :arg amplitude: the amplitude of the monitor function + :type amplitude: :class:`float` + :arg width: the width of the transition region + :type width: :class:`float` + """ + super().__init__(2, centre, radius, amplitude, width) + + +class SolutionBasedMonitorBuilder(MonitorBuilder, metaclass=abc.ABCMeta): + """ + Abstract base class for monitor factories based on solution data. + """ + + @abc.abstractmethod + def __init__(self, dim, solution): + """ + :arg dim: mesh dimension + :type dim: :class:`int` + :arg solution: solution to base the monitor on + :type solution: :class:`firedrake.function.Function` + """ + super().__init__(dim) + assert isinstance(solution, Function) + self.solution = solution + + def projection(self, mesh): + """ + Project the solution field onto the given mesh. + + :arg mesh: the mesh on which the solution is to be projected + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: the projected solution field + :rtype: :class:`firedrake.function.Function` + """ + fs = FunctionSpace(mesh, self.solution.ufl_element()) + return Function(fs).project(self.solution) + + +# TODO: Support computing gradient with Clement interpolant +class GradientMonitorBuilder(SolutionBasedMonitorBuilder): + r""" + Builder class for monitor functions based on gradients of solutions: + + .. math:: + m(\mathbf{x}) = 1 + \alpha\frac{\nabla u\cdot\nabla u} + {\max_{\mathbf{x}\in\Omega}\nabla u\cdot\nabla u}, + + where :math:`\alpha` is a scale factor and :math:`u` is the solution field of + interest. + """ + + def __init__(self, dim, solution, scale_factor): + """ + :arg dim: mesh dimension + :type dim: :class:`int` + :arg solution: solution to recover the gradient of + :type solution: :class:`firedrake.function.Function` + :arg scale_factor: scale factor + :type scale_factor: :class:`float` + """ + super().__init__(dim, solution) + assert scale_factor > 0 + self.gradient_scale_factor = Constant(scale_factor) + + def recover_gradient(self, target_space): + r""" + Recover the gradient of the solution field projected onto the current mesh. + + :arg target_space: space to recover gradient in + :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace` + :return: the recovered gradient in vector :math:`\mathbb{P}1` space + :rtype: :class:`firedrake.function.Function` + """ + mesh = target_space.mesh() + return recover_gradient_l2(self.projection(mesh), target_space=target_space) + + def _monitor(self, mesh): + """ + Monitor function based on recovered gradient. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: expression of the gradient-based monitor function on given mesh + :rtype: :class:`ufl.core.expr.Expr` + """ + g = self.recover_gradient(VectorFunctionSpace(mesh, "CG", 1)) + gg = Function(FunctionSpace(mesh, "CG", 1)).interpolate(ufl.dot(g, g)) + return Constant(1.0) + self.gradient_scale_factor * ( + gg / norm(gg, norm_type="linf") + ) + + +# TODO: Support computing Hessian with double L2 projection +class HessianMonitorBuilder(SolutionBasedMonitorBuilder): + r""" + Builder class for monitor functions based on Hessians of solutions. + + .. math:: + m(\mathbf{x}) = 1 + \beta\frac{\mathbf{H}(u):\mathbf{H}(u)} + {\max_{\mathbf{x}\in\Omega}\mathbf{H}(u):\mathbf{H}(u)}, + + where :math:`\beta` is a scale factor, :math:`u` is the solution field of interest, + and :math:`\mathbf{H}(u)` is the Hessian + """ + + def __init__(self, dim, solution, scale_factor): + """ + :arg dim: mesh dimension + :type dim: :class:`int` + :arg solution: solution to recover the Hessian of + :type solution: :class:`firedrake.function.Function` + :arg scale_factor: scale factor + :type scale_factor: :class:`float` + """ + super().__init__(dim, solution) + assert scale_factor > 0 + self.hessian_scale_factor = Constant(scale_factor) + + def recover_hessian(self, target_space): + r""" + Recover the Hessian of the solution field. + + :arg target_space: space to recover Hessian in + :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace` + :return: the recovered Hessian in tensor :math:`\mathbb{P}1` space + :rtype: :class:`firedrake.function.Function` + """ + return recover_hessian_clement(self.projection(target_space.mesh()))[1] + + def _monitor(self, mesh): + """ + Monitor function based on recovered Hessian. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: expression of the Hessian-based monitor function on given mesh + :rtype: :class:`ufl.core.expr.Expr` + """ + H = self.recover_hessian(TensorFunctionSpace(mesh, "CG", 1)) + HH = Function(FunctionSpace(mesh, "CG", 1)).interpolate(ufl.inner(H, H)) + return Constant(1.0) + self.hessian_scale_factor * ( + HH / norm(HH, norm_type="linf") + ) + + +class GradientHessianMonitorBuilder(GradientMonitorBuilder, HessianMonitorBuilder): + r""" + Builder class for monitor functions based on both gradients and Hessians of + solutions. + + .. math:: + m(\mathbf{x}) = 1 + \alpha\frac{\nabla u\cdot\nabla u} + {\max_{\mathbf{x}\in\Omega}\nabla u\cdot\nabla u} + + \beta\frac{\mathbf{H}(u):\mathbf{H}(u)} + {\max_{\mathbf{x}\in\Omega}\mathbf{H}(u):\mathbf{H}(u)}, + + where :math:`\alpha` is a scale factor for the gradient part, :math:`\beta` is a + scale factor for the Hessian part, :math:`u` is the solution field of interest, and + :math:`\mathbf{H}(u)` is the Hessian + """ + + def __init__(self, dim, solution, gradient_scale_factor, hessian_scale_factor): + """ + :arg dim: mesh dimension + :type dim: :class:`int` + :arg gradient_scale_factor: scale factor for the gradient part + :type gradient_scale_factor: :class:`float` + :arg hessian_scale_factor: scale factor for the Hessian part + :type hessian_scale_factor: :class:`float` + :arg solution: solution to recover the gradient and Hessian of + :type solution: :class:`firedrake.function.Function` + """ + SolutionBasedMonitorBuilder.__init__(self, dim, solution) + self.gradient_scale_factor = Constant(gradient_scale_factor) + self.hessian_scale_factor = Constant(hessian_scale_factor) + + def _monitor(self, mesh): + """ + Monitor function based on recovered gradient and Hessian. + + :arg mesh: the mesh on which the monitor function is to be defined + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :return: expression of the gradient- and Hessian-based monitor function + on given mesh + :rtype: :class:`ufl.core.expr.Expr` + """ + # Recover gradient + g = self.recover_gradient(VectorFunctionSpace(mesh, "CG", 1)) + gg = Function(FunctionSpace(mesh, "CG", 1)).interpolate(ufl.dot(g, g)) + + # Recover Hessian + H = self.recover_hessian(TensorFunctionSpace(mesh, "CG", 1)) + HH = Function(FunctionSpace(mesh, "CG", 1)).interpolate(ufl.inner(H, H)) + + # Combine both gradient and Hessian parts + return ( + Constant(1.0) + + self.gradient_scale_factor * (gg / norm(gg, norm_type="linf")) + + self.hessian_scale_factor * (HH / norm(HH, norm_type="linf")) + ) diff --git a/test/test_monitor.py b/test/test_monitor.py new file mode 100644 index 0000000..b3c343e --- /dev/null +++ b/test/test_monitor.py @@ -0,0 +1,152 @@ +""" +Unit tests for monitor module. +""" + +import unittest + +import numpy as np +from firedrake import SpatialCoordinate, UnitSquareMesh +from firedrake.function import Function +from firedrake.functionspace import FunctionSpace +from firedrake.norms import errornorm + +from movement.monitor import ( + BallMonitorBuilder, + ConstantMonitorBuilder, + GradientHessianMonitorBuilder, + GradientMonitorBuilder, + HessianMonitorBuilder, +) + + +class BaseClasses: + """ + Base classes for monitor factories. + """ + + class TestMonitorBuilder(unittest.TestCase): + """ + Base class for monitor factory unit tests. + """ + + def setUp(self): + self.mesh = UnitSquareMesh(5, 5) + self.P1 = FunctionSpace(self.mesh, "CG", 1) + self.solution = Function(self.P1) + + +class TestConstant(BaseClasses.TestMonitorBuilder): + """ + Unit tests for :class:`~.ConstantMonitorBuilder`. + """ + + def test_value_get_monitor(self): + mb = ConstantMonitorBuilder(self.mesh.topological_dimension()) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 1)) + + def test_value_call(self): + mb = ConstantMonitorBuilder(self.mesh.topological_dimension()) + self.assertTrue(np.allclose(mb()(self.mesh).dat.data, 1)) + + +class TestBall(BaseClasses.TestMonitorBuilder): + """ + Unit tests for :class:`~.BallMonitorBuilder`. + """ + + def test_tiny_amplitude(self): + mb = BallMonitorBuilder( + dim=2, centre=(0, 0), radius=0.1, amplitude=1e-8, width=0.1 + ) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 1)) + + +class TestGradient(BaseClasses.TestMonitorBuilder): + """ + Unit tests for :class:`~.GradientMonitorBuilder`. + """ + + def test_tiny_scale_factor(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x**2) + mb = GradientMonitorBuilder(dim=2, solution=self.solution, scale_factor=1e-8) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 1)) + + def test_linear(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x) + mb = GradientMonitorBuilder(dim=2, solution=self.solution, scale_factor=1) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 2)) + + +class TestHessian(BaseClasses.TestMonitorBuilder): + """ + Unit tests for :class:`~.HessianMonitorBuilder`. + """ + + def test_tiny_scale_factor(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x**3) + mb = HessianMonitorBuilder(dim=2, solution=self.solution, scale_factor=1e-8) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 1)) + + def test_quadratic(self): + x, y = SpatialCoordinate(self.mesh) + P2 = FunctionSpace(self.mesh, "CG", 2) + solution = Function(P2) + solution.interpolate(0.5 * x**2) + mb = HessianMonitorBuilder(dim=2, solution=solution, scale_factor=1) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 2)) + + +class TestGradientHessian(BaseClasses.TestMonitorBuilder): + """ + Unit tests for :class:`~.GradientHessianMonitorBuilder`. + """ + + def test_tiny_scale_factors(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x**3) + mb = GradientHessianMonitorBuilder( + dim=2, + solution=self.solution, + gradient_scale_factor=1e-8, + hessian_scale_factor=1e-8, + ) + self.assertTrue(np.allclose(mb.get_monitor()(self.mesh).dat.data, 1)) + + def test_tiny_hessian_scale_factor(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x**3) + mb1 = GradientHessianMonitorBuilder( + dim=2, + solution=self.solution, + gradient_scale_factor=1, + hessian_scale_factor=1e-8, + ) + mb2 = GradientMonitorBuilder( + dim=2, + solution=self.solution, + scale_factor=1, + ) + self.assertAlmostEqual( + errornorm(mb1.get_monitor()(self.mesh), mb2.get_monitor()(self.mesh)), 0 + ) + + def test_tiny_gradient_scale_factor(self): + x, y = SpatialCoordinate(self.mesh) + self.solution.interpolate(x**3) + mb1 = GradientHessianMonitorBuilder( + dim=2, + solution=self.solution, + gradient_scale_factor=1e-8, + hessian_scale_factor=1, + ) + mb2 = HessianMonitorBuilder( + dim=2, + solution=self.solution, + scale_factor=1, + ) + self.assertAlmostEqual( + errornorm(mb1.get_monitor()(self.mesh), mb2.get_monitor()(self.mesh)), 0 + )