Limiter to clip fields to zero

examples/compressible/unsaturated_bubble.py

@@ -67,6 +67,7 @@
 VCG1 = FunctionSpace(mesh, "CG", 1)
 Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element())
 Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1)
+Vt = domain.spaces("theta")
 
 u_opts = RecoveryOptions(embedding_space=Vu_DG1,
                          recovered_space=Vu_CG1,
@@ -92,10 +93,12 @@
 # Physics schemes
 # NB: can't yet use wrapper or limiter options with physics
 rainfall_method = DGUpwind(eqns, 'rain', outflow=True)
+zero_limiter = MixedFSLimiter(eqns, {'water_vapour': ZeroLimiter(Vt),
+                                     'cloud_water': ZeroLimiter(Vt)})
 physics_schemes = [(Fallout(eqns, 'rain', domain, rainfall_method), SSPRK3(domain)),
                    (Coalescence(eqns), ForwardEuler(domain)),
                    (EvaporationOfRain(eqns), ForwardEuler(domain)),
-                   (SaturationAdjustment(eqns), ForwardEuler(domain))]
+                   (SaturationAdjustment(eqns), ForwardEuler(domain, limiter=zero_limiter))]
 
 # Time stepper
 stepper = SemiImplicitQuasiNewton(eqns, io, transported_fields,
@@ -116,7 +119,6 @@
 
 # spaces
 Vu = domain.spaces("HDiv")
-Vt = domain.spaces("theta")
 Vr = domain.spaces("DG")
 x, z = SpatialCoordinate(mesh)
 quadrature_degree = (4, 4)

gusto/kernels.py

@@ -75,3 +75,40 @@ def apply(self, field_hat, field_DG1, field_old):
                   "field_DG1": (field_DG1, READ),
                   "field_old": (field_old, READ)},
                  is_loopy_kernel=True)
+
+
+class ClipZero():
+    """Clips any negative field values to be zero."""
+
+    def __init__(self, V):
+        """
+        Args:
+            V (:class:`FunctionSpace`): The space of the field to be clipped.
+        """
+        shapes = {'nDOFs': V.finat_element.space_dimension()}
+        domain = "{{[i]: 0 <= i < {nDOFs}}}".format(**shapes)
+
+        instrs = ("""
+                  for i
+                      if field_in[i] < 0.0
+                          field[i] = 0.0
+                      else
+                          field[i] = field_in[i]
+                      end
+                  end
+                  """)
+
+        self._kernel = (domain, instrs)
+
+    def apply(self, field, field_in):
+        """
+        Performs the par loop.
+
+        Args:
+            field (:class:`Function`): The field to be written to.
+            field_in (:class:`Function`): The field to be clipped.
+        """
+        par_loop(self._kernel, dx,
+                 {"field": (field, WRITE),
+                  "field_in": (field_in, READ)},
+                 is_loopy_kernel=True)

gusto/limiters.py

@@ -8,11 +8,12 @@
 from firedrake import (BrokenElement, Function, FunctionSpace, interval,
                        FiniteElement, TensorProductElement)
 from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter
-from gusto.kernels import LimitMidpoints
+from gusto.kernels import LimitMidpoints, ClipZero
 
 import numpy as np
 
-__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter", "MixedFSLimiter"]
+__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter", "ZeroLimiter",
+           "MixedFSLimiter"]
 
 
 class DG1Limiter(object):
@@ -163,6 +164,52 @@ def apply(self, field):
         field.interpolate(self.field_hat)
 
 
+class ZeroLimiter(object):
+    """
+    A simple limiter to enforce non-negativity of a field pointwise.
+
+    Negative values are simply clipped to be zero. There is also the option to
+    project the field to another function space to enforce non-negativity there.
+    """
+
+    def __init__(self, space, clipping_space=None):
+        """
+        Args:
+            space (:class:`FunctionSpace`): the space of the incoming field to
+                clip.
+            clipping_space (:class:`FunctionSpace`, optional): the space in
+                which to clip the field. If not specified, the space of the
+                input field is used.
+        """
+
+        self.space = space
+        if clipping_space is not None:
+            self.clipping_space = clipping_space
+            self.map_to_clip = True
+            self.field_to_clip = Function(self.clipping_space)
+        else:
+            self.clipping_space = space
+            self.map_to_clip = False
+
+        self._kernel = ClipZero(self.clipping_space)
+
+    def apply(self, field):
+        """
+        The application of the limiter to the field.
+
+        Args:
+            field (:class:`Function`): the field to apply the limiter to.
+         """
+
+        # Obtain field in clipping space
+        if self.map_to_clip:
+            self.field_to_clip.interpolate(field)
+            self._kernel.apply(self.field_to_clip, self.field_to_clip)
+            field.interpolate(self.field_to_clip)
+        else:
+            self._kernel.apply(field, field)
+
+
 class NoLimiter(object):
     """A blank limiter that does nothing."""
 

