Mixed order recovery

examples/compressible/dry_bryan_fritsch.py

@@ -9,10 +9,8 @@
 from firedrake import (IntervalMesh, ExtrudedMesh,
                        SpatialCoordinate, conditional, cos, pi, sqrt,
                        TestFunction, dx, TrialFunction, Constant, Function,
-                       LinearVariationalProblem, LinearVariationalSolver,
-                       FunctionSpace, VectorFunctionSpace)
+                       LinearVariationalProblem, LinearVariationalSolver)
 import sys
-
 # ---------------------------------------------------------------------------- #
 # Test case parameters
 # ---------------------------------------------------------------------------- #
@@ -61,19 +59,14 @@
 io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
 # Transport schemes -- set up options for using recovery wrapper
-VDG1 = domain.spaces("DG1_equispaced")
-VCG1 = FunctionSpace(mesh, "CG", 1)
-Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element())
-Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1)
-
-u_opts = RecoveryOptions(embedding_space=Vu_DG1,
-                         recovered_space=Vu_CG1,
-                         boundary_method=BoundaryMethod.taylor)
-rho_opts = RecoveryOptions(embedding_space=VDG1,
-                           recovered_space=VCG1,
-                           boundary_method=BoundaryMethod.taylor)
-theta_opts = RecoveryOptions(embedding_space=VDG1,
-                             recovered_space=VCG1)
+boundary_methods = {'DG': BoundaryMethod.taylor,
+                    'HDiv': BoundaryMethod.taylor}
+
+recovery_spaces = RecoverySpaces(domain, boundary_methods, use_vector_spaces=True)
+
+u_opts = recovery_spaces.HDiv_options
+rho_opts = recovery_spaces.DG_options
+theta_opts = recovery_spaces.theta_options
 
 transported_fields = [SSPRK3(domain, "rho", options=rho_opts),
                       SSPRK3(domain, "theta", options=theta_opts),

examples/compressible/unsaturated_bubble.py

@@ -6,11 +6,12 @@
 Limiters are applied to the transport of the water species.
 """
 from gusto import *
+from gusto.equations import thermodynamics
 from firedrake import (PeriodicIntervalMesh, ExtrudedMesh,
                        SpatialCoordinate, conditional, cos, pi, sqrt, exp,
                        TestFunction, dx, TrialFunction, Constant, Function,
                        LinearVariationalProblem, LinearVariationalSolver,
-                       FunctionSpace, VectorFunctionSpace, errornorm)
+                       errornorm)
 from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter
 import sys
 
@@ -61,21 +62,17 @@
 diagnostic_fields = [RelativeHumidity(eqns), Perturbation('theta'),
                      Perturbation('water_vapour'), Perturbation('rho'), Perturbation('RelativeHumidity')]
 io = IO(domain, output, diagnostic_fields=diagnostic_fields)
-
 # Transport schemes -- specify options for using recovery wrapper
+boundary_methods = {'DG': BoundaryMethod.taylor,
+                    'HDiv': BoundaryMethod.taylor}
+
+recovery_spaces = RecoverySpaces(domain, boundary_method=boundary_methods, use_vector_spaces=True)
+
+u_opts = recovery_spaces.HDiv_options
+rho_opts = recovery_spaces.DG_options
+theta_opts = recovery_spaces.theta_options
+
 VDG1 = domain.spaces("DG1_equispaced")
-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,
-                         boundary_method=BoundaryMethod.taylor)
-rho_opts = RecoveryOptions(embedding_space=VDG1,
-                           recovered_space=VCG1,
-                           boundary_method=BoundaryMethod.taylor)
-theta_opts = RecoveryOptions(embedding_space=VDG1, recovered_space=VCG1)
 limiter = VertexBasedLimiter(VDG1)
 
 transported_fields = [SSPRK3(domain, "u", options=u_opts),
@@ -92,6 +89,7 @@
 
 # Physics schemes
 # NB: can't yet use wrapper or limiter options with physics
+Vt = domain.spaces('theta')
 rainfall_method = DGUpwind(eqns, 'rain', outflow=True)
 zero_limiter = MixedFSLimiter(eqns, {'water_vapour': ZeroLimiter(Vt),
                                      'cloud_water': ZeroLimiter(Vt)})

gusto/recovery/__init__.py

@@ -1,3 +1,4 @@
 from gusto.recovery.averaging import *            # noqa
 from gusto.recovery.recovery import *             # noqa
-from gusto.recovery.reversible_recovery import *  # noqa
+from gusto.recovery.reversible_recovery import *  # noqa
+from gusto.recovery.recovery_options import *     # noqa

gusto/recovery/recovery_options.py

@@ -0,0 +1,117 @@
+"""Contains an object to generate the set of spaces used for the recovery wrapper """
+from firedrake import (interval, FiniteElement, TensorProductElement, FunctionSpace,
+                       VectorFunctionSpace)
+from gusto.core.function_spaces import DeRhamComplex
+from gusto.core.configuration import RecoveryOptions
+
+
+class RecoverySpaces(object):
+    """
+    Finds or builds necessary spaces to carry out recovery transport for lowest
+    and mixed order domains (0,0), (0,1) and  (1,0)
+    """
+
+    def __init__(self, domain, boundary_method=None, use_vector_spaces=False):
+        """
+        Args:
+            domain (:class:`Domain`): the model's domain object, containing the
+                mesh and the compatible function spaces.
+
+            boundary_method (:variable:'dict', optional): A dictionary containing the space
+                the boundary method is to be applied to along with specified method. Acceptable keys are "DG",
+                "HDiv" and "theta". acceptable values are (BoundaryMethod.taylor/hcurl/extruded),
+                passed as ('space', 'boundary method'). Defaults to None
+
+            use_vector_spaces (bool, optional):. Determines if we need to use DG / CG
+                space for the embedded and recovery space for the HDiv field instead of the usual
+                HDiv, HCurl spaces. Defaults to False
+        """
+        family = domain.family
+        mesh = domain.mesh
+        # Need spaces from current deRham and a higher order deRham
+        self.de_Rham = DeRhamComplex(mesh, family,
+                                     horizontal_degree=1,
+                                     vertical_degree=1,
+                                     complex_name='recovery_de_Rham')
+
+        valid_keys = ['DG', 'HDiv', 'theta']
+        if boundary_method is not None:
+            for key in boundary_method:
+                if key not in valid_keys:
+                    raise KeyError(f'Recovery spaces: boundary method key {key} not valid. Valid keys are DG, HDiv, theta')
+
+        # ----------------------------------------------------------------------
+        # Building theta options if on an extruded mesh
+        # ----------------------------------------------------------------------
+
+        # Check if extruded and if so builds theta spaces
+        if hasattr(mesh, "_base_mesh"):
+            # check if boundary method is present
+            if hasattr(boundary_method, 'theta'):
+                theta_boundary_method = boundary_method['theta']
+            else:
+                theta_boundary_method = None
+            cell = mesh._base_mesh.ufl_cell().cellname()
+            DG_hori_ele = FiniteElement('DG', cell, 1, variant='equispaced')
+            DG_vert_ele = FiniteElement('DG', interval, (domain.vertical_degree + 1), variant='equispaced')
+            CG_hori_ele = FiniteElement('CG', cell, 1)
+            CG_vert_ele = FiniteElement('CG', interval, (domain.vertical_degree + 1))
+
+            VDG_ele = TensorProductElement(DG_hori_ele, DG_vert_ele)
+            VCG_ele = TensorProductElement(CG_hori_ele, CG_vert_ele)
+            VDG_theta = FunctionSpace(mesh, VDG_ele)
+            VCG_theta = FunctionSpace(mesh, VCG_ele)
+
+            self.theta_options = RecoveryOptions(embedding_space=VDG_theta,
+                                                 recovered_space=VCG_theta,
+                                                 boundary_method=theta_boundary_method)
+        else:
+            cell = self.mesh.ufl_cell().cellname()
+
+        # ----------------------------------------------------------------------
+        # Building the DG options
+        # ----------------------------------------------------------------------
+        if hasattr(boundary_method, 'DG'):
+            DG_boundary_method = boundary_method['DG']
+        else:
+            DG_boundary_method = None
+
+        DG_embedding_space = domain.spaces.DG1_equispaced
+        # Recovered space needs builing manually to avoid uneccesary DoFs
+        CG_hori_ele_DG = FiniteElement('CG', cell, 1)
+        CG_vert_ele_DG = FiniteElement('CG', interval, 1)
+        VCG_ele_DG = TensorProductElement(CG_hori_ele_DG, CG_vert_ele_DG)
+        DG_recovered_space = FunctionSpace(mesh, VCG_ele_DG)
+
+        # DG_recovered_space = domain.spaces.H1
+        self.DG_options = RecoveryOptions(embedding_space=DG_embedding_space,
+                                          recovered_space=DG_recovered_space,
+                                          boundary_method=DG_boundary_method)
+        # ----------------------------------------------------------------------
+        # Building HDiv options
+        # ----------------------------------------------------------------------
+
+        if hasattr(boundary_method, 'HDiv'):
+            HDiv_boundary_method = boundary_method['HDiv']
+        else:
+            HDiv_boundary_method = None
+
+        if use_vector_spaces:
+            Vu_DG1 = VectorFunctionSpace(mesh, DG_embedding_space.ufl_element())
+            Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1)
+
+            HDiv_embedding_Space = Vu_DG1
+            HDiv_recovered_Space = Vu_CG1
+
+        else:
+
+            HDiv_embedding_Space = self.de_Rham.HDiv
+            HDiv_recovered_Space = self.de_Rham.HCurl
+
+        self.HDiv_options = RecoveryOptions(embedding_space=HDiv_embedding_Space,
+                                            recovered_space=HDiv_recovered_Space,
+                                            injection_method='recover',
+                                            project_high_method='project',
+                                            project_low_method='project',
+                                            broken_method='project',
+                                            boundary_method=HDiv_boundary_method)

unit-tests/recovery_tests/test_mixed_spaces.py

@@ -0,0 +1,32 @@
+"""
+Test that the Recovery_options object builds spaces of the correct degree
+"""
+from firedrake import UnitIntervalMesh, ExtrudedMesh
+from gusto.core import Domain
+from gusto.recovery.recovery_options import RecoverySpaces
+import pytest
+
+order_correct_degree_dict = {(0, 0): {'DG': (1, 1),
+                                      'HDiv': (2, 2),
+                                      'theta': (1, 1)},
+                             (0, 1): {'DG': (1, 1),
+                                      'HDiv': (2, 2),
+                                      'theta': (1, 2)},
+                             (1, 0): {'DG': (1, 1),
+                                      'HDiv': (2, 2),
+                                      'theta': (1, 1)}
+                             }
+
+
+@pytest.mark.parametrize("order", [(0, 0), (1, 0), (0, 1)])
+@pytest.mark.parametrize('space', ['DG', 'HDiv', 'theta'])
+def test_mixed_spaces(order, space):
+    mesh = UnitIntervalMesh(3)
+    emesh = ExtrudedMesh(mesh, 3, 0.33)
+    dt = 1
+    domain = Domain(emesh, dt, family='CG', horizontal_degree=order[0], vertical_degree=order[1])
+    recovery_spaces = RecoverySpaces(domain)
+
+    tesing_space = getattr(recovery_spaces, f'{space}_options')
+    degree = tesing_space.recovered_space.finat_element.degree
+    assert degree == order_correct_degree_dict[order][space]