From afa859c527f5747f5b9135ae0b501ced7af636b6 Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Sun, 22 Oct 2023 20:13:04 +0100 Subject: [PATCH] Restore check for explicit physics schemes and a few other tidyings --- gusto/coordinates.py | 7 +++---- gusto/diagnostics.py | 2 +- gusto/domain.py | 28 +++++++++------------------- gusto/logging.py | 1 + gusto/physics.py | 31 ++++++++++++++++++++----------- gusto/timeloop.py | 16 ++++++++++------ 6 files changed, 44 insertions(+), 41 deletions(-) diff --git a/gusto/coordinates.py b/gusto/coordinates.py index f8cf94f2d..571effbe8 100644 --- a/gusto/coordinates.py +++ b/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) diff --git a/gusto/diagnostics.py b/gusto/diagnostics.py index cfcde1a7d..024a42fa9 100644 --- a/gusto/diagnostics.py +++ b/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() diff --git a/gusto/domain.py b/gusto/domain.py index 4300920b3..8f22083c6 100644 --- a/gusto/domain.py +++ b/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 diff --git a/gusto/logging.py b/gusto/logging.py index c6eef1e5e..dc116275b 100644 --- a/gusto/logging.py +++ b/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() diff --git a/gusto/physics.py b/gusto/physics.py index 658163005..76c8ee4c3 100644 --- a/gusto/physics.py +++ b/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) diff --git a/gusto/timeloop.py b/gusto/timeloop.py index 7252d563c..c020c8bb8 100644 --- a/gusto/timeloop.py +++ b/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: