Mixed fs limiter

gusto/limiters.py

@@ -12,7 +12,7 @@
 
 import numpy as np
 
-__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter"]
+__all__ = ["DG1Limiter", "ThetaLimiter", "NoLimiter", "MixedFSLimiter"]
 
 
 class DG1Limiter(object):
@@ -178,3 +178,35 @@ def apply(self, field):
                 applied, if this was not a blank limiter.
         """
         pass
+
+
+class MixedFSLimiter(object):
+    """An object to hold a dictionary that defines limiters for transported prognostic
+    variables. Different limiters may be applied to different fields and not every transported variable needs a defined limiter.
+    """
+
+    def __init__(self, equation, sublimiters):
+        """
+        Args:
+            equation (:class: 'PrognosticEquationSet'): the prognostic equation(s)
+            sublimiters (dict): A dictionary holding limiters defined for individual prognostic variables
+        Raises:
+            ValueError: If a limiter is defined for a field that is not in the prognostic variable set
+        """
+        self.sublimiters = sublimiters
+
+        for field, sublimiter in sublimiters.items():
+            # Check that the field is in the prognostic variable set:
+            if field not in equation.field_names:
+                raise ValueError(f"The limiter defined for {field} is for a field that does not exist in the equation set")
+            else:
+                self.sublimiters[field].idx = equation.field_names.index(field)
+
+    def apply(self, fields):
+        """
+        Apply the individual limiters to specific prognostic variables
+        """
+
+        for _, sublimiter in self.sublimiters.items():
+            field = fields.subfunctions[sublimiter.idx]
+            sublimiter.apply(field)

integration-tests/transport/test_limiters.py

@@ -50,17 +50,32 @@ def setup_limiters(dirname, space):
         VDG1 = domain.spaces('DG1_equispaced')
     elif space == 'Vtheta_degree_1':
         V = domain.spaces('theta')
+    elif space == 'mixed_FS':
+        V = domain.spaces('HDiv')
+        VA = domain.spaces('DG')
+        VB = domain.spaces('DG1_equispaced')
     else:
         raise NotImplementedError
 
     Vpsi = domain.spaces('H1')
 
     # Equation
-    eqn = AdvectionEquation(domain, V, 'tracer')
+    if space == 'mixed_FS':
+        tracerA = ActiveTracer(name='tracerA', space='DG',
+                               variable_type=TracerVariableType.mixing_ratio,
+                               transport_eqn=TransportEquationType.advective)
+        tracerB = ActiveTracer(name='tracerB', space='DG1_equispaced',
+                               variable_type=TracerVariableType.mixing_ratio,
+                               transport_eqn=TransportEquationType.advective)
+        tracers = [tracerA, tracerB]
+        eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V)
+        output = OutputParameters(dirname=dirname+'/limiters', dumpfreq=1,
+                                  dumplist=['u', 'tracerA', 'tracerB', 'true_tracerA', 'true_tracerB'])
+    else:
+        eqn = AdvectionEquation(domain, V, 'tracer')
+        output = OutputParameters(dirname=dirname+'/limiters',
+                                  dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer'])
 
-    # I/O
-    output = OutputParameters(dirname=dirname+'/limiters',
-                              dumpfreq=1, dumplist=['u', 'tracer', 'true_tracer'])
     io = IO(domain, output)
 
     # ------------------------------------------------------------------------ #
@@ -84,10 +99,19 @@ def setup_limiters(dirname, space):
     elif space == 'Vtheta_degree_1':
         opts = EmbeddedDGOptions()
         transport_schemes = SSPRK3(domain, options=opts, limiter=ThetaLimiter(V))
+
+    elif space == 'mixed_FS':
+        sublimiters = {'tracerA': DG1Limiter(domain.spaces('DG')),
+                       'tracerB': VertexBasedLimiter(domain.spaces('DG1_equispaced'))}
+        MixedLimiter = MixedFSLimiter(eqn, sublimiters)
+        transport_schemes = SSPRK3(domain, limiter=MixedLimiter)
     else:
         raise NotImplementedError
 
-    transport_method = DGUpwind(eqn, "tracer")
+    if space == 'mixed_FS':
+        transport_method = [DGUpwind(eqn, 'tracerA'), DGUpwind(eqn, 'tracerB')]
+    else:
+        transport_method = DGUpwind(eqn, "tracer")
 
     # Build time stepper
     stepper = PrescribedTransport(eqn, transport_schemes, io, transport_method)
@@ -96,8 +120,14 @@ def setup_limiters(dirname, space):
     # Initial condition
     # ------------------------------------------------------------------------ #
 
-    tracer0 = stepper.fields('tracer', V)
-    true_field = stepper.fields('true_tracer', space=V)
+    if space == 'mixed_FS':
+        tracerA_0 = stepper.fields('tracerA', space=VA)
+        tracerB_0 = stepper.fields('tracerB', space=VB)
+        true_fieldA = stepper.fields('true_tracerA', space=VA)
+        true_fieldB = stepper.fields('true_tracerB', space=VB)
+    else:
+        tracer0 = stepper.fields('tracer', V)
+        true_field = stepper.fields('true_tracer', space=V)
 
     x, z = SpatialCoordinate(mesh)
 
@@ -146,9 +176,17 @@ def setup_limiters(dirname, space):
                              0.0)
 
         if i == 0:
-            tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+            if space == 'mixed_FS':
+                tracerA_0.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+                tracerB_0.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+            else:
+                tracer0.interpolate(Constant(tracer_min) + expr_1 + expr_2)
         elif i == 1:
-            true_field.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+            if space == 'mixed_FS':
+                true_fieldA.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+                true_fieldB.interpolate(Constant(tracer_min) + expr_1 + expr_2)
+            else:
+                true_field.interpolate(Constant(tracer_min) + expr_1 + expr_2)
         else:
             raise ValueError
 
@@ -181,29 +219,67 @@ def setup_limiters(dirname, space):
     gradperp = lambda v: as_vector([-v.dx(1), v.dx(0)])
     u.project(gradperp(psi))
 
-    return stepper, tmax, true_field
+    if space == 'mixed_FS':
+        return stepper, tmax, true_fieldA, true_fieldB
+    else:
+        return stepper, tmax, true_field
 
 
-@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1',
-                                   'DG0', 'DG1', 'DG1_equispaced'])
+@pytest.mark.parametrize('space', ['Vtheta_degree_0', 'Vtheta_degree_1', 'DG0',
+                                   'DG1', 'DG1_equispaced', 'mixed_FS'])
 def test_limiters(tmpdir, space):
 
     # Setup and run
     dirname = str(tmpdir)
-    stepper, tmax, true_field = setup_limiters(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'
+    if space == 'mixed_FS':
+        stepper, tmax, true_fieldA, true_fieldB = setup_limiters(dirname, space)
+    else:
+        stepper, tmax, true_field = setup_limiters(dirname, space)
+
+    stepper.run(t=0, tmax=tmax)
 
     tol = 1e-9
 
-    # Check for no new overshoots
-    assert np.max(final_field.dat.data) <= np.max(true_field.dat.data) + tol, \
-        'Application of limiter has not prevented overshoots'
+    if space == 'mixed_FS':
+        final_fieldA = stepper.fields('tracerA')
+        final_fieldB = stepper.fields('tracerB')
+
+        # Check tracer is roughly in the correct place
+        assert norm(true_fieldA - final_fieldA) / norm(true_fieldA) < 0.05, \
+            'Something is wrong with the DG space tracer using a mixed limiter'
+
+        # Check tracer is roughly in the correct place
+        assert norm(true_fieldB - final_fieldB) / norm(true_fieldB) < 0.05, \
+            'Something is wrong with the DG1 equispaced tracer using a mixed limiter'
+
+        # Check for no new overshoots in A
+        assert np.max(final_fieldA.dat.data) <= np.max(true_fieldA.dat.data) + tol, \
+            'Application of the DG space limiter in the mixed limiter has not prevented overshoots'
+
+        # Check for no new undershoots in A
+        assert np.min(final_fieldA.dat.data) >= np.min(true_fieldA.dat.data) - tol, \
+            'Application of the DG space limiter in the mixed limiter has not prevented undershoots'
+
+        # Check for no new overshoots in B
+        assert np.max(final_fieldB.dat.data) <= np.max(true_fieldB.dat.data) + tol, \
+            'Application of the DG1 equispaced limiter in the mixed limiter has not prevented overshoots'
+
+        # Check for no new undershoots in B
+        assert np.min(final_fieldB.dat.data) >= np.min(true_fieldB.dat.data) - tol, \
+            'Application of the DG1 equispaced limiter in the mixed limiter has not prevented undershoots'
+
+    else:
+        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.max(final_field.dat.data) <= np.max(true_field.dat.data) + tol, \
+            'Application of limiter has not prevented overshoots'
 
-    # Check for no new undershoots
-    assert np.min(final_field.dat.data) >= np.min(true_field.dat.data) - tol, \
-        'Application of limiter has not prevented undershoots'
+        # Check for no new undershoots
+        assert np.min(final_field.dat.data) >= np.min(true_field.dat.data) - tol, \
+            'Application of limiter has not prevented undershoots'