diff --git a/examples/boussinesq/skamarock_klemp_compressible.py b/examples/boussinesq/skamarock_klemp_compressible.py index c15daff84..5fb9ce075 100644 --- a/examples/boussinesq/skamarock_klemp_compressible.py +++ b/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 diff --git a/examples/boussinesq/skamarock_klemp_incompressible.py b/examples/boussinesq/skamarock_klemp_incompressible.py index b7bb5fd5e..8aa3d7265 100644 --- a/examples/boussinesq/skamarock_klemp_incompressible.py +++ b/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 diff --git a/examples/boussinesq/skamarock_klemp_linear.py b/examples/boussinesq/skamarock_klemp_linear.py index e8145548d..2d430b58b 100644 --- a/examples/boussinesq/skamarock_klemp_linear.py +++ b/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 diff --git a/examples/compressible_euler/dcmip_3_1_gravity_wave.py b/examples/compressible_euler/dcmip_3_1_gravity_wave.py index 6862dd94f..5953eeed9 100644 --- a/examples/compressible_euler/dcmip_3_1_gravity_wave.py +++ b/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 = ( diff --git a/examples/compressible_euler/dry_bryan_fritsch.py b/examples/compressible_euler/dry_bryan_fritsch.py index dc4c92e9e..380f182fd 100644 --- a/examples/compressible_euler/dry_bryan_fritsch.py +++ b/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] diff --git a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py b/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py index 5ec8ab12b..5ce8e5e49 100644 --- a/examples/compressible_euler/skamarock_klemp_nonhydrostatic.py +++ b/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: diff --git a/examples/compressible_euler/straka_bubble.py b/examples/compressible_euler/straka_bubble.py index fed55a1ec..d33eef5d5 100644 --- a/examples/compressible_euler/straka_bubble.py +++ b/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, diff --git a/examples/compressible_euler/unsaturated_bubble.py b/examples/compressible_euler/unsaturated_bubble.py index 59a911cb9..f7af29759 100644 --- a/examples/compressible_euler/unsaturated_bubble.py +++ b/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 diff --git a/examples/shallow_water/linear_thermal_galewsky_jet.py b/examples/shallow_water/linear_thermal_galewsky_jet.py index a4fd68631..a5622c055 100644 --- a/examples/shallow_water/linear_thermal_galewsky_jet.py +++ b/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) diff --git a/examples/shallow_water/linear_williamson_2.py b/examples/shallow_water/linear_williamson_2.py index 9d022bf11..a50b68697 100644 --- a/examples/shallow_water/linear_williamson_2.py +++ b/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) diff --git a/examples/shallow_water/moist_convective_williamson_2.py b/examples/shallow_water/moist_convective_williamson_2.py index 376a76c64..63038b9f4 100644 --- a/examples/shallow_water/moist_convective_williamson_2.py +++ b/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 diff --git a/examples/shallow_water/moist_thermal_equivb_gw.py b/examples/shallow_water/moist_thermal_equivb_gw.py index 912e9deba..f2928461c 100644 --- a/examples/shallow_water/moist_thermal_equivb_gw.py +++ b/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 diff --git a/examples/shallow_water/moist_thermal_williamson_5.py b/examples/shallow_water/moist_thermal_williamson_5.py index 376cc94ad..625ee57dc 100644 --- a/examples/shallow_water/moist_thermal_williamson_5.py +++ b/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 diff --git a/examples/shallow_water/shallow_water_1d_wave.py b/examples/shallow_water/shallow_water_1d_wave.py index 3e0b1582b..057a1f157 100644 --- a/examples/shallow_water/shallow_water_1d_wave.py +++ b/examples/shallow_water/shallow_water_1d_wave.py @@ -62,9 +62,9 @@ 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), @@ -72,7 +72,7 @@ def shallow_water_1d_wave( ] # 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 diff --git a/examples/shallow_water/thermal_williamson_2.py b/examples/shallow_water/thermal_williamson_2.py index 7adc305dc..318f44962 100644 --- a/examples/shallow_water/thermal_williamson_2.py +++ b/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( diff --git a/examples/shallow_water/williamson_2.py b/examples/shallow_water/williamson_2.py index 8c5f60917..f362605bd 100644 --- a/examples/shallow_water/williamson_2.py +++ b/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) diff --git a/examples/shallow_water/williamson_5.py b/examples/shallow_water/williamson_5.py index efe41b183..13f48a691 100644 --- a/examples/shallow_water/williamson_5.py +++ b/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 diff --git a/gusto/core/__init__.py b/gusto/core/__init__.py index 432c66b9d..3cebab236 100644 --- a/gusto/core/__init__.py +++ b/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 \ No newline at end of file +from gusto.core.meshes import * # noqa diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py index 40852033a..0a542d8cb 100644 --- a/gusto/core/configuration.py +++ b/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 diff --git a/gusto/core/conservative_projection.py b/gusto/core/conservative_projection.py index 1ab1455d6..8397cc67c 100644 --- a/gusto/core/conservative_projection.py +++ b/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) diff --git a/gusto/core/domain.py b/gusto/core/domain.py index 27addd9fa..2caedfab3 100644 --- a/gusto/core/domain.py +++ b/gusto/core/domain.py @@ -65,15 +65,16 @@ def __init__(self, mesh, dt, family, degree=None, # -------------------------------------------------------------------- # # Store central dt for use in the rest of the model + R = FunctionSpace(mesh, "R", 0) if type(dt) is Constant: - self.dt = dt + self.dt = Function(R, val=float(dt)) elif type(dt) in (float, int): - self.dt = Constant(dt) + self.dt = Function(R, val=dt) else: raise TypeError(f'dt must be a Constant, float or int, not {type(dt)}') # Make a placeholder for the time - self.t = Constant(0.0) + self.t = Function(R, val=0.0) # -------------------------------------------------------------------- # # Build compatible function spaces diff --git a/gusto/core/equation_configuration.py b/gusto/core/equation_configuration.py new file mode 100644 index 000000000..d8dca705b --- /dev/null +++ b/gusto/core/equation_configuration.py @@ -0,0 +1,152 @@ +"""Some simple tools for configuring the model.""" +from firedrake import Function, FunctionSpace, Constant +import inspect + + +__all__ = [ + "BoussinesqParameters", "CompressibleParameters", + "ShallowWaterParameters", + "SpongeLayerParameters", "DiffusionParameters", "BoundaryLayerParameters", +] + + +class EquationParameters(object): + """A base configuration object for storing equation parameters.""" + + mesh = None + + def __init__(self, mesh, **kwargs): + """ + Args: + mesh: for creating the real function space + **kwargs: attributes and their values to be stored in the object. + """ + self.mesh = mesh + typecheck = lambda val: type(val) in [float, int, Constant] + params = dict(inspect.getmembers(self, typecheck)) + params.update(kwargs.items()) + for name, value in params.items(): + self.__setattr__(name, value) + + def __setattr__(self, name, value): + """ + Sets the model configuration attributes. + + When attributes are provided as floats or integers, these are converted + to Firedrake :class:`Constant` objects, other than a handful of special + integers. + + Args: + name: the attribute's name. + value: the value to provide to the attribute. + + Raises: + AttributeError: if the :class:`Configuration` object does not have + this attribute pre-defined. + """ + if not hasattr(self, name): + raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name)) + + # Almost all parameters should be functions on the real space + # -- but there are some specific exceptions which should be + # kept as integers + if self.mesh is not None: + # This check is required so that on instantiation we do + # not hit this line while self.mesh is still None + R = FunctionSpace(self.mesh, 'R', 0) + if type(value) in [float, int, Constant]: + object.__setattr__(self, name, Function(R, val=float(value))) + else: + object.__setattr__(self, name, value) + + +class BoussinesqParameters(EquationParameters): + """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(EquationParameters): + """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(EquationParameters): + """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 SpongeLayerParameters(EquationParameters): + """Specifies parameters describing a 'sponge' (damping) layer.""" + + H = None + z_level = None + mubar = None + + +class DiffusionParameters(EquationParameters): + """Parameters for a diffusion term with an interior penalty method.""" + + kappa = None + mu = None + + +class BoundaryLayerParameters(EquationParameters): + """ + 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(EquationParameters): + """ + 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 diff --git a/gusto/core/io.py b/gusto/core/io.py index bd80f86d5..1e69a8651 100644 --- a/gusto/core/io.py +++ b/gusto/core/io.py @@ -8,7 +8,8 @@ from gusto.diagnostics import Diagnostics, CourantNumber from gusto.core.meshes import get_flat_latlon_mesh from firedrake import (Function, functionspaceimpl, Constant, COMM_WORLD, - DumbCheckpoint, FILE_CREATE, FILE_READ, CheckpointFile) + DumbCheckpoint, FILE_CREATE, FILE_READ, CheckpointFile, + MeshGeometry) from firedrake.output import VTKFile from pyop2.mpi import MPI import numpy as np @@ -257,7 +258,7 @@ def log_parameters(self, equation): """ if hasattr(equation, 'parameters') and equation.parameters is not None: logger.info("Physical parameters that take non-default values:") - logger.info(", ".join("%s: %s" % (k, float(v)) for (k, v) in vars(equation.parameters).items())) + logger.info(", ".join("%s: %s" % (k, float(v)) for (k, v) in vars(equation.parameters).items() if type(v) is not MeshGeometry)) def setup_log_courant(self, state_fields, name='u', component="whole", expression=None): diff --git a/gusto/equations/boussinesq_equations.py b/gusto/equations/boussinesq_equations.py index b9823f110..d94e2a038 100644 --- a/gusto/equations/boussinesq_equations.py +++ b/gusto/equations/boussinesq_equations.py @@ -168,10 +168,12 @@ def __init__(self, domain, parameters, # -------------------------------------------------------------------- # if compressible: cs = parameters.cs + # On assuming ``cs`` as a constant, it is right keep it out of the + # integration. linear_div_form = divergence(subject( - prognostic(cs**2 * phi * div(u_trial) * dx, 'p'), self.X)) + prognostic(cs**2 * (phi * div(u_trial) * dx), 'p'), self.X)) divergence_form = divergence(linearisation( - subject(prognostic(cs**2 * phi * div(u) * dx, 'p'), self.X), + subject(prognostic(cs**2 * (phi * div(u) * dx), 'p'), self.X), linear_div_form)) else: # This enforces that div(u) = 0 diff --git a/gusto/equations/compressible_euler_equations.py b/gusto/equations/compressible_euler_equations.py index 73bc76d45..d512432cf 100644 --- a/gusto/equations/compressible_euler_equations.py +++ b/gusto/equations/compressible_euler_equations.py @@ -17,7 +17,6 @@ ) from gusto.equations.active_tracers import Phases, TracerVariableType from gusto.equations.prognostic_equations import PrognosticEquationSet - __all__ = ["CompressibleEulerEquations", "HydrostaticCompressibleEulerEquations"] @@ -45,7 +44,7 @@ def __init__(self, domain, parameters, sponge_options=None, Args: domain (:class:`Domain`): the model's domain object, containing the mesh and the compatible function spaces. - parameters (:class:`Configuration`, optional): an object containing + x (:class:`Configuration`, optional): an object containing the model's physical parameters. sponge_options (:class:`SpongeLayerParameters`, optional): any parameters for applying a sponge layer to the upper boundary. diff --git a/gusto/equations/shallow_water_equations.py b/gusto/equations/shallow_water_equations.py index fd61cedf7..0e9aebe2a 100644 --- a/gusto/equations/shallow_water_equations.py +++ b/gusto/equations/shallow_water_equations.py @@ -163,8 +163,9 @@ def _setup_residual(self, fexpr, topog_expr, u_transport_option): # -------------------------------------------------------------------- # # Pressure Gradient Term # -------------------------------------------------------------------- # + # On assuming ``g``, it is right to keep it out of the integral. pressure_gradient_form = pressure_gradient( - subject(prognostic(-g*div(w)*D*dx, 'u'), self.X)) + subject(prognostic(-g*(div(w)*D*dx), 'u'), self.X)) residual = (mass_form + adv_form + pressure_gradient_form) @@ -418,6 +419,7 @@ def _setup_residual(self, fexpr, topog_expr, u_transport_option): # provide linearisation if self.equivalent_buoyancy: beta2 = self.parameters.beta2 + qsat_expr = self.compute_saturation(self.X) qv = conditional(qt < qsat_expr, qt, qsat_expr) qvbar = conditional(qtbar < qsat_expr, qtbar, qsat_expr) diff --git a/gusto/physics/boundary_and_turbulence.py b/gusto/physics/boundary_and_turbulence.py index dc19c7d75..ea4575c2e 100644 --- a/gusto/physics/boundary_and_turbulence.py +++ b/gusto/physics/boundary_and_turbulence.py @@ -8,7 +8,7 @@ SpatialCoordinate, dS_v ) from firedrake.fml import subject -from gusto.core.configuration import BoundaryLayerParameters +from gusto.core.equation_configuration import BoundaryLayerParameters from gusto.recovery import Recoverer, BoundaryMethod from gusto.equations import CompressibleEulerEquations from gusto.core.labels import prognostic, source_label diff --git a/gusto/physics/held_suarez_forcing.py b/gusto/physics/held_suarez_forcing.py index d7573d009..e572eb578 100644 --- a/gusto/physics/held_suarez_forcing.py +++ b/gusto/physics/held_suarez_forcing.py @@ -7,7 +7,7 @@ from gusto.physics.physics_parametrisation import PhysicsParametrisation from gusto.core.labels import prognostic from gusto.equations import thermodynamics -from gusto.core.configuration import HeldSuarezParameters +from gusto.core.equation_configuration import HeldSuarezParameters from gusto.core import logger @@ -26,7 +26,7 @@ def __init__(self, equation, variable_name, parameters, hs_parameters=None): """ label_name = f'relaxation_{variable_name}' if hs_parameters is None: - hs_parameters = HeldSuarezParameters() + hs_parameters = HeldSuarezParameters(equation.domain.mesh) logger.warning('Using default Held-Suarez parameters') super().__init__(equation, label_name, hs_parameters) @@ -121,7 +121,7 @@ def __init__(self, equation, hs_parameters=None): """ label_name = 'rayleigh_friction' if hs_parameters is None: - hs_parameters = HeldSuarezParameters() + hs_parameters = HeldSuarezParameters(equation.domain.mesh) logger.warning('Using default Held-Suarez parameters') super().__init__(equation, label_name, hs_parameters) diff --git a/gusto/physics/shallow_water_microphysics.py b/gusto/physics/shallow_water_microphysics.py index 29b3637bd..f30685fa0 100644 --- a/gusto/physics/shallow_water_microphysics.py +++ b/gusto/physics/shallow_water_microphysics.py @@ -3,7 +3,8 @@ """ from firedrake import ( - conditional, Function, dx, min_value, max_value, Constant, assemble, split + conditional, Function, dx, min_value, max_value, FunctionSpace, + assemble, split ) from firedrake.__future__ import interpolate from firedrake.fml import subject @@ -97,13 +98,14 @@ def __init__(self, equation, saturation_curve, self.source_expr = split(self.source)[self.Vv_idx] self.source_int = self.source.subfunctions[self.Vv_idx] + R = FunctionSpace(equation.domain.mesh, "R", 0) # tau is the timescale for conversion (may or may not be the timestep) if tau is not None: self.set_tau_to_dt = False - self.tau = tau + self.tau = Function(R).assign(tau) else: self.set_tau_to_dt = True - self.tau = Constant(0) + self.tau = Function(R) logger.info("Timescale for rain conversion has been set to dt. If this is not the intention then provide a tau parameter as an argument to InstantRain.") if self.time_varying_saturation: @@ -280,12 +282,13 @@ def __init__(self, equation, saturation_curve, V_idxs.append(self.Vb_idx) # tau is the timescale for condensation/evaporation (may or may not be the timestep) + R = FunctionSpace(equation.domain.mesh, "R", 0) if tau is not None: self.set_tau_to_dt = False - self.tau = tau + self.tau = Function(R).assign(tau) else: self.set_tau_to_dt = True - self.tau = Constant(0) + self.tau = Function(R) logger.info("Timescale for moisture conversion between vapour and cloud has been set to dt. If this is not the intention then provide a tau parameter as an argument to SWSaturationAdjustment.") if self.time_varying_saturation: diff --git a/gusto/spatial_methods/diffusion_methods.py b/gusto/spatial_methods/diffusion_methods.py index 768b7532e..bc33da4df 100644 --- a/gusto/spatial_methods/diffusion_methods.py +++ b/gusto/spatial_methods/diffusion_methods.py @@ -161,4 +161,5 @@ def __init__(self, equation, variable, diffusion_parameters): kappa = diffusion_parameters.kappa self.form = diffusion(kappa * self.test.dx(0) * self.field.dx(0) * dx) else: - raise NotImplementedError("CG diffusion only implemented in 1D") + kappa = diffusion_parameters.kappa + self.form = diffusion(kappa * inner(grad(self.test), grad(self.field)) * dx) diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py index db71412ca..7b5eee9a6 100644 --- a/gusto/time_discretisation/time_discretisation.py +++ b/gusto/time_discretisation/time_discretisation.py @@ -9,8 +9,8 @@ import math from firedrake import (Function, TestFunction, TestFunctions, DirichletBC, - Constant, NonlinearVariationalProblem, - NonlinearVariationalSolver) + NonlinearVariationalProblem, NonlinearVariationalSolver, + FunctionSpace) from firedrake.fml import (replace_subject, replace_test_function, Term, all_terms, drop) from firedrake.formmanipulation import split_form @@ -88,10 +88,10 @@ def __init__(self, domain, field_name=None, subcycling_options=None, self.domain = domain self.field_name = field_name self.equation = None - - self.dt = Constant(0.0) + R = FunctionSpace(domain.mesh, "R", 0) + self.dt = Function(R, val=0.0) self.dt.assign(domain.dt) - self.original_dt = Constant(0.0) + self.original_dt = Function(R, val=0.0) self.original_dt.assign(self.dt) self.options = options self.limiter = limiter diff --git a/gusto/timestepping/semi_implicit_quasi_newton.py b/gusto/timestepping/semi_implicit_quasi_newton.py index 51da84f2c..c0ed20fb1 100644 --- a/gusto/timestepping/semi_implicit_quasi_newton.py +++ b/gusto/timestepping/semi_implicit_quasi_newton.py @@ -5,7 +5,7 @@ from firedrake import ( Function, Constant, TrialFunctions, DirichletBC, div, assemble, - LinearVariationalProblem, LinearVariationalSolver + LinearVariationalProblem, LinearVariationalSolver, FunctionSpace ) from firedrake.fml import drop, replace_subject from firedrake.__future__ import interpolate @@ -37,7 +37,7 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, auxiliary_equations_and_schemes=None, linear_solver=None, diffusion_schemes=None, physics_schemes=None, slow_physics_schemes=None, fast_physics_schemes=None, - alpha=Constant(0.5), off_centred_u=False, + alpha=0.5, off_centred_u=False, num_outer=2, num_inner=2, accelerator=False, predictor=None, reference_update_freq=None, spinup_steps=0): @@ -73,9 +73,9 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, (:class:`PhysicsParametrisation`, :class:`TimeDiscretisation`). These schemes are evaluated within the outer loop. Defaults to None. - alpha (`ufl.Constant`, optional): the semi-implicit off-centering + alpha (`float, optional): the semi-implicit off-centering parameter. A value of 1 corresponds to fully implicit, while 0 - corresponds to fully explicit. Defaults to Constant(0.5). + corresponds to fully explicit. Defaults to 0.5. off_centred_u (bool, optional): option to offcentre the transporting velocity. Defaults to False, in which case transporting velocity is centred. If True offcentring uses value of alpha. @@ -114,12 +114,14 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, self.num_outer = num_outer self.num_inner = num_inner - self.alpha = Constant(alpha) + mesh = equation_set.domain.mesh + R = FunctionSpace(mesh, "R", 0) + self.alpha = Function(R, val=float(alpha)) self.predictor = predictor self.accelerator = accelerator # Options relating to reference profiles and spin-up - self._alpha_original = Constant(alpha) + self._alpha_original = float(alpha) # float so as to not upset adjoint self.reference_update_freq = reference_update_freq self.to_update_ref_profile = False self.spinup_steps = spinup_steps @@ -131,7 +133,7 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods, # default is to not offcentre transporting velocity but if it # is offcentred then use the same value as alpha - self.alpha_u = Constant(alpha) if off_centred_u else Constant(0.5) + self.alpha_u = Function(R, val=float(alpha)) if off_centred_u else Function(R, val=0.5) self.spatial_methods = spatial_methods @@ -539,7 +541,7 @@ def __init__(self, equation, alpha): Args: equation (:class:`PrognosticEquationSet`): the prognostic equations containing the forcing terms. - alpha (:class:`Constant`): semi-implicit off-centering factor. An + alpha (:class:`Function`): semi-implicit off-centering factor. An alpha of 0 corresponds to fully explicit, while a factor of 1 corresponds to fully implicit. """ @@ -568,7 +570,8 @@ def __init__(self, equation, alpha): map_if_false=drop) # the explicit forms are multiplied by (1-alpha) and moved to the rhs - L_explicit = -(Constant(1)-alpha)*dt*residual.label_map( + one_minus_alpha = Function(alpha.function_space(), val=1-alpha) + L_explicit = -one_minus_alpha*dt*residual.label_map( lambda t: any(t.has_label(time_derivative, hydrostatic, *implicit_terms, return_tuple=True)), diff --git a/gusto/timestepping/timestepper.py b/gusto/timestepping/timestepper.py index 24a2058f7..f4f5bba0a 100644 --- a/gusto/timestepping/timestepper.py +++ b/gusto/timestepping/timestepper.py @@ -235,7 +235,7 @@ def run(self, t, tmax, pick_up=False): self.timestep() - self.t.assign(self.t + self.dt) + self.t.assign(float(self.t) + float(self.dt)) self.step += 1 with timed_stage("Dump output"): diff --git a/integration-tests/adjoints/test_diffusion_sensitivity.py b/integration-tests/adjoints/test_diffusion_sensitivity.py new file mode 100644 index 000000000..7c353eef9 --- /dev/null +++ b/integration-tests/adjoints/test_diffusion_sensitivity.py @@ -0,0 +1,78 @@ +import pytest +import numpy as np + +from firedrake import * +from firedrake.adjoint import * +from pyadjoint import get_working_tape +from gusto import * + + +@pytest.fixture(autouse=True) +def handle_taping(): + yield + tape = get_working_tape() + tape.clear_tape() + + +@pytest.fixture(autouse=True, scope="module") +def handle_annotation(): + from firedrake.adjoint import annotate_tape, continue_annotation + if not annotate_tape(): + continue_annotation() + yield + # Ensure annotation is paused when we finish. + annotate = annotate_tape() + if annotate: + pause_annotation() + + +@pytest.mark.parametrize("nu_is_control", [True, False]) +def test_diffusion_sensitivity(nu_is_control, tmpdir): + assert get_working_tape()._blocks == [] + n = 30 + mesh = PeriodicUnitSquareMesh(n, n) + output = OutputParameters(dirname=str(tmpdir)) + dt = 0.01 + domain = Domain(mesh, 10*dt, family="BDM", degree=1) + io = IO(domain, output) + + V = VectorFunctionSpace(mesh, "CG", 2) + domain.spaces.add_space("vecCG", V) + + R = FunctionSpace(mesh, "R", 0) + # We need to define nu as a function in order to have a control variable. + nu = Function(R, val=0.0001) + diffusion_params = DiffusionParameters(mesh, kappa=nu) + eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params) + + diffusion_scheme = BackwardEuler(domain) + diffusion_methods = [CGDiffusion(eqn, "f", diffusion_params)] + timestepper = Timestepper(eqn, diffusion_scheme, io, spatial_methods=diffusion_methods) + + x = SpatialCoordinate(mesh) + fexpr = as_vector((sin(2*pi*x[0]), cos(2*pi*x[1]))) + timestepper.fields("f").interpolate(fexpr) + + end = 0.1 + timestepper.run(0., end) + + u = timestepper.fields("f") + J = assemble(inner(u, u)*dx) + + if nu_is_control: + control = Control(nu) + h = Function(R, val=0.0001) # the direction of the perturbation + else: + control = Control(u) + # the direction of the perturbation + h = Function(V).interpolate(fexpr * np.random.rand()) + + # the functional as a pure function of nu + Jhat = ReducedFunctional(J, control) + + if nu_is_control: + assert np.allclose(J, Jhat(nu)) + assert taylor_test(Jhat, nu, h) > 1.95 + else: + assert np.allclose(J, Jhat(u)) + assert taylor_test(Jhat, u, h) > 1.95 diff --git a/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py b/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py new file mode 100644 index 000000000..17e527d6f --- /dev/null +++ b/integration-tests/adjoints/test_moist_thermal_williamson_5_sensitivity.py @@ -0,0 +1,218 @@ +""" +The moist thermal form of Test Case 5 (flow over a mountain) of Williamson et +al, 1992: +``A standard test set for numerical approximations to the shallow water +equations in spherical geometry'', JCP. + +The initial conditions are taken from Zerroukat & Allen, 2015: +``A moist Boussinesq shallow water equations set for testing atmospheric +models'', JCP. + +Three moist variables (vapour, cloud liquid and rain) are used. This set of +equations involves an active buoyancy field. + +The example here uses the icosahedral sphere mesh and degree 1 spaces. An +explicit RK4 timestepper is used. +""" +import pytest +import numpy as np + +from firedrake import ( + SpatialCoordinate, as_vector, pi, sqrt, min_value, exp, cos, sin, assemble, dx, inner, Function +) +from firedrake.adjoint import * +from pyadjoint import get_working_tape +from gusto import ( + Domain, IO, OutputParameters, Timestepper, RK4, DGUpwind, + ShallowWaterParameters, ThermalShallowWaterEquations, lonlatr_from_xyz, + InstantRain, SWSaturationAdjustment, WaterVapour, CloudWater, + Rain, GeneralIcosahedralSphereMesh +) + + +@pytest.fixture(autouse=True) +def handle_taping(): + yield + tape = get_working_tape() + tape.clear_tape() + + +@pytest.fixture(autouse=True, scope="module") +def handle_annotation(): + from firedrake.adjoint import annotate_tape, continue_annotation + if not annotate_tape(): + continue_annotation() + yield + # Ensure annotation is paused when we finish. + annotate = annotate_tape() + if annotate: + pause_annotation() + + +def test_moist_thermal_williamson_5_sensitivity( + tmpdir, ncells_per_edge=8, dt=600, tmax=50.*24.*60.*60., + dumpfreq=2880 +): + assert get_working_tape()._blocks == [] + # ------------------------------------------------------------------------ # + # Parameters for test case + # ------------------------------------------------------------------------ # + radius = 6371220. # planetary radius (m) + mean_depth = 5960 # reference depth (m) + g = 9.80616 # acceleration due to gravity (m/s^2) + u_max = 20. # max amplitude of the zonal wind (m/s) + epsilon = 1/300 # linear air expansion coeff (1/K) + theta_SP = -40*epsilon # value of theta at south pole (no units) + theta_EQ = 30*epsilon # value of theta at equator (no units) + theta_NP = -20*epsilon # value of theta at north pole (no units) + mu1 = 0.05 # scaling for theta with longitude (no units) + mu2 = 0.98 # proportion of qsat to make init qv (no units) + q0 = 135 # qsat scaling, gives init q_v of ~0.02, (kg/kg) + beta2 = 10*g # buoyancy-vaporisation factor (m/s^2) + nu = 20. # qsat factor in exponent (no units) + qprecip = 1e-4 # cloud to rain conversion threshold (kg/kg) + gamma_r = 1e-3 # rain-coalescence implicit factor + mountain_height = 2000. # height of mountain (m) + R0 = pi/9. # radius of mountain (rad) + lamda_c = -pi/2. # longitudinal centre of mountain (rad) + phi_c = pi/6. # latitudinal centre of mountain (rad) + + # ------------------------------------------------------------------------ # + # Our settings for this set up + # ------------------------------------------------------------------------ # + + element_order = 1 + u_eqn_type = 'vector_invariant_form' + + # ------------------------------------------------------------------------ # + # Set up model objects + # ------------------------------------------------------------------------ # + + # Domain + mesh = GeneralIcosahedralSphereMesh(radius, ncells_per_edge, degree=2) + domain = Domain(mesh, dt, "BDM", element_order) + x, y, z = SpatialCoordinate(mesh) + lamda, phi, _ = lonlatr_from_xyz(x, y, z) + + # Equation: coriolis + parameters = ShallowWaterParameters(mesh, H=mean_depth, g=g) + Omega = parameters.Omega + fexpr = 2*Omega*z/radius + + # Equation: topography + rsq = min_value(R0**2, (lamda - lamda_c)**2 + (phi - phi_c)**2) + r = sqrt(rsq) + tpexpr = mountain_height * (1 - r/R0) + + # Equation: moisture + tracers = [ + WaterVapour(space='DG'), CloudWater(space='DG'), Rain(space='DG') + ] + eqns = ThermalShallowWaterEquations( + domain, parameters, fexpr=fexpr, topog_expr=tpexpr, + active_tracers=tracers, u_transport_option=u_eqn_type + ) + + # I/O + output = OutputParameters( + dirname=str(tmpdir), dumplist_latlon=['D'], dumpfreq=dumpfreq, + dump_vtus=True, dump_nc=False, log_courant=False, + dumplist=['D', 'b', 'water_vapour', 'cloud_water'] + ) + io = IO(domain, output) + + # Physics ------------------------------------------------------------------ + # Saturation function -- first define simple expression + def q_sat(b, D): + return (q0/(g*D + g*tpexpr)) * exp(nu*(1 - b/g)) + + # Function to pass to physics (takes mixed function as argument) + def phys_sat_func(x_in): + D = x_in.sub(1) + b = x_in.sub(2) + return q_sat(b, D) + + # Feedback proportionality is dependent on D and b + def gamma_v(x_in): + D = x_in.sub(1) + b = x_in.sub(2) + return 1.0 / (1.0 + nu*beta2/g*q_sat(b, D)) + + SWSaturationAdjustment( + eqns, phys_sat_func, time_varying_saturation=True, + parameters=parameters, thermal_feedback=True, + beta2=beta2, gamma_v=gamma_v, time_varying_gamma_v=True + ) + + InstantRain( + eqns, qprecip, vapour_name="cloud_water", rain_name="rain", + gamma_r=gamma_r + ) + + transport_methods = [ + DGUpwind(eqns, field_name) for field_name in eqns.field_names + ] + + # Timestepper + stepper = Timestepper( + eqns, RK4(domain), io, spatial_methods=transport_methods + ) + + # ------------------------------------------------------------------------ # + # Initial conditions + # ------------------------------------------------------------------------ # + + u0 = stepper.fields("u") + D0 = stepper.fields("D") + b0 = stepper.fields("b") + v0 = stepper.fields("water_vapour") + c0 = stepper.fields("cloud_water") + r0 = stepper.fields("rain") + + uexpr = as_vector([-u_max*y/radius, u_max*x/radius, 0.0]) + + Dexpr = ( + mean_depth - tpexpr + - (radius * Omega * u_max + 0.5*u_max**2)*(z/radius)**2/g + ) + + # Expression for initial buoyancy - note the bracket around 1-mu + theta_expr = ( + 2/(pi**2) * ( + phi*(phi - pi/2)*theta_SP + - 2*(phi + pi/2) * (phi - pi/2)*(1 - mu1)*theta_EQ + + phi*(phi + pi/2)*theta_NP + ) + + mu1*theta_EQ*cos(phi)*sin(lamda) + ) + bexpr = g * (1 - theta_expr) + + # Expression for initial vapour depends on initial saturation + vexpr = mu2 * q_sat(bexpr, Dexpr) + + # Initialise (cloud and rain initially zero) + u0.project(uexpr) + D0.interpolate(Dexpr) + b0.interpolate(bexpr) + v0.interpolate(vexpr) + c0.assign(0.0) + r0.assign(0.0) + + # ----------------------------------------------------------------- # + # Run + # ----------------------------------------------------------------- # + stepper.run(t=0, tmax=dt) + + u_tf = stepper.fields('u') # Final velocity field + D_tf = stepper.fields('D') # Final depth field + + J = assemble(0.5*inner(u_tf, u_tf)*dx + 0.5*g*D_tf**2*dx) + + Jhat = ReducedFunctional(J, Control(D0)) + + assert np.allclose(Jhat(D0), J) + with stop_annotating(): + # Stop annotation to perform the Taylor test + h0 = Function(D0.function_space()) + h0.assign(D0 * np.random.rand()) + assert taylor_test(Jhat, D0, h0) > 1.95 diff --git a/integration-tests/adjoints/test_shallow_water_sensitivity.py b/integration-tests/adjoints/test_shallow_water_sensitivity.py new file mode 100644 index 000000000..22ffb9681 --- /dev/null +++ b/integration-tests/adjoints/test_shallow_water_sensitivity.py @@ -0,0 +1,100 @@ +import pytest +import numpy as np + +from firedrake import * +from firedrake.adjoint import * +from pyadjoint import get_working_tape +from gusto import * + + +@pytest.fixture(autouse=True) +def handle_taping(): + yield + tape = get_working_tape() + tape.clear_tape() + + +@pytest.fixture(autouse=True, scope="module") +def handle_annotation(): + from firedrake.adjoint import annotate_tape, continue_annotation + if not annotate_tape(): + continue_annotation() + yield + # Ensure annotation is paused when we finish. + annotate = annotate_tape() + if annotate: + pause_annotation() + + +def test_shallow_water(tmpdir): + assert get_working_tape()._blocks == [] + # setup shallow water parameters + R = 6371220. + H = 5960. + dt = 900. + + # Domain + mesh = IcosahedralSphereMesh(radius=R, refinement_level=3, degree=2) + x = SpatialCoordinate(mesh) + domain = Domain(mesh, dt, 'BDM', 1) + parameters = ShallowWaterParameters(mesh, H=H) + + # Equation + Omega = parameters.Omega + fexpr = 2*Omega*x[2]/R + lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2]) + R0 = pi/9. + R0sq = R0**2 + lamda_c = -pi/2. + lsq = (lamda - lamda_c)**2 + theta_c = pi/6. + thsq = (theta - theta_c)**2 + rsq = min_value(R0sq, lsq+thsq) + r = sqrt(rsq) + bexpr = 2000 * (1 - r/R0) + eqn = ShallowWaterEquations(domain, parameters, fexpr=fexpr, topog_expr=bexpr) + + # I/O + output = OutputParameters(dirname=str(tmpdir), log_courant=False) + io = IO(domain, output) + + # Transport schemes + transported_fields = [TrapeziumRule(domain, "u"), SSPRK3(domain, "D")] + transport_methods = [DGUpwind(eqn, "u"), DGUpwind(eqn, "D")] + + # Time stepper + stepper = SemiImplicitQuasiNewton( + eqn, io, transported_fields, transport_methods + ) + + u0 = stepper.fields('u') + D0 = stepper.fields('D') + u_max = 20. # Maximum amplitude of the zonal wind (m/s) + uexpr = as_vector([-u_max*x[1]/R, u_max*x[0]/R, 0.0]) + g = parameters.g + Rsq = R**2 + Dexpr = H - ((R * Omega * u_max + 0.5*u_max**2)*x[2]**2/Rsq)/g - bexpr + + u0.project(uexpr) + D0.interpolate(Dexpr) + + Dbar = Function(D0.function_space()).assign(H) + stepper.set_reference_profiles([('D', Dbar)]) + + stepper.run(0., 5*dt) + + u_tf = stepper.fields('u') # Final velocity field + D_tf = stepper.fields('D') # Final depth field + + J = assemble(0.5*inner(u_tf, u_tf)*dx + 0.5*g*D_tf**2*dx) + + control = [Control(D0), Control(u0)] # Control variables + J_hat = ReducedFunctional(J, control) + assert np.isclose(J_hat([D0, u0]), J, rtol=1e-10) + with stop_annotating(): + # Stop annotation to perform the Taylor test + h0 = Function(D0.function_space()) + h1 = Function(u0.function_space()) + h0.assign(D0 * np.random.rand()) + h1.assign(u0 * np.random.rand()) + assert taylor_test(J_hat, [D0, u0], [h0, h1]) > 1.95 diff --git a/integration-tests/balance/test_boussinesq_balance.py b/integration-tests/balance/test_boussinesq_balance.py index 3fae46620..21e55d9e4 100644 --- a/integration-tests/balance/test_boussinesq_balance.py +++ b/integration-tests/balance/test_boussinesq_balance.py @@ -31,7 +31,7 @@ def setup_balance(dirname): domain = Domain(mesh, dt, "CG", 1) # Equation - parameters = BoussinesqParameters() + parameters = BoussinesqParameters(mesh) eqns = BoussinesqEquations(domain, parameters) # I/O diff --git a/integration-tests/balance/test_compressible_balance.py b/integration-tests/balance/test_compressible_balance.py index a9726dd5a..29edcfa8d 100644 --- a/integration-tests/balance/test_compressible_balance.py +++ b/integration-tests/balance/test_compressible_balance.py @@ -31,7 +31,7 @@ def setup_balance(dirname): domain = Domain(mesh, dt, "CG", 1) # Equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqns = CompressibleEulerEquations(domain, parameters) # I/O diff --git a/integration-tests/balance/test_saturated_balance.py b/integration-tests/balance/test_saturated_balance.py index 2c89da089..42124bcb6 100644 --- a/integration-tests/balance/test_saturated_balance.py +++ b/integration-tests/balance/test_saturated_balance.py @@ -39,7 +39,7 @@ def setup_saturated(dirname, recovered): u_transport_option = "vector_advection_form" else: u_transport_option = "vector_invariant_form" - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqns = CompressibleEulerEquations( domain, parameters, u_transport_option=u_transport_option, active_tracers=tracers) diff --git a/integration-tests/balance/test_unsaturated_balance.py b/integration-tests/balance/test_unsaturated_balance.py index c6e3f33a9..d3a6463f8 100644 --- a/integration-tests/balance/test_unsaturated_balance.py +++ b/integration-tests/balance/test_unsaturated_balance.py @@ -42,7 +42,7 @@ def setup_unsaturated(dirname, recovered): u_transport_option = "vector_advection_form" else: u_transport_option = "vector_invariant_form" - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqns = CompressibleEulerEquations( domain, parameters, u_transport_option=u_transport_option, active_tracers=tracers) diff --git a/integration-tests/diffusion/test_diffusion.py b/integration-tests/diffusion/test_diffusion.py index fd570107d..149664d69 100644 --- a/integration-tests/diffusion/test_diffusion.py +++ b/integration-tests/diffusion/test_diffusion.py @@ -4,7 +4,8 @@ """ from gusto import * -from firedrake import (VectorFunctionSpace, Constant, as_vector, errornorm) +from firedrake import (VectorFunctionSpace, as_vector, errornorm, + TensorFunctionSpace) import pytest @@ -32,7 +33,7 @@ def test_scalar_diffusion(tmpdir, DG, tracer_setup): mu = 5. - diffusion_params = DiffusionParameters(kappa=kappa, mu=mu) + diffusion_params = DiffusionParameters(domain.mesh, kappa=kappa, mu=mu) eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params) diffusion_scheme = BackwardEuler(domain) diffusion_methods = [InteriorPenaltyDiffusion(eqn, "f", diffusion_params)] @@ -52,21 +53,23 @@ def test_vector_diffusion(tmpdir, DG, tracer_setup): f_init = setup.f_init tmax = setup.tmax tol = 3.e-2 - kappa = 1. f_end_expr = (1/(1+4*tmax))*f_init**(1/(1+4*tmax)) - kappa = Constant([[kappa, 0.], [0., kappa]]) + VR = TensorFunctionSpace(domain.mesh, "R", 0) + kappa = Function(VR, val=((1., 0.), (0., 1.))) + if DG: V = VectorFunctionSpace(domain.mesh, "DG", 1) else: V = domain.spaces("HDiv") + f_init = as_vector([f_init, 0.]) f_end_expr = as_vector([f_end_expr, 0.]) mu = 5. - diffusion_params = DiffusionParameters(kappa=kappa, mu=mu) + diffusion_params = DiffusionParameters(domain.mesh, kappa=kappa, mu=mu) eqn = DiffusionEquation(domain, V, "f", diffusion_parameters=diffusion_params) diffusion_scheme = BackwardEuler(domain) diffusion_methods = [InteriorPenaltyDiffusion(eqn, "f", diffusion_params)] diff --git a/integration-tests/equations/test_advection_diffusion.py b/integration-tests/equations/test_advection_diffusion.py index 778eb9f73..aac3ab2ee 100644 --- a/integration-tests/equations/test_advection_diffusion.py +++ b/integration-tests/equations/test_advection_diffusion.py @@ -22,7 +22,7 @@ def run_advection_diffusion(tmpdir): domain = Domain(mesh, dt, "CG", 1) # Equation - diffusion_params = DiffusionParameters(kappa=0.75, mu=5) + diffusion_params = DiffusionParameters(mesh, kappa=0.75, mu=5) V = domain.spaces("DG") Vu = VectorFunctionSpace(mesh, "CG", 1) diff --git a/integration-tests/equations/test_boussinesq.py b/integration-tests/equations/test_boussinesq.py index 2655d87d3..8de65aef6 100644 --- a/integration-tests/equations/test_boussinesq.py +++ b/integration-tests/equations/test_boussinesq.py @@ -31,7 +31,7 @@ def run_boussinesq(tmpdir, compressible): domain = Domain(mesh, dt, "CG", 1) # Equation - parameters = BoussinesqParameters() + parameters = BoussinesqParameters(mesh) eqn = BoussinesqEquations(domain, parameters, compressible=compressible) # I/O diff --git a/integration-tests/equations/test_dry_compressible.py b/integration-tests/equations/test_dry_compressible.py index 846d38bbf..f1d24c014 100644 --- a/integration-tests/equations/test_dry_compressible.py +++ b/integration-tests/equations/test_dry_compressible.py @@ -30,7 +30,7 @@ def run_dry_compressible(tmpdir): domain = Domain(mesh, dt, "CG", 1) # Equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters) # I/O @@ -98,7 +98,8 @@ def run_dry_compressible(tmpdir): checkpoint=True) check_mesh = pick_up_mesh(check_output, mesh_name) check_domain = Domain(check_mesh, dt, "CG", 1) - check_eqn = CompressibleEulerEquations(check_domain, parameters) + check_parameters = CompressibleParameters(check_mesh) + check_eqn = CompressibleEulerEquations(check_domain, check_parameters) check_io = IO(check_domain, check_output) check_stepper = SemiImplicitQuasiNewton(check_eqn, check_io, [], []) check_stepper.io.pick_up_from_checkpoint(check_stepper.fields) diff --git a/integration-tests/equations/test_linear_sw_wave.py b/integration-tests/equations/test_linear_sw_wave.py index 95eecdca4..e7cd59d74 100644 --- a/integration-tests/equations/test_linear_sw_wave.py +++ b/integration-tests/equations/test_linear_sw_wave.py @@ -28,7 +28,7 @@ def run_linear_sw_wave(tmpdir): domain = Domain(mesh, dt, 'BDM', 1) # Equation - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) fexpr = Constant(1) eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr) diff --git a/integration-tests/equations/test_moist_compressible.py b/integration-tests/equations/test_moist_compressible.py index f2909f75e..6b8d937f7 100644 --- a/integration-tests/equations/test_moist_compressible.py +++ b/integration-tests/equations/test_moist_compressible.py @@ -30,7 +30,7 @@ def run_moist_compressible(tmpdir): domain = Domain(mesh, dt, "CG", 1) # Equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) tracers = [WaterVapour(name='vapour_mixing_ratio'), CloudWater(name='cloud_liquid_mixing_ratio')] eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers) @@ -104,7 +104,9 @@ def run_moist_compressible(tmpdir): checkpoint=True) check_mesh = pick_up_mesh(check_output, mesh_name) check_domain = Domain(check_mesh, dt, "CG", 1) - check_eqn = CompressibleEulerEquations(check_domain, parameters, active_tracers=tracers) + check_parameters = CompressibleParameters(check_mesh) + check_eqn = CompressibleEulerEquations(check_domain, check_parameters, + active_tracers=tracers) check_io = IO(check_domain, output=check_output) check_stepper = SemiImplicitQuasiNewton(check_eqn, check_io, [], []) check_stepper.io.pick_up_from_checkpoint(check_stepper.fields) diff --git a/integration-tests/equations/test_sw_fplane.py b/integration-tests/equations/test_sw_fplane.py index 3d0acff5f..711adb7bb 100644 --- a/integration-tests/equations/test_sw_fplane.py +++ b/integration-tests/equations/test_sw_fplane.py @@ -23,7 +23,7 @@ def run_sw_fplane(tmpdir): # Equation H = 2 g = 50 - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) f0 = 10 fexpr = Constant(f0) eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr) diff --git a/integration-tests/equations/test_sw_linear_triangle.py b/integration-tests/equations/test_sw_linear_triangle.py index d533b3286..e34bb122f 100644 --- a/integration-tests/equations/test_sw_linear_triangle.py +++ b/integration-tests/equations/test_sw_linear_triangle.py @@ -30,7 +30,7 @@ def setup_sw(dirname): domain = Domain(mesh, dt, "BDM", degree=1) # Equation - parameters = ShallowWaterParameters(H=H) + parameters = ShallowWaterParameters(mesh, H=H) Omega = parameters.Omega fexpr = 2*Omega*x[2]/R eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr) diff --git a/integration-tests/equations/test_sw_triangle.py b/integration-tests/equations/test_sw_triangle.py index 455595c31..bc25b2e63 100644 --- a/integration-tests/equations/test_sw_triangle.py +++ b/integration-tests/equations/test_sw_triangle.py @@ -36,7 +36,7 @@ def setup_sw(dirname, dt, u_transport_option): x = SpatialCoordinate(mesh) # Equation - parameters = ShallowWaterParameters(H=H) + parameters = ShallowWaterParameters(mesh, H=H) Omega = parameters.Omega fexpr = 2*Omega*x[2]/R eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr, diff --git a/integration-tests/equations/test_thermal_sw.py b/integration-tests/equations/test_thermal_sw.py index f201e78c6..09b706fe9 100644 --- a/integration-tests/equations/test_thermal_sw.py +++ b/integration-tests/equations/test_thermal_sw.py @@ -37,7 +37,7 @@ def setup_sw(dirname, dt, u_transport_option): x = SpatialCoordinate(mesh) # Equation - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) Omega = parameters.Omega fexpr = 2*Omega*x[2]/R eqns = ThermalShallowWaterEquations(domain, parameters, fexpr=fexpr, diff --git a/integration-tests/model/test_checkpointing.py b/integration-tests/model/test_checkpointing.py index dc72fb1de..09bb5ff16 100644 --- a/integration-tests/model/test_checkpointing.py +++ b/integration-tests/model/test_checkpointing.py @@ -15,7 +15,7 @@ def set_up_model_objects(mesh, dt, output, stepper_type, ref_update_freq): domain = Domain(mesh, dt, "CG", 1) - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqns = CompressibleEulerEquations(domain, parameters) # Have two diagnostic fields that depend on initial values -- check if diff --git a/integration-tests/model/test_passive_tracer.py b/integration-tests/model/test_passive_tracer.py index 41bd5608e..2482722b0 100644 --- a/integration-tests/model/test_passive_tracer.py +++ b/integration-tests/model/test_passive_tracer.py @@ -26,7 +26,7 @@ def run_tracer(setup): x = SpatialCoordinate(mesh) H = 0.1 - parameters = ShallowWaterParameters(H=H) + parameters = ShallowWaterParameters(mesh, H=H) Omega = parameters.Omega g = parameters.g umax = setup.umax diff --git a/integration-tests/model/test_simultaneous_SIQN.py b/integration-tests/model/test_simultaneous_SIQN.py index ea3a12327..ddc471997 100644 --- a/integration-tests/model/test_simultaneous_SIQN.py +++ b/integration-tests/model/test_simultaneous_SIQN.py @@ -62,7 +62,7 @@ def run_simult_SIQN(tmpdir, order): density_name='rho')] # Equation - params = CompressibleParameters() + params = CompressibleParameters(mesh) eqns = CompressibleEulerEquations( domain, params, active_tracers=tracers, u_transport_option=u_eqn_type ) @@ -212,7 +212,7 @@ def run_simult_SIQN(tmpdir, order): # wind initially zero u0.project(as_vector( - [Constant(0.0, domain=mesh), Constant(0.0, domain=mesh)] + [Constant(0.0), Constant(0.0)] )) stepper.set_reference_profiles( @@ -238,7 +238,8 @@ def run_simult_SIQN(tmpdir, order): checkpoint=True) check_mesh = pick_up_mesh(check_output, mesh_name) check_domain = Domain(check_mesh, dt, "CG", order) - check_eqn = CompressibleEulerEquations(check_domain, params, active_tracers=tracers, u_transport_option=u_eqn_type) + check_params = CompressibleParameters(check_mesh) + check_eqn = CompressibleEulerEquations(check_domain, check_params, active_tracers=tracers, u_transport_option=u_eqn_type) check_io = IO(check_domain, check_output) check_stepper = SemiImplicitQuasiNewton(check_eqn, check_io, [], []) check_stepper.io.pick_up_from_checkpoint(check_stepper.fields) diff --git a/integration-tests/model/test_split_timestepper.py b/integration-tests/model/test_split_timestepper.py index 5f675cfd6..3f3eff94c 100644 --- a/integration-tests/model/test_split_timestepper.py +++ b/integration-tests/model/test_split_timestepper.py @@ -25,7 +25,7 @@ def run_split_timestepper_adv_diff_physics(tmpdir, timestepper): domain = Domain(mesh, dt, "CG", 1) # Equation - diffusion_params = DiffusionParameters(kappa=0.75, mu=5) + diffusion_params = DiffusionParameters(mesh, kappa=0.75, mu=5) V = domain.spaces("DG") Vu = VectorFunctionSpace(mesh, "CG", 1) diff --git a/integration-tests/physics/test_boundary_layer_mixing.py b/integration-tests/physics/test_boundary_layer_mixing.py index cc76a7371..e7203fd56 100644 --- a/integration-tests/physics/test_boundary_layer_mixing.py +++ b/integration-tests/physics/test_boundary_layer_mixing.py @@ -33,7 +33,7 @@ def run_boundary_layer_mixing(dirname, field_name, recovered, semi_implicit): _, z = SpatialCoordinate(mesh) # Set up equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters) # I/O @@ -43,7 +43,7 @@ def run_boundary_layer_mixing(dirname, field_name, recovered, semi_implicit): io = IO(domain, output) # Physics scheme - surf_params = BoundaryLayerParameters() + surf_params = BoundaryLayerParameters(mesh) physics_parametrisation = BoundaryLayerMixing(eqn, field_name, surf_params) if recovered: diff --git a/integration-tests/physics/test_held_suarez_friction.py b/integration-tests/physics/test_held_suarez_friction.py index dfe7e34f9..b1b44b9b0 100644 --- a/integration-tests/physics/test_held_suarez_friction.py +++ b/integration-tests/physics/test_held_suarez_friction.py @@ -29,7 +29,7 @@ def run_apply_rayleigh_friction(dirname): domain = Domain(mesh, dt, "CG", 0) # Set up equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters) # I/O diff --git a/integration-tests/physics/test_held_suarez_relaxation.py b/integration-tests/physics/test_held_suarez_relaxation.py index 615ba8d29..17f550dc6 100644 --- a/integration-tests/physics/test_held_suarez_relaxation.py +++ b/integration-tests/physics/test_held_suarez_relaxation.py @@ -30,7 +30,7 @@ def run_held_suarez_relaxation(dirname, temp): domain = Domain(mesh, dt, "CG", 0) # Set up equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters) # I/O diff --git a/integration-tests/physics/test_instant_rain.py b/integration-tests/physics/test_instant_rain.py index 3cb42d839..1f7b31ce3 100644 --- a/integration-tests/physics/test_instant_rain.py +++ b/integration-tests/physics/test_instant_rain.py @@ -39,7 +39,7 @@ def run_instant_rain(dirname, physics_coupling): rain = Rain(name="rain", space="DG", transport_eqn=TransportEquationType.no_transport) - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr, active_tracers=[vapour, rain]) diff --git a/integration-tests/physics/test_saturation_adjustment.py b/integration-tests/physics/test_saturation_adjustment.py index b39ee0345..c36c09e37 100644 --- a/integration-tests/physics/test_saturation_adjustment.py +++ b/integration-tests/physics/test_saturation_adjustment.py @@ -38,7 +38,7 @@ def run_cond_evap(dirname, process, physics_coupling): # Set up equation tracers = [WaterVapour(), CloudWater()] - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers) # I/O diff --git a/integration-tests/physics/test_static_adjustment.py b/integration-tests/physics/test_static_adjustment.py index 0803d9b2c..560fc6ebf 100644 --- a/integration-tests/physics/test_static_adjustment.py +++ b/integration-tests/physics/test_static_adjustment.py @@ -36,7 +36,7 @@ def run_static_adjustment(dirname, theta_variable): # Set up equation tracers = [WaterVapour()] - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers) # I/O diff --git a/integration-tests/physics/test_suppress_vertical_wind.py b/integration-tests/physics/test_suppress_vertical_wind.py index fee46b84a..6f70d7fe0 100644 --- a/integration-tests/physics/test_suppress_vertical_wind.py +++ b/integration-tests/physics/test_suppress_vertical_wind.py @@ -35,7 +35,7 @@ def run_suppress_vertical_wind(dirname, physics_coupling): # Set up equation tracers = [WaterVapour()] - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers) # I/O diff --git a/integration-tests/physics/test_surface_fluxes.py b/integration-tests/physics/test_surface_fluxes.py index 14990941c..fee7dabee 100644 --- a/integration-tests/physics/test_surface_fluxes.py +++ b/integration-tests/physics/test_surface_fluxes.py @@ -37,7 +37,7 @@ def run_surface_fluxes(dirname, moist, implicit_formulation, physics_coupling): # Set up equation tracers = [WaterVapour()] if moist else None vapour_name = 'water_vapour' if moist else None - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters, active_tracers=tracers) # I/O @@ -47,7 +47,7 @@ def run_surface_fluxes(dirname, moist, implicit_formulation, physics_coupling): io = IO(domain, output) # Physics scheme - surf_params = BoundaryLayerParameters() + surf_params = BoundaryLayerParameters(mesh) T_surf = Constant(300.0) if physics_coupling == "split": diff --git a/integration-tests/physics/test_sw_saturation_adjustment.py b/integration-tests/physics/test_sw_saturation_adjustment.py index 5b430a45c..ecdfdf2f5 100644 --- a/integration-tests/physics/test_sw_saturation_adjustment.py +++ b/integration-tests/physics/test_sw_saturation_adjustment.py @@ -43,7 +43,7 @@ def run_sw_cond_evap(dirname, process, physics_coupling): sat = 100 # Equation - parameters = ShallowWaterParameters(H=H) + parameters = ShallowWaterParameters(mesh, H=H) Omega = parameters.Omega fexpr = 2*Omega*x[2]/R diff --git a/integration-tests/physics/test_wind_drag.py b/integration-tests/physics/test_wind_drag.py index 4214ac1e5..8720ec69a 100644 --- a/integration-tests/physics/test_wind_drag.py +++ b/integration-tests/physics/test_wind_drag.py @@ -33,7 +33,7 @@ def run_wind_drag(dirname, implicit_formulation, physics_coupling): _, z = SpatialCoordinate(mesh) # Set up equation - parameters = CompressibleParameters() + parameters = CompressibleParameters(mesh) eqn = CompressibleEulerEquations(domain, parameters) # I/O @@ -43,7 +43,7 @@ def run_wind_drag(dirname, implicit_formulation, physics_coupling): io = IO(domain, output) # Physics scheme - surf_params = BoundaryLayerParameters() + surf_params = BoundaryLayerParameters(mesh) physics_parametrisation = WindDrag(eqn, implicit_formulation, surf_params) if physics_coupling == "split": diff --git a/integration-tests/rexi/test_linear_compressible_boussinesq.py b/integration-tests/rexi/test_linear_compressible_boussinesq.py index 89a3c85d7..abcc9d8f6 100644 --- a/integration-tests/rexi/test_linear_compressible_boussinesq.py +++ b/integration-tests/rexi/test_linear_compressible_boussinesq.py @@ -31,7 +31,7 @@ def run_rexi_linear_boussinesq(tmpdir): domain = Domain(mesh, tmax, 'CG', 1) # Equation - parameters = BoussinesqParameters(cs=300) + parameters = BoussinesqParameters(mesh, cs=300) eqns = LinearBoussinesqEquations(domain, parameters) # ---------------------------------------------------------------------- # diff --git a/integration-tests/rexi/test_linear_sw.py b/integration-tests/rexi/test_linear_sw.py index 77d5eafe2..57214eb7a 100644 --- a/integration-tests/rexi/test_linear_sw.py +++ b/integration-tests/rexi/test_linear_sw.py @@ -45,7 +45,7 @@ def run_rexi_sw(tmpdir, ensemble=None): domain = Domain(mesh, tmax, 'BDM', 1) # Equation - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) fexpr = Constant(f) eqns = LinearShallowWaterEquations(domain, parameters, fexpr=fexpr) diff --git a/integration-tests/transport/test_zero_limiter.py b/integration-tests/transport/test_zero_limiter.py index 0604dd3d3..5ece6ca1d 100644 --- a/integration-tests/transport/test_zero_limiter.py +++ b/integration-tests/transport/test_zero_limiter.py @@ -43,7 +43,7 @@ def setup_zero_limiter(dirname, limiter=False, rain=False): x = SpatialCoordinate(mesh) # Equations - parameters = ShallowWaterParameters(H=H, g=g) + parameters = ShallowWaterParameters(mesh, H=H, g=g) Omega = parameters.Omega fexpr = 2*Omega*x[2]/R diff --git a/unit-tests/diagnostic_tests/test_cumulative_sum.py b/unit-tests/diagnostic_tests/test_cumulative_sum.py index 295e522ca..918447fb9 100644 --- a/unit-tests/diagnostic_tests/test_cumulative_sum.py +++ b/unit-tests/diagnostic_tests/test_cumulative_sum.py @@ -18,7 +18,7 @@ def test_dewpoint(): mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers) domain = Domain(mesh, 0.1, 'CG', 1) - params = CompressibleParameters() + params = CompressibleParameters(mesh) active_tracers = [WaterVapour()] eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers) prog_fields = TimeLevelFields(eqn) diff --git a/unit-tests/diagnostic_tests/test_dewpoint_temperature.py b/unit-tests/diagnostic_tests/test_dewpoint_temperature.py index eaf1fbe55..1b4fb2980 100644 --- a/unit-tests/diagnostic_tests/test_dewpoint_temperature.py +++ b/unit-tests/diagnostic_tests/test_dewpoint_temperature.py @@ -18,7 +18,7 @@ def test_dewpoint(): mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers) domain = Domain(mesh, 0.1, 'CG', 1) - params = CompressibleParameters() + params = CompressibleParameters(mesh) active_tracers = [WaterVapour()] eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers) prog_fields = TimeLevelFields(eqn) diff --git a/unit-tests/diagnostic_tests/test_wetbulb_temperature.py b/unit-tests/diagnostic_tests/test_wetbulb_temperature.py index 5909fe710..1410a7b30 100644 --- a/unit-tests/diagnostic_tests/test_wetbulb_temperature.py +++ b/unit-tests/diagnostic_tests/test_wetbulb_temperature.py @@ -18,7 +18,7 @@ def test_wetbulb_temperature(): mesh = ExtrudedMesh(m, layers=nlayers, layer_height=H/nlayers) domain = Domain(mesh, 0.1, 'CG', 1) - params = CompressibleParameters() + params = CompressibleParameters(mesh) active_tracers = [WaterVapour()] eqn = CompressibleEulerEquations(domain, params, active_tracers=active_tracers) prog_fields = TimeLevelFields(eqn)