Restore check for explicit physics schemes and a few other tidyings

gusto/coordinates.py

@@ -7,6 +7,7 @@
 from gusto.logging import logger
 from firedrake import SpatialCoordinate, Function
 import numpy as np
+import pandas as pd
 
 
 class Coordinates(object):
@@ -170,8 +171,6 @@ def get_column_data(self, field, domain):
                 ordered column data.
         """
 
-        import pandas as pd
-
         space_name = field.function_space().name
         if space_name not in self.chi_coords.keys():
             self.register_space(domain, space_name)
@@ -189,13 +188,13 @@ def get_column_data(self, field, domain):
         # Work out digits to round to, based on number of points and range of coords
         num_points = np.size(coords_X)
         data_range = np.max(coords_X) - np.min(coords_X)
-        if data_range > 1e-16:
+        if data_range > np.finfo(type(data_range)).tiny:
             digits = int(np.floor(-np.log10(data_range / num_points)) + 3)
             coords_X = coords_X.round(digits)
 
         if data_is_3d:
             data_range = np.max(coords_Y) - np.min(coords_Y)
-            if data_range > 1e-16:
+            if data_range > np.finfo(type(data_range)).tiny:
                 # Only round if there is already some range
                 digits = int(np.floor(-np.log10(data_range / num_points)) + 3)
                 coords_Y = coords_Y.round(digits)

gusto/diagnostics.py

@@ -222,7 +222,7 @@ def setup(self, domain, state_fields, space=None):
     def compute(self):
         """Compute the diagnostic field from the current state."""
 
-        logger.info(f'Computing diagnostic {self.name} with {self.method} method')
+        logger.debug(f'Computing diagnostic {self.name} with {self.method} method')
 
         if self.method == 'interpolate':
             self.evaluator.interpolate()

gusto/domain.py

@@ -4,7 +4,6 @@
 model's time interval.
 """
 
