Enable adjoint

examples/boussinesq/skamarock_klemp_compressible.py

@@ -66,7 +66,7 @@ def skamarock_klemp_compressible_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters(cs=cs)
+    parameters = BoussinesqParameters(mesh, cs=cs)
     eqns = BoussinesqEquations(domain, parameters)
 
     # I/O

examples/boussinesq/skamarock_klemp_incompressible.py

@@ -65,7 +65,7 @@ def skamarock_klemp_incompressible_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters()
+    parameters = BoussinesqParameters(mesh)
     eqns = BoussinesqEquations(domain, parameters, compressible=False)
 
     # I/O

examples/boussinesq/skamarock_klemp_linear.py

@@ -63,7 +63,7 @@ def skamarock_klemp_linear_bouss(
     domain = Domain(mesh, dt, 'CG', element_order)
 
     # Equation
-    parameters = BoussinesqParameters(cs=cs)
+    parameters = BoussinesqParameters(mesh, cs=cs)
     eqns = LinearBoussinesqEquations(domain, parameters)
 
     # I/O

examples/compressible_euler/dcmip_3_1_gravity_wave.py

@@ -43,16 +43,10 @@ def dcmip_3_1_gravity_wave(
     # Test case parameters
     # ------------------------------------------------------------------------ #
 
-    parameters = CompressibleParameters()
     a_ref = 6.37122e6               # Radius of the Earth (m)
     X = 125.0                       # Reduced-size Earth reduction factor
     a = a_ref/X                     # Scaled radius of planet (m)
-    g = parameters.g                # Acceleration due to gravity (m/s^2)
     N = 0.01                        # Brunt-Vaisala frequency (1/s)
-    p_0 = parameters.p_0            # Reference pressure (Pa, not hPa)
-    c_p = parameters.cp             # SHC of dry air at constant pressure (J/kg/K)
-    R_d = parameters.R_d            # Gas constant for dry air (J/kg/K)
-    kappa = parameters.kappa        # R_d/c_p
     T_eq = 300.0                    # Isothermal atmospheric temperature (K)
     p_eq = 1000.0 * 100.0           # Reference surface pressure at the equator
     u_max = 20.0                    # Maximum amplitude of the zonal wind (m/s)
@@ -83,6 +77,7 @@ def dcmip_3_1_gravity_wave(
     domain = Domain(mesh, dt, "RTCF", element_order)
 
     # Equation
+    parameters = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, parameters, u_transport_option=u_eqn_type
     )
@@ -137,6 +132,13 @@ def dcmip_3_1_gravity_wave(
     # Initial conditions with u0
     uexpr = as_vector([-u_max*y/a, u_max*x/a, 0.0])
 
+    # Parameters required for initialising
+    g = parameters.g                # Acceleration due to gravity (m/s^2)
+    p_0 = parameters.p_0            # Reference pressure (Pa, not hPa)
+    c_p = parameters.cp             # SHC of dry air at constant pressure (J/kg/K)
+    R_d = parameters.R_d            # Gas constant for dry air (J/kg/K)
+    kappa = parameters.kappa        # R_d/c_p
+
     # Surface temperature
     G = g**2/(N**2*c_p)
     Ts_expr = (

examples/compressible_euler/dry_bryan_fritsch.py

@@ -66,7 +66,7 @@ def dry_bryan_fritsch(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     eqns = CompressibleEulerEquations(
         domain, params, u_transport_option=u_eqn_type,
         no_normal_flow_bc_ids=[1, 2]

examples/compressible_euler/skamarock_klemp_nonhydrostatic.py

@@ -78,7 +78,7 @@ def skamarock_klemp_nonhydrostatic(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    parameters = CompressibleParameters()
+    parameters = CompressibleParameters(mesh)
     if hydrostatic:
         eqns = HydrostaticCompressibleEulerEquations(domain, parameters)
     else:

examples/compressible_euler/straka_bubble.py

@@ -70,8 +70,8 @@ def straka_bubble(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    parameters = CompressibleParameters()
-    diffusion_params = DiffusionParameters(kappa=kappa, mu=mu0/delta)
+    parameters = CompressibleParameters(mesh)
+    diffusion_params = DiffusionParameters(mesh, kappa=kappa, mu=mu0/delta)
     diffusion_options = [("u", diffusion_params), ("theta", diffusion_params)]
     eqns = CompressibleEulerEquations(
         domain, parameters, u_transport_option=u_eqn_type,

examples/compressible_euler/unsaturated_bubble.py

@@ -84,7 +84,7 @@ def unsaturated_bubble(
     domain = Domain(mesh, dt, "CG", element_order)
 
     # Equation
-    params = CompressibleParameters()
+    params = CompressibleParameters(mesh)
     tracers = [WaterVapour(), CloudWater(), Rain()]
     eqns = CompressibleEulerEquations(
         domain, params, active_tracers=tracers, u_transport_option=u_eqn_type

examples/shallow_water/linear_thermal_galewsky_jet.py

@@ -64,7 +64,7 @@ def linear_thermal_galewsky_jet(
     domain = Domain(mesh, dt, 'BDM', element_order)
 
     # Equation
-    parameters = ShallowWaterParameters(H=H)
+    parameters = ShallowWaterParameters(mesh, H=H)
     Omega = parameters.Omega
     fexpr = 2*Omega*xyz[2]/R
     eqns = LinearThermalShallowWaterEquations(domain, parameters, fexpr=fexpr)

examples/shallow_water/linear_williamson_2.py

@@ -57,7 +57,7 @@ def linear_williamson_2(
     domain = Domain(mesh, dt, 'BDM', element_order)
 
     # Equation
-    parameters = ShallowWaterParameters(H=mean_depth)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
     eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr)

examples/shallow_water/moist_convective_williamson_2.py

@@ -76,7 +76,7 @@ def moist_convect_williamson_2(
     _, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equations
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 

examples/shallow_water/moist_thermal_equivb_gw.py

@@ -59,7 +59,8 @@ def moist_thermal_gw(
     xyz = SpatialCoordinate(mesh)
 
     # Equation parameters
-    parameters = ShallowWaterParameters(H=mean_depth, q0=q0, nu=nu, beta2=beta2)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, q0=q0,
+                                        nu=nu, beta2=beta2)
     Omega = parameters.Omega
     fexpr = 2*Omega*xyz[2]/radius
 

examples/shallow_water/moist_thermal_williamson_5.py

@@ -86,7 +86,7 @@ def moist_thermal_williamson_5(
     lamda, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equation: coriolis
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 

examples/shallow_water/shallow_water_1d_wave.py

@@ -62,17 +62,17 @@ def shallow_water_1d_wave(
 
     # Diffusion
     delta = domain_length / ncells
-    u_diffusion_opts = DiffusionParameters(kappa=kappa)
-    v_diffusion_opts = DiffusionParameters(kappa=kappa, mu=10/delta)
-    D_diffusion_opts = DiffusionParameters(kappa=kappa, mu=10/delta)
+    u_diffusion_opts = DiffusionParameters(mesh, kappa=kappa)
+    v_diffusion_opts = DiffusionParameters(mesh, kappa=kappa, mu=10/delta)
+    D_diffusion_opts = DiffusionParameters(mesh, kappa=kappa, mu=10/delta)
     diffusion_options = [
         ("u", u_diffusion_opts),
         ("v", v_diffusion_opts),
         ("D", D_diffusion_opts)
     ]
 
     # Equation
-    parameters = ShallowWaterParameters(H=1/epsilon, g=1/epsilon)
+    parameters = ShallowWaterParameters(mesh, H=1/epsilon, g=1/epsilon)
     eqns = ShallowWaterEquations_1d(
         domain, parameters, fexpr=Constant(1/epsilon),
         diffusion_options=diffusion_options

examples/shallow_water/thermal_williamson_2.py

@@ -68,7 +68,7 @@ def thermal_williamson_2(
     x, y, z = SpatialCoordinate(mesh)
 
     # Equations
-    params = ShallowWaterParameters(H=mean_depth, g=g)
+    params = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = params.Omega
     fexpr = 2*Omega*z/radius
     eqns = ThermalShallowWaterEquations(

examples/shallow_water/williamson_2.py

@@ -62,7 +62,7 @@ def williamson_2(
     xyz = SpatialCoordinate(mesh)
 
     # Equation
-    parameters = ShallowWaterParameters(H=mean_depth)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth)
     Omega = parameters.Omega
     _, lat, _ = rotated_lonlatr_coords(xyz, rotate_pole_to)
     e_lon, _, _ = rotated_lonlatr_vectors(xyz, rotate_pole_to)

examples/shallow_water/williamson_5.py

@@ -65,7 +65,7 @@ def williamson_5(
     lamda, phi, _ = lonlatr_from_xyz(x, y, z)
 
     # Equation: coriolis
-    parameters = ShallowWaterParameters(H=mean_depth, g=g)
+    parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g)
     Omega = parameters.Omega
     fexpr = 2*Omega*z/radius
 

gusto/core/__init__.py

@@ -3,10 +3,11 @@
 from gusto.core.coordinates import *               # noqa
 from gusto.core.coord_transforms import *          # noqa
 from gusto.core.domain import *                    # noqa
+from gusto.core.equation_configuration import *    # noqa
 from gusto.core.fields import *                    # noqa
 from gusto.core.function_spaces import *           # noqa
 from gusto.core.io import *                        # noqa
 from gusto.core.kernels import *                   # noqa
 from gusto.core.labels import *                    # noqa
 from gusto.core.logging import *                   # noqa
-from gusto.core.meshes import *                    # noqa
+from gusto.core.meshes import *                    # noqa

gusto/core/configuration.py

@@ -1,17 +1,14 @@
 """Some simple tools for configuring the model."""
 from abc import ABCMeta, abstractproperty
 from enum import Enum
-from firedrake import sqrt, Constant
+from firedrake import sqrt
 
 
 __all__ = [
     "Configuration",
     "IntegrateByParts", "TransportEquationType", "OutputParameters",
-    "BoussinesqParameters", "CompressibleParameters",
-    "ShallowWaterParameters",
     "EmbeddedDGOptions", "ConservativeEmbeddedDGOptions", "RecoveryOptions",
     "ConservativeRecoveryOptions", "SUPGOptions", "MixedFSOptions",
-    "SpongeLayerParameters", "DiffusionParameters", "BoundaryLayerParameters",
     "SubcyclingOptions"
 ]
 
@@ -64,7 +61,7 @@ def __setattr__(self, name, value):
 
         When attributes are provided as floats or integers, these are converted
         to Firedrake :class:`Constant` objects, other than a handful of special
-        integers (dumpfreq, pddumpfreq, chkptfreq and log_level).
+        integers.
 
         Args:
             name: the attribute's name.
@@ -75,18 +72,19 @@ def __setattr__(self, name, value):
                 this attribute pre-defined.
         """
         if not hasattr(self, name):
-            raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name))
+            raise AttributeError(f"{type(self).__name__} object has no attribute {name}.")
 
-        # Almost all parameters should be Constants -- but there are some
-        # specific exceptions which should be kept as integers
+        # Almost all parameters should be functions on the real space
+        # -- but there are some specific exceptions which should be
+        # kept as integers
         non_constants = [
             'dumpfreq', 'pddumpfreq', 'chkptfreq',
             'fixed_subcycles', 'max_subcycles', 'subcycle_by_courant'
         ]
         if type(value) in [float, int] and name not in non_constants:
-            object.__setattr__(self, name, Constant(value))
-        else:
-            object.__setattr__(self, name, value)
+            raise AttributeError(f"Attribute {name} requires a mesh.")
+
+        object.__setattr__(self, name, value)
 
 
 class OutputParameters(Configuration):
@@ -115,55 +113,6 @@ class OutputParameters(Configuration):
     tolerance = None
 
 
-class BoussinesqParameters(Configuration):
-    """Physical parameters for the Boussinesq equations."""
-
-    g = 9.810616
-    N = 0.01  # Brunt-Vaisala frequency (1/s)
-    cs = 340  # sound speed (for compressible case) (m/s)
-    Omega = None
-
-
-class CompressibleParameters(Configuration):
-    """Physical parameters for the Compressible Euler equations."""
-
-    g = 9.810616
-    N = 0.01  # Brunt-Vaisala frequency (1/s)
-    cp = 1004.5  # SHC of dry air at const. pressure (J/kg/K)
-    R_d = 287.  # Gas constant for dry air (J/kg/K)
-    kappa = 2.0/7.0  # R_d/c_p
-    p_0 = 1000.0*100.0  # reference pressure (Pa, not hPa)
-    cv = 717.5  # SHC of dry air at const. volume (J/kg/K)
-    c_pl = 4186.  # SHC of liq. wat. at const. pressure (J/kg/K)
-    c_pv = 1885.  # SHC of wat. vap. at const. pressure (J/kg/K)
-    c_vv = 1424.  # SHC of wat. vap. at const. pressure (J/kg/K)
-    R_v = 461.  # gas constant of water vapour
-    L_v0 = 2.5e6  # ref. value for latent heat of vap. (J/kg)
-    T_0 = 273.15  # ref. temperature
-    w_sat1 = 380.3  # first const. in Teten's formula (Pa)
-    w_sat2 = -17.27  # second const. in Teten's formula (no units)
-    w_sat3 = 35.86  # third const. in Teten's formula (K)
-    w_sat4 = 610.9  # fourth const. in Teten's formula (Pa)
-    Omega = None    # Rotation rate
-
-
-class ShallowWaterParameters(Configuration):
-    """Physical parameters for the shallow-water equations."""
-
-    g = 9.80616
-    Omega = 7.292e-5  # rotation rate
-    H = None  # mean depth
-    # Factor that multiplies the vapour in the equivalent buoyancy
-    # formulation of the thermal shallow water equations
-    beta2 = None
-    # Scaling factor for the saturation function in the equivalent buoyancy
-    # formulation of the thermal shallow water equations
-    nu = None
-    # Scaling factor for the saturation function in the equivalent buoyancy
-    # formulation of the thermal shallow water equations
-    q0 = None
-
-
 class WrapperOptions(Configuration, metaclass=ABCMeta):
     """Base class for specifying options for a transport scheme."""
 
@@ -239,49 +188,6 @@ class MixedFSOptions(WrapperOptions):
     suboptions = None
 
 
-class SpongeLayerParameters(Configuration):
-    """Specifies parameters describing a 'sponge' (damping) layer."""
-
-    H = None
-    z_level = None
-    mubar = None
-
-
-class DiffusionParameters(Configuration):
-    """Parameters for a diffusion term with an interior penalty method."""
-
-    kappa = None
-    mu = None
-
-
-class BoundaryLayerParameters(Configuration):
-    """
-    Parameters for the idealised wind drag, surface flux and boundary layer
-    mixing schemes.
-    """
-
-    coeff_drag_0 = 7e-4         # Zeroth drag coefficient (dimensionless)
-    coeff_drag_1 = 6.5e-5       # First drag coefficient (s/m)
-    coeff_drag_2 = 2e-3         # Second drag coefficient (dimensionless)
-    coeff_heat = 1.1e-3         # Dimensionless surface sensible heat coefficient
-    coeff_evap = 1.1e-3         # Dimensionless surface evaporation coefficient
-    height_surface_layer = 75.  # Height (m) of surface level (usually lowest level)
-    mu = 100.                   # Interior penalty coefficient for vertical diffusion
-
-
-class HeldSuarezParameters(Configuration):
-    """
-    Parameters used in the default configuration for the Held Suarez test case.
-    """
-    T0stra = 200               # Stratosphere temp
-    T0surf = 315               # Surface temperature at equator
-    T0horiz = 60               # Equator to pole temperature difference
-    T0vert = 10                # Stability parameter
-    sigmab = 0.7               # Height of the boundary layer
-    tau_d = 40 * 24 * 60 * 60  # 40 day time scale
-    tau_u = 4 * 24 * 60 * 60   # 4 day timescale
-
-
 class SubcyclingOptions(Configuration):
     """
     Describes the process of subcycling a time discretisation, by dividing the

gusto/core/conservative_projection.py

@@ -68,10 +68,11 @@ def __init__(self, rho_source, rho_target, m_source, m_target,
         self.m_target = m_target
 
         V = self.m_target.function_space()
-        mesh = V.mesh()
 
-        self.m_mean = Constant(0.0, domain=mesh)
-        self.volume = assemble(Constant(1.0, domain=mesh)*dx)
+        self.m_mean = Constant(0.0)
+        self.volume = assemble(
+            1*dx(domain=ufl.domain.extract_unique_domain(self.m_target))
+        )
 
         test = TestFunction(V)
         m_trial = TrialFunction(V)