From b9bbdb1ec1b525957e424944e22c60420f386e69 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 12:04:28 -0600 Subject: [PATCH 01/17] Add verify_properties to step and ocean_model_step --- polaris/ocean/conservation.py | 117 ++++++++++++++++++++++++ polaris/ocean/model/ocean_model_step.py | 64 +++++++++++++ polaris/ocean/ocean.cfg | 4 + polaris/run/serial.py | 14 +++ polaris/step.py | 15 ++- 5 files changed, 213 insertions(+), 1 deletion(-) create mode 100644 polaris/ocean/conservation.py diff --git a/polaris/ocean/conservation.py b/polaris/ocean/conservation.py new file mode 100644 index 0000000000..f48114d94b --- /dev/null +++ b/polaris/ocean/conservation.py @@ -0,0 +1,117 @@ +from mpas_tools.cime.constants import constants + +# Should match config_density0 +rho_sw = 1026.0 + +cp_sw = constants['SHR_CONST_CPSW'] + + +def compute_total_mass(ds_mesh, ds): + ds = reduce_dataset_time_dim(ds) + area_cell = ds_mesh.areaCell + layer_thickness = ds.layerThickness + total_mass = rho_sw * ( + area_cell * layer_thickness.sum(dim='nVertLevels') + ).sum(dim='nCells') + return total_mass + + +def compute_total_mass_nonboussinesq(ds_mesh, ds): + ds = reduce_dataset_time_dim(ds) + area_cell = ds_mesh.areaCell + layer_thickness = ds.layerThickness + rho = ds.density + total_mass = rho * ( + area_cell * (layer_thickness * rho).sum(dim='nVertLevels') + ).sum(dim='nCells') + return total_mass + + +def compute_total_energy(ds_mesh, ds): + ds = reduce_dataset_time_dim(ds) + area_cell = ds_mesh.areaCell + layer_thickness = ds.layerThickness + temperature = ds.temperature + total_energy = ( + rho_sw + * cp_sw + * ( + area_cell * (layer_thickness * temperature).sum(dim='nVertLevels') + ).sum(dim='nCells') + ) + return total_energy + + +def compute_total_salt(ds_mesh, ds): + ds = reduce_dataset_time_dim(ds) + area_cell = ds_mesh.areaCell + layer_thickness = ds.layerThickness + salinity = ds.salinity + total_energy = ( + area_cell * (layer_thickness * salinity).sum(dim='nVertLevels') + ).sum(dim='nCells') + return total_energy + + +def reduce_dataset_time_dim(ds): + for time_dim in ['time', 'Time']: + if time_dim in ds.dims: + if ds.sizes[time_dim] > 1: + raise ValueError( + 'compute_total_energy is designed to work on a dataset ' + 'with one time slice' + ) + else: + ds = ds.isel(time_dim=0) + return ds + + +# def compute_net_mass_flux(ds): +# netMassFlux = & +# + accumulatedRainFlux & +# + accumulatedSnowFlux & +# + accumulatedEvaporationFlux & +# + accumulatedSeaIceFlux & +# + accumulatedRiverRunoffFlux & +# + accumulatedIceRunoffFlux & +# + accumulatedIcebergFlux & +# + accumulatedFrazilFlux & +# + accumulatedLandIceFlux +# note, accumulatedLandIceFrazilFlux not added because already in +# accumulatedFrazilFlux +# if (config_subglacial_runoff_mode == 'data'): +# netMassFlux = netMassFlux + accumulatedSubglacialRunoffFlux +# return net_mass_flux +# +# def compute_net_salt_flux(ds): +# netSaltFlux = accumulatedSeaIceSalinityFlux & +# + accumulatedFrazilSalinityFlux +# +# if (config_subglacial_runoff_mode == 'data'): +# netSaltFlux = netSaltFlux + accumulatedSubglacialRunoffSalinityFlux +# return net_salt_flux +# +# def compute_net_energy_flux(ds): +# netEnergyFlux = & +# accumulatedSeaIceHeatFlux & +# + accumulatedShortWaveHeatFlux & +# + accumulatedLongWaveHeatFluxUp & +# + accumulatedLongWaveHeatFluxDown & +# + accumulatedLatentHeatFlux & +# + accumulatedSensibleHeatFlux & +# + accumulatedMeltingSnowHeatFlux & +# + accumulatedMeltingIceRunoffHeatFlux & +# + accumulatedIcebergHeatFlux & +# + accumulatedFrazilHeatFlux & +# + accumulatedLandIceHeatFlux & +# + accumulatedRainTemperatureFlux *rho_sw*cp_sw & +# + accumulatedEvapTemperatureFlux *rho_sw*cp_sw & +# + accumulatedSeaIceTemperatureFlux *rho_sw*cp_sw & +# + accumulatedRiverRunoffTemperatureFlux *rho_sw*cp_sw & +# + accumulatedIcebergTemperatureFlux*rho_sw*cp_sw +# note, accumulatedLandIceFrazilHeatFlux not added because already in +# accumulatedFrazilHeatFlux +# if (config_subglacial_runoff_mode == 'data'): +# netEnergyFlux += accumulatedSubglacialRunoffTemperatureFlux * \ +# rho_sw*cp_sw +# return net_energy_flux diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index ac777a9e1f..eb8f39c024 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -1,10 +1,12 @@ import importlib.resources as imp_res +import os from types import ModuleType from typing import Any, Dict, List, Optional, Tuple, Union from ruamel.yaml import YAML from polaris.model_step import ModelStep +from polaris.ocean.conservation import compute_total_mass from polaris.tasks.ocean import Ocean OptionValue = Union[str, int, float, bool] @@ -362,6 +364,68 @@ def update_namelist_eos(self) -> None: options=replacements, config_model='ocean' ) + def verify_properties(self): + checked = False + success = True + if self.work_dir is None: + raise ValueError( + 'The work directory must be set before the step ' + 'output properties can be verified.' + ) + failed_properties = [] + for filename, properties in self.properties_to_verify.items(): + filename = str(filename) + mesh_filename = os.path.join(self.work_dir, 'mesh.nc') + this_filename = os.path.join(self.work_dir, filename) + ds_mesh = self.component.open_model_dataset(mesh_filename) + ds = self.component.open_model_dataset(this_filename) + for output_property in properties: + if output_property == 'mass conservation': + tol = self.config.getfloat( + 'ocean', 'mass_conservation_tolerance' + ) + init_mass = compute_total_mass(ds_mesh, ds.isel(Time=0)) + final_mass = compute_total_mass(ds_mesh, ds.isel(Time=-1)) + # mass_flux = compute_total_mass(ds.isel(Time=-1)) + result = abs(init_mass - final_mass) < tol + # absMassError = netMassFlux * dtAvg - massChange + # relMassError = absoluteMassError / (finalmass + 1.) + # absSaltError = netSaltFlux * dtAvg - saltChange + # relSaltError = absoluteSaltError / (finalSalt - 1.) + # absEnergyError = netEnergyFlux * dtAvg - energyChange + # relEnergyError = absoluteEnergyError / (finalEnergy - 1.) + else: + raise ValueError( + 'Could not find method to execute property check ' + f'{output_property}' + ) + + success = success and result + checked = True + if not result: + failed_properties.extend(output_property) + + if checked and success: + log_filename = os.path.join( + self.work_dir, 'property_check_passed.log' + ) + with open(log_filename, 'w') as result_log_file: + result_log_file.write( + f'Output file {filename} passed property checks.\n' + ) + elif checked and not success: + log_filename = os.path.join( + self.work_dir, 'property_check_failed.log' + ) + failed_properties_str = '\n '.join(failed_properties) + with open(log_filename, 'w') as result_log_file: + result_log_file.write( + f'Property checks on {filename} failed for:\n ' + f'{failed_properties_str}\n' + ) + + return checked, success + def _update_ntasks(self) -> None: """ Update ``ntasks`` and ``min_tasks`` for the step based on the estimated diff --git a/polaris/ocean/ocean.cfg b/polaris/ocean/ocean.cfg index d57eb21033..c7c9357031 100644 --- a/polaris/ocean/ocean.cfg +++ b/polaris/ocean/ocean.cfg @@ -3,6 +3,10 @@ # Options related the ocean component [ocean] + +# Tolerance for mass conservation +mass_conservation_tolerance = 1e-8 + # Which model, MPAS-Ocean or Omega, is used, or "detect" to auto-detect model = detect diff --git a/polaris/run/serial.py b/polaris/run/serial.py index dc9e4d32bd..f7eab3a637 100644 --- a/polaris/run/serial.py +++ b/polaris/run/serial.py @@ -543,6 +543,20 @@ def _run_task(task, available_resources): step_time = time.time() - step_start step_time_str = str(timedelta(seconds=round(step_time))) + compared, status = step.verify_properties() + if compared: + if status: + property_str = pass_str + else: + property_str = fail_str + _print_to_stdout( + task, f' property verification: {property_str}' + ) + if property_passed is None: + property_passed = status + elif not status: + property_passed = False + compared, status = step.validate_baselines() if compared: if status: diff --git a/polaris/step.py b/polaris/step.py index 361099bc6e..2006abb8c4 100644 --- a/polaris/step.py +++ b/polaris/step.py @@ -262,6 +262,7 @@ def __init__( # may be set during setup if there is a baseline for comparison self.baseline_dir = None self.validate_vars = dict() + self.properties_to_verify = dict() self.setup_complete = False # these will be set before running the step, dummy placeholders for now @@ -477,7 +478,9 @@ def add_input_file( ) ) - def add_output_file(self, filename, validate_vars=None): + def add_output_file( + self, filename, validate_vars=None, verify_properties=None + ): """ Add the output file that must be produced by this step and may be made available as an input to steps, perhaps in other tasks. This file @@ -499,6 +502,8 @@ def add_output_file(self, filename, validate_vars=None): self.outputs.append(filename) if validate_vars is not None: self.validate_vars[filename] = validate_vars + if verify_properties is not None: + self.properties_to_verify[filename] = verify_properties def add_dependency(self, step, name=None): """ @@ -534,6 +539,14 @@ def add_dependency(self, step, name=None): target = f'{step.path}/step_after_run.pickle' self.add_input_file(filename=filename, work_dir_target=target) + def verify_properties(self): + """ + This method should be overridden to verify properties of step outputs + """ + checked = False + success = True + return checked, success + def validate_baselines(self): """ Compare variables between output files in this step and in the same From eb95bc7f923b9b03dbdacd15bff3b4a1c102cfb7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 12:06:45 -0600 Subject: [PATCH 02/17] Add conservation-related utilities --- polaris/ocean/conservation.py | 65 ++++------------------------------- 1 file changed, 7 insertions(+), 58 deletions(-) diff --git a/polaris/ocean/conservation.py b/polaris/ocean/conservation.py index f48114d94b..dc1fe52b13 100644 --- a/polaris/ocean/conservation.py +++ b/polaris/ocean/conservation.py @@ -7,7 +7,7 @@ def compute_total_mass(ds_mesh, ds): - ds = reduce_dataset_time_dim(ds) + ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness total_mass = rho_sw * ( @@ -17,7 +17,7 @@ def compute_total_mass(ds_mesh, ds): def compute_total_mass_nonboussinesq(ds_mesh, ds): - ds = reduce_dataset_time_dim(ds) + ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness rho = ds.density @@ -28,7 +28,7 @@ def compute_total_mass_nonboussinesq(ds_mesh, ds): def compute_total_energy(ds_mesh, ds): - ds = reduce_dataset_time_dim(ds) + ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness temperature = ds.temperature @@ -43,17 +43,17 @@ def compute_total_energy(ds_mesh, ds): def compute_total_salt(ds_mesh, ds): - ds = reduce_dataset_time_dim(ds) + ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness salinity = ds.salinity - total_energy = ( + total_salt = ( area_cell * (layer_thickness * salinity).sum(dim='nVertLevels') ).sum(dim='nCells') - return total_energy + return total_salt -def reduce_dataset_time_dim(ds): +def _reduce_dataset_time_dim(ds): for time_dim in ['time', 'Time']: if time_dim in ds.dims: if ds.sizes[time_dim] > 1: @@ -64,54 +64,3 @@ def reduce_dataset_time_dim(ds): else: ds = ds.isel(time_dim=0) return ds - - -# def compute_net_mass_flux(ds): -# netMassFlux = & -# + accumulatedRainFlux & -# + accumulatedSnowFlux & -# + accumulatedEvaporationFlux & -# + accumulatedSeaIceFlux & -# + accumulatedRiverRunoffFlux & -# + accumulatedIceRunoffFlux & -# + accumulatedIcebergFlux & -# + accumulatedFrazilFlux & -# + accumulatedLandIceFlux -# note, accumulatedLandIceFrazilFlux not added because already in -# accumulatedFrazilFlux -# if (config_subglacial_runoff_mode == 'data'): -# netMassFlux = netMassFlux + accumulatedSubglacialRunoffFlux -# return net_mass_flux -# -# def compute_net_salt_flux(ds): -# netSaltFlux = accumulatedSeaIceSalinityFlux & -# + accumulatedFrazilSalinityFlux -# -# if (config_subglacial_runoff_mode == 'data'): -# netSaltFlux = netSaltFlux + accumulatedSubglacialRunoffSalinityFlux -# return net_salt_flux -# -# def compute_net_energy_flux(ds): -# netEnergyFlux = & -# accumulatedSeaIceHeatFlux & -# + accumulatedShortWaveHeatFlux & -# + accumulatedLongWaveHeatFluxUp & -# + accumulatedLongWaveHeatFluxDown & -# + accumulatedLatentHeatFlux & -# + accumulatedSensibleHeatFlux & -# + accumulatedMeltingSnowHeatFlux & -# + accumulatedMeltingIceRunoffHeatFlux & -# + accumulatedIcebergHeatFlux & -# + accumulatedFrazilHeatFlux & -# + accumulatedLandIceHeatFlux & -# + accumulatedRainTemperatureFlux *rho_sw*cp_sw & -# + accumulatedEvapTemperatureFlux *rho_sw*cp_sw & -# + accumulatedSeaIceTemperatureFlux *rho_sw*cp_sw & -# + accumulatedRiverRunoffTemperatureFlux *rho_sw*cp_sw & -# + accumulatedIcebergTemperatureFlux*rho_sw*cp_sw -# note, accumulatedLandIceFrazilHeatFlux not added because already in -# accumulatedFrazilHeatFlux -# if (config_subglacial_runoff_mode == 'data'): -# netEnergyFlux += accumulatedSubglacialRunoffTemperatureFlux * \ -# rho_sw*cp_sw -# return net_energy_flux From 9ba4df2d89fbb9fd6fb66c338f1664e4aa13d19b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 16:35:25 -0600 Subject: [PATCH 03/17] Test mass conservation for btrgyre --- polaris/tasks/ocean/barotropic_gyre/forward.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/polaris/tasks/ocean/barotropic_gyre/forward.py b/polaris/tasks/ocean/barotropic_gyre/forward.py index 6997ee36cf..e00246d696 100644 --- a/polaris/tasks/ocean/barotropic_gyre/forward.py +++ b/polaris/tasks/ocean/barotropic_gyre/forward.py @@ -99,11 +99,15 @@ def __init__( # make sure output is double precision self.add_yaml_file('polaris.ocean.config', 'output.yaml') + self.add_input_file( + filename='mesh.nc', target='../init/culled_mesh.nc' + ) self.add_input_file(filename='init.nc', target='../init/init.nc') self.add_output_file( filename='output.nc', validate_vars=['layerThickness', 'normalVelocity'], + verify_properties=['mass conservation'], ) self.package = 'polaris.tasks.ocean.barotropic_gyre' From 625d404a9894d60c61259a51794f657db0874687 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 16:37:58 -0600 Subject: [PATCH 04/17] Add other property tests to ocean_model_step --- polaris/ocean/model/ocean_model_step.py | 40 ++++++++++++++++++------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index eb8f39c024..88fbd0f3a7 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -6,7 +6,11 @@ from ruamel.yaml import YAML from polaris.model_step import ModelStep -from polaris.ocean.conservation import compute_total_mass +from polaris.ocean.conservation import ( + compute_total_energy, + compute_total_mass, + compute_total_salt, +) from polaris.tasks.ocean import Ocean OptionValue = Union[str, int, float, bool] @@ -386,14 +390,28 @@ def verify_properties(self): ) init_mass = compute_total_mass(ds_mesh, ds.isel(Time=0)) final_mass = compute_total_mass(ds_mesh, ds.isel(Time=-1)) - # mass_flux = compute_total_mass(ds.isel(Time=-1)) - result = abs(init_mass - final_mass) < tol - # absMassError = netMassFlux * dtAvg - massChange - # relMassError = absoluteMassError / (finalmass + 1.) - # absSaltError = netSaltFlux * dtAvg - saltChange - # relSaltError = absoluteSaltError / (finalSalt - 1.) - # absEnergyError = netEnergyFlux * dtAvg - energyChange - # relEnergyError = absoluteEnergyError / (finalEnergy - 1.) + mass_change = final_mass - init_mass + result = abs(mass_change) / (final_mass + 1.0) < tol + elif output_property == 'salt conservation': + tol = self.config.getfloat( + 'ocean', 'salt_conservation_tolerance' + ) + init_salt = compute_total_salt(ds_mesh, ds.isel(Time=0)) + final_salt = compute_total_salt(ds_mesh, ds.isel(Time=-1)) + salt_change = final_salt - init_salt + result = abs(salt_change) / (final_salt - 1.0) < tol + elif output_property == 'energy conservation': + tol = self.config.getfloat( + 'ocean', 'energy_conservation_tolerance' + ) + init_energy = compute_total_energy( + ds_mesh, ds.isel(Time=0) + ) + final_energy = compute_total_energy( + ds_mesh, ds.isel(Time=-1) + ) + energy_change = final_energy - init_energy + result = abs(energy_change) / (final_energy - 1.0) < tol else: raise ValueError( 'Could not find method to execute property check ' @@ -403,8 +421,8 @@ def verify_properties(self): success = success and result checked = True if not result: - failed_properties.extend(output_property) - + failed_properties.append(output_property) + print(failed_properties) if checked and success: log_filename = os.path.join( self.work_dir, 'property_check_passed.log' From 0d096eaea67281d0bb45fe081fcf023c66383741 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 16:39:25 -0600 Subject: [PATCH 05/17] Set conservation tolerances in ocean cfg section --- polaris/ocean/ocean.cfg | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/polaris/ocean/ocean.cfg b/polaris/ocean/ocean.cfg index c7c9357031..9827d3d1d9 100644 --- a/polaris/ocean/ocean.cfg +++ b/polaris/ocean/ocean.cfg @@ -4,9 +4,15 @@ # Options related the ocean component [ocean] -# Tolerance for mass conservation +# Tolerance for mass conservation, normalized by total mass mass_conservation_tolerance = 1e-8 +# Tolerance for salt conservation, normalized by total salt +salt_conservation_tolerance = 1e-8 + +# Tolerance for thermal energy conservation, normalized by total energy +energy_conservation_tolerance = 1e-8 + # Which model, MPAS-Ocean or Omega, is used, or "detect" to auto-detect model = detect From 335e69b46ee08a47ac40ad29014ee5925f09a61f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 16:47:10 -0600 Subject: [PATCH 06/17] Check that properties have been verified before running --- polaris/run/serial.py | 39 +++++++++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 4 deletions(-) diff --git a/polaris/run/serial.py b/polaris/run/serial.py index f7eab3a637..e9aeec4f25 100644 --- a/polaris/run/serial.py +++ b/polaris/run/serial.py @@ -466,6 +466,29 @@ def _read_baseline_status_from_logs(step_work_dir: str) -> Optional[bool]: return None +def _read_property_status_from_logs(step_work_dir: str) -> Optional[bool]: + """Get property check status from existing log markers. + + Returns + ------- + Optional[bool] + True if ``property_check_passed.log`` exists, False if + ``property_check_failed.log`` exists, otherwise None. + """ + property_check_pass_filename = os.path.join( + step_work_dir, 'property_check_passed.log' + ) + property_check_fail_filename = os.path.join( + step_work_dir, 'property_check_failed.log' + ) + + if os.path.exists(property_check_pass_filename): + return True + if os.path.exists(property_check_fail_filename): + return False + return None + + def _accumulate_baselines( baselines_passed: Optional[bool], status: bool ) -> Optional[bool]: @@ -486,6 +509,7 @@ def _run_task(task, available_resources): logger = task.logger cwd = os.getcwd() baselines_passed = None + property_passed = None for step_name in task.steps_to_run: step = task.steps[step_name] complete_filename = os.path.join( @@ -506,6 +530,16 @@ def _run_task(task, available_resources): baselines_passed = _accumulate_baselines( baselines_passed, baseline_status ) + property_status = None + property_status = _read_property_status_from_logs(step.work_dir) + if property_status is not None: + property_str = pass_str if property_status else fail_str + _print_to_stdout( + task, f' property comp.: {property_str}' + ) + property_passed = _accumulate_baselines( + property_passed, property_status + ) continue if step.cached: _print_to_stdout(task, ' cached') @@ -552,10 +586,7 @@ def _run_task(task, available_resources): _print_to_stdout( task, f' property verification: {property_str}' ) - if property_passed is None: - property_passed = status - elif not status: - property_passed = False + property_passed = _accumulate_baselines(property_passed, status) compared, status = step.validate_baselines() if compared: From 68e30010e2d09018b2e244f193def866e33c0086 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 10:15:01 -0600 Subject: [PATCH 07/17] Consistently use 'check' properties instead of 'verify' --- polaris/ocean/model/ocean_model_step.py | 6 +++--- polaris/run/serial.py | 4 ++-- polaris/step.py | 12 ++++++------ polaris/tasks/ocean/barotropic_gyre/forward.py | 2 +- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index 88fbd0f3a7..5bfe44024c 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -368,16 +368,16 @@ def update_namelist_eos(self) -> None: options=replacements, config_model='ocean' ) - def verify_properties(self): + def check_properties(self): checked = False success = True if self.work_dir is None: raise ValueError( 'The work directory must be set before the step ' - 'output properties can be verified.' + 'output properties can be checked.' ) failed_properties = [] - for filename, properties in self.properties_to_verify.items(): + for filename, properties in self.properties_to_check.items(): filename = str(filename) mesh_filename = os.path.join(self.work_dir, 'mesh.nc') this_filename = os.path.join(self.work_dir, filename) diff --git a/polaris/run/serial.py b/polaris/run/serial.py index e9aeec4f25..27b97ece7b 100644 --- a/polaris/run/serial.py +++ b/polaris/run/serial.py @@ -577,14 +577,14 @@ def _run_task(task, available_resources): step_time = time.time() - step_start step_time_str = str(timedelta(seconds=round(step_time))) - compared, status = step.verify_properties() + compared, status = step.check_properties() if compared: if status: property_str = pass_str else: property_str = fail_str _print_to_stdout( - task, f' property verification: {property_str}' + task, f' property checks: {property_str}' ) property_passed = _accumulate_baselines(property_passed, status) diff --git a/polaris/step.py b/polaris/step.py index 2006abb8c4..9b9f13193c 100644 --- a/polaris/step.py +++ b/polaris/step.py @@ -262,7 +262,7 @@ def __init__( # may be set during setup if there is a baseline for comparison self.baseline_dir = None self.validate_vars = dict() - self.properties_to_verify = dict() + self.properties_to_check = dict() self.setup_complete = False # these will be set before running the step, dummy placeholders for now @@ -479,7 +479,7 @@ def add_input_file( ) def add_output_file( - self, filename, validate_vars=None, verify_properties=None + self, filename, validate_vars=None, check_properties=None ): """ Add the output file that must be produced by this step and may be made @@ -502,8 +502,8 @@ def add_output_file( self.outputs.append(filename) if validate_vars is not None: self.validate_vars[filename] = validate_vars - if verify_properties is not None: - self.properties_to_verify[filename] = verify_properties + if check_properties is not None: + self.properties_to_check[filename] = check_properties def add_dependency(self, step, name=None): """ @@ -539,9 +539,9 @@ def add_dependency(self, step, name=None): target = f'{step.path}/step_after_run.pickle' self.add_input_file(filename=filename, work_dir_target=target) - def verify_properties(self): + def check_properties(self): """ - This method should be overridden to verify properties of step outputs + This method should be overridden to check properties of step outputs """ checked = False success = True diff --git a/polaris/tasks/ocean/barotropic_gyre/forward.py b/polaris/tasks/ocean/barotropic_gyre/forward.py index e00246d696..9efe0dcb56 100644 --- a/polaris/tasks/ocean/barotropic_gyre/forward.py +++ b/polaris/tasks/ocean/barotropic_gyre/forward.py @@ -107,7 +107,7 @@ def __init__( self.add_output_file( filename='output.nc', validate_vars=['layerThickness', 'normalVelocity'], - verify_properties=['mass conservation'], + check_properties=['mass conservation'], ) self.package = 'polaris.tasks.ocean.barotropic_gyre' From 4ac8032426c770a69e94885be91f4390b9cbb93f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 17 Dec 2025 17:03:22 -0600 Subject: [PATCH 08/17] Update the docs and doc strings --- docs/developers_guide/framework/validation.md | 48 +++++++++++++++++++ .../adding_forward.md | 3 ++ polaris/ocean/conservation.py | 14 ++++++ polaris/step.py | 3 ++ 4 files changed, 68 insertions(+) diff --git a/docs/developers_guide/framework/validation.md b/docs/developers_guide/framework/validation.md index c4047152a8..6e47c07d16 100644 --- a/docs/developers_guide/framework/validation.md +++ b/docs/developers_guide/framework/validation.md @@ -211,3 +211,51 @@ which need to have their variables renamed to the MPAS-Ocean names for use in Polaris. The `ds1` and `ds2` keyword arguments are used to supply datasets corresponding to `filename1` and `filename2`, respectively, in such circumstances. + +# Property checks + +For some output files, you may wish to run checks of certain properties such as +conservation of mass or energy. Currently, only conservation checks for the +ocean are available. + +To run property checks, add a list of properties +the keyword argument `verify_properties` to the +:py:meth:`polaris.Step.add_output_file()` method. As an example: + +```python +from polaris import Step + + +class Forward(OceanModelStep): + def __init__(self, task): + super().__init__(task=task, name='forward') + + self.add_input_file( + filename='mesh.nc', target='../init/culled_mesh.nc' + ) + self.add_output_file('output.nc', + verify_properties=['mass conservation']) +``` + +## Conservation checks + +Mass, salt and energy conservation checks are available for ocean model output. +The checks will fail if the relative change in the total quantity exceeds the +tolerance given in + +```cfg +# Options related the ocean component +[ocean] + +# Tolerance for mass conservation, normalized by total mass +mass_conservation_tolerance = 1e-8 + +# Tolerance for salt conservation, normalized by total salt +salt_conservation_tolerance = 1e-8 + +# Tolerance for thermal energy conservation, normalized by total energy +energy_conservation_tolerance = 1e-8 +``` + +As shown in the previous example, we have added a mesh file with the name +'mesh.nc' because conservation checks require the area of cells. diff --git a/docs/tutorials/dev_add_category_of_tasks/adding_forward.md b/docs/tutorials/dev_add_category_of_tasks/adding_forward.md index e93b2f825b..4ad2d17369 100644 --- a/docs/tutorials/dev_add_category_of_tasks/adding_forward.md +++ b/docs/tutorials/dev_add_category_of_tasks/adding_forward.md @@ -139,6 +139,9 @@ class Forward(OceanModelStep): 'normalVelocity', 'temperature', ], + check_properties=[ + 'mass conservation', + ], ) ``` diff --git a/polaris/ocean/conservation.py b/polaris/ocean/conservation.py index dc1fe52b13..4279a3c13a 100644 --- a/polaris/ocean/conservation.py +++ b/polaris/ocean/conservation.py @@ -7,6 +7,10 @@ def compute_total_mass(ds_mesh, ds): + """ + Compute the total mass in an ocean model output file using a constant + density + """ ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness @@ -17,6 +21,10 @@ def compute_total_mass(ds_mesh, ds): def compute_total_mass_nonboussinesq(ds_mesh, ds): + """ + Compute the total mass in an ocean model output file using the density + field given in the output file + """ ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness @@ -28,6 +36,9 @@ def compute_total_mass_nonboussinesq(ds_mesh, ds): def compute_total_energy(ds_mesh, ds): + """ + Compute the total heat content an ocean model output file + """ ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness @@ -43,6 +54,9 @@ def compute_total_energy(ds_mesh, ds): def compute_total_salt(ds_mesh, ds): + """ + Compute the total salt in an ocean model output file + """ ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness diff --git a/polaris/step.py b/polaris/step.py index 9b9f13193c..6112de1c2e 100644 --- a/polaris/step.py +++ b/polaris/step.py @@ -124,6 +124,9 @@ class Step: comparison should be performed if a baseline run has been provided. The baseline validation is performed after the step has run. + properties_to_check: dict of list + A list of properties to check for each output file. + logger : logging.Logger A logger for output from the step From cf845c203c7c38b2f08a798b8496d165039730c0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 10:27:38 -0600 Subject: [PATCH 09/17] Placeholder non-boussinesq mass conservation check for omega --- polaris/ocean/model/ocean_model_step.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index 5bfe44024c..8f32d420d4 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -9,6 +9,7 @@ from polaris.ocean.conservation import ( compute_total_energy, compute_total_mass, + # compute_total_mass_nonboussinesq, # Add when Omega EOS is used compute_total_salt, ) from polaris.tasks.ocean import Ocean @@ -388,6 +389,15 @@ def check_properties(self): tol = self.config.getfloat( 'ocean', 'mass_conservation_tolerance' ) + # Add when Omega EOS is used + # if self.config.get('ocean', 'model') == 'omega': + # init_mass = compute_total_mass_nonboussinesq( + # ds_mesh, ds.isel(Time=0) + # ) + # final_mass = compute_total_mass_nonboussinesq( + # ds_mesh, ds.isel(Time=-1) + # ) + # else: init_mass = compute_total_mass(ds_mesh, ds.isel(Time=0)) final_mass = compute_total_mass(ds_mesh, ds.isel(Time=-1)) mass_change = final_mass - init_mass From 34b123b4146946255fa0a6faee707051f988059a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 14:04:18 -0600 Subject: [PATCH 10/17] Fixup unrelated docs --- docs/developers_guide/api.md | 1 + docs/developers_guide/ocean/api.md | 15 +++++++++++++++ .../ocean/tasks/barotropic_channel.md | 4 ++-- .../ocean/tasks/customizable_viz.md | 4 ++-- 4 files changed, 20 insertions(+), 4 deletions(-) diff --git a/docs/developers_guide/api.md b/docs/developers_guide/api.md index 158622b82d..36b5435662 100644 --- a/docs/developers_guide/api.md +++ b/docs/developers_guide/api.md @@ -161,6 +161,7 @@ seaice/api Step.add_output_file Step.add_dependency Step.validate_baselines + Step.check_properties Step.set_shared_config ``` diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md index 229433602c..4e7872d67d 100644 --- a/docs/developers_guide/ocean/api.md +++ b/docs/developers_guide/ocean/api.md @@ -461,6 +461,20 @@ ## Ocean Framework +### Conservation utilities + +```{eval-rst} +.. currentmodule:: polaris.ocean.conservation + +.. autosummary:: + :toctree: generated/ + + compute_total_mass + compute_total_mass_nonboussinesq + compute_total_salt + compute_total_energy +``` + ### Convergence Tests ```{eval-rst} @@ -533,6 +547,7 @@ OceanModelStep OceanModelStep.setup + OceanModelStep.check_properties OceanModelStep.constrain_resources OceanModelStep.compute_cell_count OceanModelStep.map_yaml_options diff --git a/docs/developers_guide/ocean/tasks/barotropic_channel.md b/docs/developers_guide/ocean/tasks/barotropic_channel.md index f585d692e7..0160354b48 100644 --- a/docs/developers_guide/ocean/tasks/barotropic_channel.md +++ b/docs/developers_guide/ocean/tasks/barotropic_channel.md @@ -32,7 +32,7 @@ stream is also generated with wind stress fields. The class {py:class}`polaris.tasks.ocean.barotropic_channel.forward.Forward` defines a step for running the ocean from the initial condition produced in the `init` step. Namelist and streams files are updated in -{py:meth}`polaris.tasks.ocean.overflow.forward.Forward.dynamic_model_config()`. +{py:meth}`polaris.tasks.ocean.barotropic_channel.forward.Forward.dynamic_model_config()`. The number of cells is approximated from config options in {py:meth}`polaris.tasks.ocean.barotropic_channel.forward.Forward.compute_cell_count()` so that this can be used to constrain the number of MPI tasks that Polaris @@ -47,7 +47,7 @@ The {py:class}`polaris.tasks.ocean.barotropic_channel.viz.Viz` plots the initial final velocity components, relative vorticity, and circulation at the bottommost vertical layer. -(dev-ocean-overflow-default)= +(dev-ocean-barotropic-channel-default)= ## default diff --git a/docs/developers_guide/ocean/tasks/customizable_viz.md b/docs/developers_guide/ocean/tasks/customizable_viz.md index faab757f9d..73f796a07e 100644 --- a/docs/developers_guide/ocean/tasks/customizable_viz.md +++ b/docs/developers_guide/ocean/tasks/customizable_viz.md @@ -1,4 +1,4 @@ -(dev-ocean-customizable_viz)= +(dev-ocean-customizable-viz)= # customizable_viz @@ -10,7 +10,7 @@ controlled via the task configuration. ## framework The config options for the `customizable_viz` tests are described in -{ref}`ocean-customizable_viz` in the User's Guide. +{ref}`ocean-customizable-viz` in the User's Guide. The test also makes use of `polaris.viz.helper` and `polaris.viz.mplstyle`. From 1898d844947739b6b0ac033138cbe874614498eb Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 16:41:56 -0600 Subject: [PATCH 11/17] Generalize salt conservation check to all tracers --- polaris/ocean/conservation.py | 17 +++++++++++----- polaris/ocean/model/ocean_model_step.py | 27 +++++++++++++++++++++++++ polaris/ocean/ocean.cfg | 3 +++ 3 files changed, 42 insertions(+), 5 deletions(-) diff --git a/polaris/ocean/conservation.py b/polaris/ocean/conservation.py index 4279a3c13a..26697c2a65 100644 --- a/polaris/ocean/conservation.py +++ b/polaris/ocean/conservation.py @@ -53,18 +53,25 @@ def compute_total_energy(ds_mesh, ds): return total_energy -def compute_total_salt(ds_mesh, ds): +def compute_total_tracer(ds_mesh, ds, tracer_name='tracer1'): """ Compute the total salt in an ocean model output file """ ds = _reduce_dataset_time_dim(ds) area_cell = ds_mesh.areaCell layer_thickness = ds.layerThickness - salinity = ds.salinity - total_salt = ( - area_cell * (layer_thickness * salinity).sum(dim='nVertLevels') + tracer = ds[tracer_name] + total_tracer = ( + area_cell * (layer_thickness * tracer).sum(dim='nVertLevels') ).sum(dim='nCells') - return total_salt + return total_tracer + + +def compute_total_salt(ds_mesh, ds): + """ + Compute the total salt in an ocean model output file + """ + return compute_total_tracer(ds_mesh, ds, tracer_name='salinity') def _reduce_dataset_time_dim(ds): diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index 8f32d420d4..9859260a1b 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -11,6 +11,7 @@ compute_total_mass, # compute_total_mass_nonboussinesq, # Add when Omega EOS is used compute_total_salt, + compute_total_tracer, ) from polaris.tasks.ocean import Ocean @@ -410,6 +411,32 @@ def check_properties(self): final_salt = compute_total_salt(ds_mesh, ds.isel(Time=-1)) salt_change = final_salt - init_salt result = abs(salt_change) / (final_salt - 1.0) < tol + elif output_property == 'tracer conservation': + tol = self.config.getfloat( + 'ocean', 'tracer_conservation_tolerance' + ) + # Loop through all tracers in mpaso_to_omega.yaml + tracers_to_check = [ + 'temperature', + 'salinity', + 'tracer1', + 'tracer2', + 'tracer3', + ] + for tracer in tracers_to_check: + if tracer in ds.keys(): + init_tracer = compute_total_tracer( + ds_mesh, ds.isel(Time=0), tracer_name=tracer + ) + final_tracer = compute_total_tracer( + ds_mesh, ds.isel(Time=-1), tracer_name=tracer + ) + tracer_change = final_tracer - init_tracer + result = result and ( + abs(tracer_change) / (final_tracer - 1.0) < tol + ) + if result: + continue elif output_property == 'energy conservation': tol = self.config.getfloat( 'ocean', 'energy_conservation_tolerance' diff --git a/polaris/ocean/ocean.cfg b/polaris/ocean/ocean.cfg index 9827d3d1d9..6c68b5f648 100644 --- a/polaris/ocean/ocean.cfg +++ b/polaris/ocean/ocean.cfg @@ -10,6 +10,9 @@ mass_conservation_tolerance = 1e-8 # Tolerance for salt conservation, normalized by total salt salt_conservation_tolerance = 1e-8 +# Tolerance for tracer conservation, normalized by total tracer value +tracer_conservation_tolerance = 1e-8 + # Tolerance for thermal energy conservation, normalized by total energy energy_conservation_tolerance = 1e-8 From b7b19b46ceb069c5acfda9a7dfb9cf935cc4297c Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 16:46:57 -0600 Subject: [PATCH 12/17] Support conservation checks for convergence tests --- polaris/ocean/convergence/forward.py | 8 +++++++- polaris/tasks/ocean/cosine_bell/forward.py | 1 + polaris/tasks/ocean/sphere_transport/forward.py | 1 + 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py index 58923a8a7a..8355572cef 100644 --- a/polaris/ocean/convergence/forward.py +++ b/polaris/ocean/convergence/forward.py @@ -38,6 +38,7 @@ def __init__( graph_target=None, output_filename='output.nc', validate_vars=None, + check_properties=None, refinement='both', ): """ @@ -108,9 +109,14 @@ def __init__( self.add_input_file( filename='init.nc', work_dir_target=f'{init.path}/initial_state.nc' ) + self.add_input_file( + filename='mesh.nc', work_dir_target=f'{init.path}/mesh.nc' + ) self.add_output_file( - filename=output_filename, validate_vars=validate_vars + filename=output_filename, + validate_vars=validate_vars, + check_properties=check_properties, ) def compute_cell_count(self): diff --git a/polaris/tasks/ocean/cosine_bell/forward.py b/polaris/tasks/ocean/cosine_bell/forward.py index 989ee83bb4..bef99c4887 100644 --- a/polaris/tasks/ocean/cosine_bell/forward.py +++ b/polaris/tasks/ocean/cosine_bell/forward.py @@ -65,6 +65,7 @@ def __init__( yaml_filename='forward.yaml', output_filename='output.nc', validate_vars=validate_vars, + check_properties=['mass conservation', 'tracer conservation'], graph_target=f'{mesh.path}/graph.info', refinement_factor=refinement_factor, refinement=refinement, diff --git a/polaris/tasks/ocean/sphere_transport/forward.py b/polaris/tasks/ocean/sphere_transport/forward.py index 695386bdac..e784aeafb6 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.py +++ b/polaris/tasks/ocean/sphere_transport/forward.py @@ -71,6 +71,7 @@ def __init__( yaml_filename='forward.yaml', output_filename='output.nc', validate_vars=validate_vars, + check_properties=['mass conservation', 'tracer conservation'], options=namelist_options, graph_target=f'{base_mesh.path}/graph.info', refinement_factor=refinement_factor, From b688c76ba32d4fcbdc5156b8ee2c3dd50238cac1 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 16:47:32 -0600 Subject: [PATCH 13/17] Check mass conservation for manufactured_solution test --- polaris/tasks/ocean/manufactured_solution/forward.py | 1 + 1 file changed, 1 insertion(+) diff --git a/polaris/tasks/ocean/manufactured_solution/forward.py b/polaris/tasks/ocean/manufactured_solution/forward.py index 493e4e0cab..ee325a2bc4 100644 --- a/polaris/tasks/ocean/manufactured_solution/forward.py +++ b/polaris/tasks/ocean/manufactured_solution/forward.py @@ -86,6 +86,7 @@ def __init__( graph_target=f'{init.path}/culled_graph.info', output_filename='output.nc', validate_vars=['layerThickness', 'normalVelocity'], + check_properties=['mass conservation'], ) self.del2 = del2 self.del4 = del4 From 275728d405dea36eac980fb0131550838c237c04 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 17:31:17 -0600 Subject: [PATCH 14/17] Add more information to property check logs --- polaris/ocean/model/ocean_model_step.py | 42 ++++++++++++++++++++----- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index 9859260a1b..34e428693d 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -378,6 +378,7 @@ def check_properties(self): 'The work directory must be set before the step ' 'output properties can be checked.' ) + passed_properties = [] failed_properties = [] for filename, properties in self.properties_to_check.items(): filename = str(filename) @@ -402,7 +403,8 @@ def check_properties(self): init_mass = compute_total_mass(ds_mesh, ds.isel(Time=0)) final_mass = compute_total_mass(ds_mesh, ds.isel(Time=-1)) mass_change = final_mass - init_mass - result = abs(mass_change) / (final_mass + 1.0) < tol + relative_error = abs(mass_change) / (final_mass + 1.0) + result = relative_error < tol elif output_property == 'salt conservation': tol = self.config.getfloat( 'ocean', 'salt_conservation_tolerance' @@ -410,7 +412,8 @@ def check_properties(self): init_salt = compute_total_salt(ds_mesh, ds.isel(Time=0)) final_salt = compute_total_salt(ds_mesh, ds.isel(Time=-1)) salt_change = final_salt - init_salt - result = abs(salt_change) / (final_salt - 1.0) < tol + relative_error = abs(salt_change) / (final_salt - 1.0) + result = relative_error < tol elif output_property == 'tracer conservation': tol = self.config.getfloat( 'ocean', 'tracer_conservation_tolerance' @@ -432,9 +435,20 @@ def check_properties(self): ds_mesh, ds.isel(Time=-1), tracer_name=tracer ) tracer_change = final_tracer - init_tracer - result = result and ( - abs(tracer_change) / (final_tracer - 1.0) < tol + relative_error = abs(tracer_change) / ( + final_tracer - 1.0 ) + result = relative_error < tol + if not result: + failed_properties.append( + f'{tracer} conservation relative error ' + f'{relative_error:.3e} exceeds {tol}' + ) + else: + passed_properties.append( + f'{tracer} conservation relative error ' + f'{relative_error:.3e}' + ) if result: continue elif output_property == 'energy conservation': @@ -448,7 +462,8 @@ def check_properties(self): ds_mesh, ds.isel(Time=-1) ) energy_change = final_energy - init_energy - result = abs(energy_change) / (final_energy - 1.0) < tol + relative_error = abs(energy_change) / (final_energy - 1.0) + result = relative_error < tol else: raise ValueError( 'Could not find method to execute property check ' @@ -457,16 +472,27 @@ def check_properties(self): success = success and result checked = True - if not result: - failed_properties.append(output_property) - print(failed_properties) + # We already appended log strings for tracer conservation + if output_property != 'tracer conservation': + if not result: + failed_properties.append( + f'{output_property} relative error ' + f'{relative_error:.3e} exceeds {tol}' + ) + else: + passed_properties.append( + f'{output_property} relative error ' + f'{relative_error:.3e}' + ) if checked and success: log_filename = os.path.join( self.work_dir, 'property_check_passed.log' ) + passed_properties_str = '\n '.join(passed_properties) with open(log_filename, 'w') as result_log_file: result_log_file.write( f'Output file {filename} passed property checks.\n' + f'{passed_properties_str}\n' ) elif checked and not success: log_filename = os.path.join( From bb4547aa4526af4c54efc949598671e5aa70d75d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 18 Dec 2025 17:32:11 -0600 Subject: [PATCH 15/17] Make tolerances less permissive --- polaris/ocean/ocean.cfg | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/polaris/ocean/ocean.cfg b/polaris/ocean/ocean.cfg index 6c68b5f648..faa5bce94b 100644 --- a/polaris/ocean/ocean.cfg +++ b/polaris/ocean/ocean.cfg @@ -5,16 +5,16 @@ [ocean] # Tolerance for mass conservation, normalized by total mass -mass_conservation_tolerance = 1e-8 +mass_conservation_tolerance = 1e-14 # Tolerance for salt conservation, normalized by total salt -salt_conservation_tolerance = 1e-8 +salt_conservation_tolerance = 1e-14 # Tolerance for tracer conservation, normalized by total tracer value -tracer_conservation_tolerance = 1e-8 +tracer_conservation_tolerance = 1e-14 # Tolerance for thermal energy conservation, normalized by total energy -energy_conservation_tolerance = 1e-8 +energy_conservation_tolerance = 1e-14 # Which model, MPAS-Ocean or Omega, is used, or "detect" to auto-detect model = detect From 48e72c8596d7f4f5a6857bb3612519ae18e9cf7a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 19 Dec 2025 15:11:00 -0600 Subject: [PATCH 16/17] Clean-up ocean conservation --- polaris/ocean/conservation.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/polaris/ocean/conservation.py b/polaris/ocean/conservation.py index 26697c2a65..d34f75fecf 100644 --- a/polaris/ocean/conservation.py +++ b/polaris/ocean/conservation.py @@ -20,21 +20,6 @@ def compute_total_mass(ds_mesh, ds): return total_mass -def compute_total_mass_nonboussinesq(ds_mesh, ds): - """ - Compute the total mass in an ocean model output file using the density - field given in the output file - """ - ds = _reduce_dataset_time_dim(ds) - area_cell = ds_mesh.areaCell - layer_thickness = ds.layerThickness - rho = ds.density - total_mass = rho * ( - area_cell * (layer_thickness * rho).sum(dim='nVertLevels') - ).sum(dim='nCells') - return total_mass - - def compute_total_energy(ds_mesh, ds): """ Compute the total heat content an ocean model output file From 3fbdd8c446d4bed4ed3f0f94afa62da8d82188f0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 19 Dec 2025 15:10:34 -0600 Subject: [PATCH 17/17] Clean-up ocean model step --- polaris/ocean/model/ocean_model_step.py | 113 ++++++++++++------------ 1 file changed, 55 insertions(+), 58 deletions(-) diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index 34e428693d..46e2bea1a9 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -386,6 +386,24 @@ def check_properties(self): this_filename = os.path.join(self.work_dir, filename) ds_mesh = self.component.open_model_dataset(mesh_filename) ds = self.component.open_model_dataset(this_filename) + if 'tracer conservation' in properties: + # All tracers in mpaso_to_omega.yaml + tracers_to_check = [ + 'temperature', + 'salinity', + 'tracer1', + 'tracer2', + 'tracer3', + ] + # Expand 'tracer conservation' into list of tracers to check + properties = [ + item + for item in properties + if item != 'tracer conservation' + ] + for tracer in tracers_to_check: + if tracer in ds.keys(): + properties.append(f'tracer conservation-{tracer}') for output_property in properties: if output_property == 'mass conservation': tol = self.config.getfloat( @@ -393,83 +411,48 @@ def check_properties(self): ) # Add when Omega EOS is used # if self.config.get('ocean', 'model') == 'omega': - # init_mass = compute_total_mass_nonboussinesq( - # ds_mesh, ds.isel(Time=0) - # ) - # final_mass = compute_total_mass_nonboussinesq( - # ds_mesh, ds.isel(Time=-1) - # ) + # relative_error = self._compute_rel_err( + # compute_total_mass_nonboussinesq, ds_mesh, ds + # ) # else: - init_mass = compute_total_mass(ds_mesh, ds.isel(Time=0)) - final_mass = compute_total_mass(ds_mesh, ds.isel(Time=-1)) - mass_change = final_mass - init_mass - relative_error = abs(mass_change) / (final_mass + 1.0) - result = relative_error < tol + relative_error = self._compute_rel_err( + compute_total_mass, ds_mesh=ds_mesh, ds=ds + ) elif output_property == 'salt conservation': tol = self.config.getfloat( 'ocean', 'salt_conservation_tolerance' ) - init_salt = compute_total_salt(ds_mesh, ds.isel(Time=0)) - final_salt = compute_total_salt(ds_mesh, ds.isel(Time=-1)) - salt_change = final_salt - init_salt - relative_error = abs(salt_change) / (final_salt - 1.0) - result = relative_error < tol - elif output_property == 'tracer conservation': + relative_error = self._compute_rel_err( + compute_total_mass, ds_mesh, ds + ) + relative_error = self._compute_rel_err( + compute_total_salt, ds_mesh, ds + ) + elif output_property.split('-')[0] == 'tracer conservation': tol = self.config.getfloat( 'ocean', 'tracer_conservation_tolerance' ) - # Loop through all tracers in mpaso_to_omega.yaml - tracers_to_check = [ - 'temperature', - 'salinity', - 'tracer1', - 'tracer2', - 'tracer3', - ] - for tracer in tracers_to_check: - if tracer in ds.keys(): - init_tracer = compute_total_tracer( - ds_mesh, ds.isel(Time=0), tracer_name=tracer - ) - final_tracer = compute_total_tracer( - ds_mesh, ds.isel(Time=-1), tracer_name=tracer - ) - tracer_change = final_tracer - init_tracer - relative_error = abs(tracer_change) / ( - final_tracer - 1.0 - ) - result = relative_error < tol - if not result: - failed_properties.append( - f'{tracer} conservation relative error ' - f'{relative_error:.3e} exceeds {tol}' - ) - else: - passed_properties.append( - f'{tracer} conservation relative error ' - f'{relative_error:.3e}' - ) - if result: - continue + tracer = output_property.split('-')[1] + relative_error = self._compute_rel_err( + compute_total_tracer, + ds_mesh, + ds, + tracer_name=tracer, + ) elif output_property == 'energy conservation': tol = self.config.getfloat( 'ocean', 'energy_conservation_tolerance' ) - init_energy = compute_total_energy( - ds_mesh, ds.isel(Time=0) - ) - final_energy = compute_total_energy( - ds_mesh, ds.isel(Time=-1) + relative_error = self._compute_rel_err( + compute_total_energy, ds_mesh, ds ) - energy_change = final_energy - init_energy - relative_error = abs(energy_change) / (final_energy - 1.0) - result = relative_error < tol else: raise ValueError( 'Could not find method to execute property check ' f'{output_property}' ) + result = relative_error < tol success = success and result checked = True # We already appended log strings for tracer conservation @@ -691,6 +674,20 @@ def _map_mpaso_to_omega_section_option( return out_sections, out_option, out_value + def _compute_rel_err( + self, + func, + ds_mesh, + ds, + time_index_start=0, + time_index_end=-1, + **kwargs, + ): + init_val = func(ds_mesh, ds.isel(Time=time_index_start), **kwargs) + final_val = func(ds_mesh, ds.isel(Time=time_index_end), **kwargs) + val_change = final_val - init_val + return abs(val_change) / (final_val + 1.0) + @staticmethod def _warn_not_found(not_found: List[str]) -> None: """Warn about options that were not found in the map"""