-from gusto.logging import logger
 from gusto.coordinates import Coordinates
 from gusto.function_spaces import Spaces, check_degree_args
 from firedrake import (Constant, SpatialCoordinate, sqrt, CellNormal, cross,
@@ -176,24 +175,15 @@ def set_height_above_surface(self):
         height_above_surface = Function(CG1)
 
         # Turn height into columnwise data
-        try:
-            columnwise_height, index_data = self.coords.get_column_data(CG1_height, self)
-
-            # Find minimum height in each column
-            surface_height_1d = np.min(columnwise_height, axis=1)
-            height_above_surface_data = columnwise_height - surface_height_1d[:, None]
-
-            self.coords.set_field_from_column_data(height_above_surface,
-                                                   height_above_surface_data,
-                                                   index_data)
-        except ModuleNotFoundError:
-            # If no pandas, then don't take orography into account
-            logger.warning(
-                'Pandas not found, so using height as height above surface. '
-                + 'This will be incorrect in the presence of orography')
-
-            if self.on_sphere:
-                height_above_surface.assign(CG1_height - self.radius)
+        columnwise_height, index_data = self.coords.get_column_data(CG1_height, self)
+
+        # Find minimum height in each column
+        surface_height_1d = np.min(columnwise_height, axis=1)
+        height_above_surface_data = columnwise_height - surface_height_1d[:, None]
+
+        self.coords.set_field_from_column_data(height_above_surface,
+                                               height_above_surface_data,
+                                               index_data)
 
         self.height_above_surface = height_above_surface
 

gusto/logging.py

@@ -186,6 +186,7 @@ def update_logfile_location(new_path, comm):
         fh.close()
         logger.removeHandler(fh)
 
+        os.makedirs(new_path, exist_ok=True)
         if parallel_log in ["FILE", "BOTH"]:
             # If all ranks are logging wait here in case a directory is being created
             comm.Barrier()

gusto/physics.py

@@ -200,6 +200,7 @@ def __init__(self, equation, vapour_name='water_vapour',
         """
 
         label_name = 'saturation_adjustment'
+        self.explicit_only = True
         super().__init__(equation, label_name, parameters=parameters)
 
         # TODO: make a check on the variable type of the active tracers
@@ -262,7 +263,7 @@ def __init__(self, equation, vapour_name='water_vapour',
                 liquid_water += self.X.subfunctions[liq_idx]
 
         # define some parameters as attributes
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
         R_d = parameters.R_d
         cp = parameters.cp
         cv = parameters.cv
@@ -439,6 +440,7 @@ def __init__(self, equation, rain_name, domain, transport_method,
             terminal_velocity = Constant(5)  # in m/s
             v.project(-terminal_velocity*domain.k)
         elif moments == AdvectedMoments.M3:
+            self.explicit_only = True
             # this advects the third moment M3 of the raindrop
             # distribution, which corresponds to the mean mass
             rho_idx = equation.field_names.index('rho')
@@ -524,6 +526,7 @@ def __init__(self, equation, cloud_name='cloud_water', rain_name='rain',
                 process in the parametrisation. Defaults to True.
         """
 
+        self.explicit_only = True
         label_name = "coalescence"
         if accretion:
             label_name += "_accretion"
@@ -549,7 +552,7 @@ def __init__(self, equation, cloud_name='cloud_water', rain_name='rain',
         self.source = Function(Vt)
 
         # define some parameters as attributes
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
         # TODO: should these parameters be hard-coded or configurable?
         k_1 = Constant(0.001)  # accretion rate in 1/s
         k_2 = Constant(2.2)  # accumulation rate in 1/s
@@ -635,6 +638,7 @@ def __init__(self, equation, rain_name='rain', vapour_name='water_vapour',
                 CompressibleEulerEquations.
         """
 
+        self.explicit_only = True
         label_name = 'evaporation_of_rain'
         super().__init__(equation, label_name, parameters=None)
 
@@ -694,7 +698,7 @@ def __init__(self, equation, rain_name='rain', vapour_name='water_vapour',
                 liquid_water += self.X.subfunctions[liq_idx]
 
         # define some parameters as attributes
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
         R_d = parameters.R_d
         cp = parameters.cp
         cv = parameters.cv
@@ -819,6 +823,7 @@ def __init__(self, equation, saturation_curve,
                 parameters are obtained from the equation.
         """
 
+        self.explicit_only = True
         label_name = 'instant_rain'
         super().__init__(equation, label_name, parameters=parameters)
 
@@ -986,6 +991,7 @@ def __init__(self, equation, saturation_curve, L_v=None,
                 parameters are obtained from the equation.
         """
 
+        self.explicit_only = True
         label_name = 'saturation_adjustment'
         super().__init__(equation, label_name, parameters=parameters)
 
@@ -1144,8 +1150,8 @@ def __init__(self, equation, T_surface_expr, vapour_name=None,
             implicit_formulation (bool, optional): whether the scheme is already
                 put into a Backwards Euler formulation (which allows this scheme
                 to actually be used with a Forwards Euler or other explicit time
-                discretisation. Otherwise, this is formulated more generally and
-                can be used with any time stepper. Defaults to False.
+                discretisation). Otherwise, this is formulated more generally
+                and can be used with any time stepper. Defaults to False.
             parameters (:class:`BoundaryLayerParameters`): configuration object
                 giving the parameters for boundary and surface level schemes.
                 Defaults to None, in which case default values are used.
@@ -1169,7 +1175,7 @@ def __init__(self, equation, T_surface_expr, vapour_name=None,
 
         self.implicit_formulation = implicit_formulation
         self.X = Function(equation.X.function_space())
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
 
         # -------------------------------------------------------------------- #
         # Extract prognostic variables
@@ -1324,8 +1330,8 @@ def __init__(self, equation, implicit_formulation=False, parameters=None):
             implicit_formulation (bool, optional): whether the scheme is already
                 put into a Backwards Euler formulation (which allows this scheme
                 to actually be used with a Forwards Euler or other explicit time
-                discretisation. Otherwise, this is formulated more generally and
-                can be used with any time stepper. Defaults to False.
+                discretisation). Otherwise, this is formulated more generally
+                and can be used with any time stepper. Defaults to False.
             parameters (:class:`BoundaryLayerParameters`): configuration object
                 giving the parameters for boundary and surface level schemes.
                 Defaults to None, in which case default values are used.
@@ -1346,7 +1352,7 @@ def __init__(self, equation, implicit_formulation=False, parameters=None):
         k = equation.domain.k
         self.implicit_formulation = implicit_formulation
         self.X = Function(equation.X.function_space())
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
 
         # -------------------------------------------------------------------- #
         # Extract prognostic variables
@@ -1441,6 +1447,7 @@ def __init__(self, equation, theta_variable='theta_vd'):
                 to "theta_vd".
         """
 
+        self.explicit_only = True
         from functools import partial
 
         # -------------------------------------------------------------------- #
@@ -1457,7 +1464,7 @@ def __init__(self, equation, theta_variable='theta_vd'):
         super().__init__(equation, label_name, parameters=equation.parameters)
 
         self.X = Function(equation.X.function_space())
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
 
         # -------------------------------------------------------------------- #
         # Extract prognostic variables
@@ -1540,6 +1547,8 @@ def __init__(self, equation, spin_up_period):
                 period.
         """
 
+        self.explicit_only = True
+
         # -------------------------------------------------------------------- #
         # Check arguments and generic initialisation
         # -------------------------------------------------------------------- #
@@ -1554,7 +1563,7 @@ def __init__(self, equation, spin_up_period):
         super().__init__(equation, label_name, parameters=equation.parameters)
 
         self.X = Function(equation.X.function_space())
-        self.dt = Constant(1.0)
+        self.dt = Constant(0.0)
         self.t = domain.t
         self.spin_up_period = Constant(spin_up_period)
         self.strength = Constant(1.0)

gusto/timeloop.py

@@ -347,9 +347,11 @@ def __init__(self, equation, scheme, io, spatial_methods=None,
             self.physics_schemes = []
 
         for parametrisation, phys_scheme in self.physics_schemes:
-            # check that the supplied schemes for physics are explicit
-            # assert isinstance(phys_scheme, ExplicitTimeDiscretisation), \
-            #     "Only explicit time discretisations can be used for physics"
+            # check that the supplied schemes for physics are valid
+            if hasattr(parametrisation, "explicit_only") and parametrisation.explicit_only:
+                assert isinstance(phys_scheme, ExplicitTimeDiscretisation), \
+                    ("Only explicit time discretisations can be used with "
+                     + f"physics scheme {parametrisation.label.label}")
             apply_bcs = False
             phys_scheme.setup(equation, apply_bcs, parametrisation.label)
 
@@ -460,10 +462,12 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods,
                                     + self.fast_physics_schemes
                                     + self.final_physics_schemes)
 
-        for _, scheme in self.all_physics_schemes:
+        for parametrisation, scheme in self.all_physics_schemes:
             assert scheme.nlevels == 1, "multilevel schemes not supported as part of this timestepping loop"
-            # assert isinstance(scheme, ExplicitTimeDiscretisation), \
-            #     "Only explicit time discretisations can be used for physics"
+            if hasattr(parametrisation, "explicit_only") and parametrisation.explicit_only:
+                assert isinstance(scheme, ExplicitTimeDiscretisation), \
+                    ("Only explicit time discretisations can be used with "
+                     + f"physics scheme {parametrisation.label.label}")
 
         self.active_transport = []
         for scheme in transport_schemes: