set quadrature degree for forms in equation residual

examples/shallow_water/moist_convective_williamson2.py

@@ -82,7 +82,7 @@
 
 # define saturation function
 def sat_func(x_in):
-    h = x_in.split()[1]
+    h = x_in.subfunctions[1]
     lamda, phi, _ = lonlatr_from_xyz(x[0], x[1], x[2])
     numerator = theta_0 + sigma*((cos(phi))**2) * ((w + sigma)*(cos(phi))**2 + 2*(phi_0 - w - sigma))
     denominator = phi_0**2 + (w + sigma)**2*(sin(phi))**4 - 2*phi_0*(w + sigma)*(sin(phi))**2

examples/shallow_water/moist_thermal_williamson5.py

@@ -84,15 +84,13 @@
 
 # Saturation function
 def sat_func(x_in):
-    h = x_in.split()[1]
-    b = x_in.split()[2]
+    _, h, b = x_in.subfunctions
     return (q0/(g*h + g*tpexpr)) * exp(20*(1 - b/g))
 
 
 # Feedback proportionality is dependent on h and b
 def gamma_v(x_in):
-    h = x_in.split()[1]
-    b = x_in.split()[2]
+    _, h, b = x_in.subfunctions
     return (1 + beta2*(20*q0/(g*h + g*tpexpr) * exp(20*(1 - b/g))))**(-1)
 
 

gusto/diagnostics/diagnostics.py

@@ -6,9 +6,10 @@
                        LinearVariationalProblem, LinearVariationalSolver,
                        ds_b, ds_v, ds_t, dS_h, dS_v, ds, dS, div, avg, pi,
                        TensorFunctionSpace, SpatialCoordinate, as_vector,
-                       Projector, Interpolator, FunctionSpace, FiniteElement,
+                       Projector, FunctionSpace, FiniteElement,
                        TensorProductElement)
 from firedrake.assign import Assigner
+from firedrake.__future__ import Interpolator
 from ufl.domain import extract_unique_domain
 
 from abc import ABCMeta, abstractmethod, abstractproperty

gusto/equations/compressible_euler_equations.py

@@ -40,7 +40,7 @@ def __init__(self, domain, parameters, sponge_options=None,
                  u_transport_option="vector_invariant_form",
                  diffusion_options=None,
                  no_normal_flow_bc_ids=None,
-                 active_tracers=None):
+                 active_tracers=None, max_quad_deg=5):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
@@ -75,6 +75,8 @@ def __init__(self, domain, parameters, sponge_options=None,
             active_tracers (list, optional): a list of `ActiveTracer` objects
                 that encode the metadata for any active tracers to be included
                 in the equations.. Defaults to None.
+            max_quad_deg (int): maximum quadrature degree for any form. Defaults
+                to 5.
 
         Raises:
             NotImplementedError: only mixing ratio tracers are implemented.
@@ -98,7 +100,8 @@ def __init__(self, domain, parameters, sponge_options=None,
         super().__init__(field_names, domain, space_names,
                          linearisation_map=linearisation_map,
                          no_normal_flow_bc_ids=no_normal_flow_bc_ids,
-                         active_tracers=active_tracers)
+                         active_tracers=active_tracers,
+                         max_quad_deg=max_quad_deg)
 
         self.parameters = parameters
         g = parameters.g
@@ -284,7 +287,8 @@ def __init__(self, domain, parameters, sponge=None,
                  u_transport_option="vector_invariant_form",
                  diffusion_options=None,
                  no_normal_flow_bc_ids=None,
-                 active_tracers=None):
+                 active_tracers=None,
+                 max_quad_deg=5):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
@@ -317,7 +321,9 @@ def __init__(self, domain, parameters, sponge=None,
                 None.
             active_tracers (list, optional): a list of `ActiveTracer` objects
                 that encode the metadata for any active tracers to be included
-                in the equations.. Defaults to None.
+                in the equations. Defaults to None.
+            max_quad_deg (int): maximum quadrature degree for any form. Defaults
+                to 5.
 
         Raises:
             NotImplementedError: only mixing ratio tracers are implemented.
@@ -329,7 +335,8 @@ def __init__(self, domain, parameters, sponge=None,
                          u_transport_option=u_transport_option,
                          diffusion_options=diffusion_options,
                          no_normal_flow_bc_ids=no_normal_flow_bc_ids,
-                         active_tracers=active_tracers)
+                         active_tracers=active_tracers,
+                         max_quad_deg=max_quad_deg)
 
         self.residual = self.residual.label_map(
             lambda t: t.has_label(time_derivative),

gusto/equations/prognostic_equations.py

@@ -24,20 +24,23 @@
 class PrognosticEquation(object, metaclass=ABCMeta):
     """Base class for prognostic equations."""
 
-    def __init__(self, domain, function_space, field_name):
+    def __init__(self, domain, function_space, field_name, max_quad_deg=5):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
                 mesh and the compatible function spaces.
             function_space (:class:`FunctionSpace`): the function space that the
                 equation's prognostic is defined on.
             field_name (str): name of the prognostic field.
+            max_quad_deg (int): maximum quadrature degree for any form. Defaults
+                to 5.
         """
 
         self.domain = domain
         self.function_space = function_space
         self.X = Function(function_space)
         self.field_name = field_name
+        self.max_quad_deg = max_quad_deg
         self.bcs = {}
         self.prescribed_fields = PrescribedFields()
 
@@ -76,7 +79,7 @@ class PrognosticEquationSet(PrognosticEquation, metaclass=ABCMeta):
 
     def __init__(self, field_names, domain, space_names,
                  linearisation_map=None, no_normal_flow_bc_ids=None,
-                 active_tracers=None):
+                 active_tracers=None, max_quad_deg=5):
         """
         Args:
             field_names (list): a list of strings for names of the prognostic
@@ -94,12 +97,15 @@ def __init__(self, field_names, domain, space_names,
             active_tracers (list, optional): a list of `ActiveTracer` objects
                 that encode the metadata for any active tracers to be included
                 in the equations.. Defaults to None.
+            max_quad_deg (int): maximum quadrature degree for any form. Defaults
+                to 5.
         """
 
         self.field_names = field_names
         self.space_names = space_names
         self.active_tracers = active_tracers
         self.linearisation_map = lambda t: False if linearisation_map is None else linearisation_map(t)
+        self.max_quad_deg = max_quad_deg
 
         # Build finite element spaces
         self.spaces = [domain.spaces(space_name) for space_name in

gusto/timestepping/timestepper.py

@@ -2,7 +2,7 @@
 
 from abc import ABCMeta, abstractmethod, abstractproperty
 from firedrake import Function, Projector, split
-from firedrake.fml import drop, Term
+from firedrake.fml import drop, Term, all_terms
 from pyop2.profiling import timed_stage
 from gusto.equations import PrognosticEquationSet
 from gusto.core import TimeLevelFields, StateFields
@@ -12,6 +12,7 @@
 from gusto.spatial_methods.transport_methods import TransportMethod
 import ufl
 
+
 __all__ = ["BaseTimestepper", "Timestepper", "PrescribedTransport"]
 
 
@@ -108,6 +109,39 @@ def setup_equation(self, equation):
             for method in self.spatial_methods:
                 method.replace_form(equation)
 
+        # -------------------------------------------------------------------- #
+        # Ensure the quadrature degree is not excessive for any integrals
+        # -------------------------------------------------------------------- #
+        def update_quadrature_degree(t):
+            # This function takes in a term and returns a new term
+            # with the same form and labels, the only difference being
+            # that any integrals in the form with no quadrature_degree
+            # set will have their quadrature degree set to max_quad_deg
+
+            # This list will hold the updated integrals
+            new_itgs = []
+
+            # Loop over integrals in this term's form
+            for itg in t.form._integrals:
+                # check if the quadrature degree is not set
+                if 'quadrature_degree' not in itg._metadata.keys():
+                    # create new integral with updated quadrature degree
+                    new_itg = itg.reconstruct(
+                        metadata={'quadrature_degree': equation.max_quad_deg})
+                    new_itgs.append(new_itg)
+                else:
+                    # if quadrature degree was already set, just keep
+                    # this integral
+                    new_itgs.append(itg)
+
+            new_form = ufl.form.Form(new_itgs)
+
+            return Term(new_form, t.labels)
+
+        equation.residual = equation.residual.label_map(
+            all_terms,
+            lambda t: update_quadrature_degree(t))
+
     def setup_transporting_velocity(self, scheme):
         """
         Set up the time discretisation by replacing the transporting velocity