integration-tests/transport/test_zero_limiter.py

@@ -0,0 +1,174 @@
+"""
+This tests the ZeroLimiter, which enforces non-negativity.
+A sharp bubble of warm air is generated in a vertical slice and then transported
+by a prescribed transport scheme. If the limiter is working, the transport
+should have produced no negative values.
+"""
+
+from gusto import *
+from firedrake import (as_vector, PeriodicIntervalMesh, pi, SpatialCoordinate,
+                       ExtrudedMesh, FunctionSpace, Function, norm,
+                       conditional, sqrt)
+import numpy as np
+import pytest
+
+
+def setup_zero_limiter(dirname, clipping_space):
+
+    # ------------------------------------------------------------------------ #
+    # Parameters for test case
+    # ------------------------------------------------------------------------ #
+
+    Ld = 1.
+    tmax = 0.2
+    dt = tmax / 40
+    rotations = 0.25
+
+    # ------------------------------------------------------------------------ #
+    # Build model objects
+    # ------------------------------------------------------------------------ #
+
+    # Domain
+    m = PeriodicIntervalMesh(20, Ld)
+    mesh = ExtrudedMesh(m, layers=20, layer_height=(Ld/20))
+    degree = 1
+
+    domain = Domain(mesh, dt, family="CG", degree=degree)
+
+    DG1 = FunctionSpace(mesh, 'DG', 1)
+    DG1_equispaced = domain.spaces('DG1_equispaced')
+
+    Vpsi = domain.spaces('H1')
+
+    eqn = AdvectionEquation(domain, DG1, 'tracer')
+    output = OutputParameters(dirname=dirname+'/limiters',
+                              dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer'])
+
+    io = IO(domain, output)
+
+    # ------------------------------------------------------------------------ #
+    # Set up transport scheme
+    # ------------------------------------------------------------------------ #
+
+    if clipping_space is None:
+        limiter = ZeroLimiter(DG1)
+    elif clipping_space == 'equispaced':
+        limiter = ZeroLimiter(DG1, clipping_space=DG1_equispaced)
+
+    transport_schemes = SSPRK3(domain, limiter=limiter)
+    transport_method = DGUpwind(eqn, "tracer")
+
+    # Build time stepper
+    stepper = PrescribedTransport(eqn, transport_schemes, io, transport_method)
+
+    # ------------------------------------------------------------------------ #
+    # Initial condition
+    # ------------------------------------------------------------------------ #
+
+    tracer0 = stepper.fields('tracer', DG1)
+    true_field = stepper.fields('true_tracer', space=DG1)
+
+    x, z = SpatialCoordinate(mesh)
+
+    tracer_min = 12.6
+    dtracer = 3.2
+
+    # First time do initial conditions, second time do final conditions
+    for i in range(2):
+
+        if i == 0:
+            x1_lower = 2 * Ld / 5
+            x1_upper = 3 * Ld / 5
+            z1_lower = 6 * Ld / 10
+            z1_upper = 8 * Ld / 10
+            x2_lower = 6 * Ld / 10
+            x2_upper = 8 * Ld / 10
+            z2_lower = 2 * Ld / 5
+            z2_upper = 3 * Ld / 5
+        elif i == 1:
+            # Rotated anti-clockwise by 90 degrees (x -> z, z -> -x)
+            x1_lower = 2 * Ld / 10
+            x1_upper = 4 * Ld / 10
+            z1_lower = 2 * Ld / 5
+            z1_upper = 3 * Ld / 5
+            x2_lower = 2 * Ld / 5
+            x2_upper = 3 * Ld / 5
+            z2_lower = 6 * Ld / 10
+            z2_upper = 8 * Ld / 10
+        else:
+            raise ValueError
+
+        expr_1 = conditional(x > x1_lower,
+                             conditional(x < x1_upper,
+                                         conditional(z > z1_lower,
+                                                     conditional(z < z1_upper, dtracer, 0.0),
+                                                     0.0),
+                                         0.0),
+                             0.0)
+
+        expr_2 = conditional(x > x2_lower,
+                             conditional(x < x2_upper,
+                                         conditional(z > z2_lower,
+                                                     conditional(z < z2_upper, dtracer, 0.0),
+                                                     0.0),
+                                         0.0),
+                             0.0)
+
+        if i == 0:
+            tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+        elif i == 1:
+            true_field.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+        else:
+            raise ValueError
+
+    # ------------------------------------------------------------------------ #
+    # Velocity profile
+    # ------------------------------------------------------------------------ #
+
+    psi = Function(Vpsi)
+    u = stepper.fields('u')
+
+    # set up solid body rotation for transport
+    # we do this slightly complicated stream function to make the velocity 0 at edges
+    # thus we avoid odd effects at boundaries
+    xc = Ld / 2
+    zc = Ld / 2
+    r = sqrt((x - xc) ** 2 + (z - zc) ** 2)
+    omega = rotations * 2 * pi / tmax
+    r_out = 9 * Ld / 20
+    r_in = 2 * Ld / 5
+    A = omega * r_in / (2 * (r_in - r_out))
+    B = - omega * r_in * r_out / (r_in - r_out)
+    C = omega * r_in ** 2 * r_out / (r_in - r_out) / 2
+    psi_expr = conditional(r < r_in,
+                           omega * r ** 2 / 2,
+                           conditional(r < r_out,
+                                       A * r ** 2 + B * r + C,
+                                       A * r_out ** 2 + B * r_out + C))
+    psi.interpolate(psi_expr)
+
+    gradperp = lambda v: as_vector([-v.dx(1), v.dx(0)])
+    u.project(gradperp(psi))
+
+    return stepper, tmax, true_field
+
+
+@pytest.mark.parametrize('space', [None, 'equispaced'])
+def test_zero_limiter(tmpdir, space):
+
+    # Setup and run
+    dirname = str(tmpdir)
+
+    stepper, tmax, true_field = setup_zero_limiter(dirname, space)
+
+    stepper.run(t=0, tmax=tmax)
+
+    final_field = stepper.fields('tracer')
+
+    # Check tracer is roughly in the correct place
+    assert norm(true_field - final_field) / norm(true_field) < 0.05, \
+        'Something appears to have gone wrong with transport of tracer using a limiter'
+
+    # Check for no new overshoots
+    assert np.min(final_field.dat.data) >= 0.0, \
+        'Application of limiter has not prevented negative values'

unit-tests/kernel_tests/test_clip_zero_kernel.py

@@ -0,0 +1,40 @@
+"""
+A test of the ClipZero kernel, which is used to enforce non-negativity.
+"""
+
+from firedrake import UnitSquareMesh, Function, FunctionSpace
+from gusto import kernels
+import numpy as np
+
+
+def test_clip_zero():
+
+    # ------------------------------------------------------------------------ #
+    # Set up meshes and spaces
+    # ------------------------------------------------------------------------ #
+
+    mesh = UnitSquareMesh(3, 3)
+
+    DG1 = FunctionSpace(mesh, "DG", 1)
+
+    field = Function(DG1)
+
+    # ------------------------------------------------------------------------ #
+    # Initial conditions
+    # ------------------------------------------------------------------------ #
+
+    field.dat.data[:] = -3.0
+
+    # ------------------------------------------------------------------------ #
+    # Apply kernel
+    # ------------------------------------------------------------------------ #
+
+    kernel = kernels.ClipZero(DG1)
+    kernel.apply(field, field)
+
+    # ------------------------------------------------------------------------ #
+    # Check values
+    # ------------------------------------------------------------------------ #
+
+    assert np.all(field.dat.data >= 0.0), \
+        'ClipZero kernel is not enforcing non-negativity'