From 85817036517b8802eb89254c93983659e16bd6d8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 12 Jan 2026 17:14:31 -0500 Subject: [PATCH 01/27] Use config option to determine whether to add coeffs_reconstruct to the stream --- polaris/ocean/config/coeffs_reconstruct.yaml | 9 +++++++++ polaris/ocean/model/ocean_model_step.py | 16 ++++++++++++++++ polaris/ocean/mpas_ocean.cfg | 6 ++++++ 3 files changed, 31 insertions(+) create mode 100644 polaris/ocean/config/coeffs_reconstruct.yaml diff --git a/polaris/ocean/config/coeffs_reconstruct.yaml b/polaris/ocean/config/coeffs_reconstruct.yaml new file mode 100644 index 0000000000..947508c07d --- /dev/null +++ b/polaris/ocean/config/coeffs_reconstruct.yaml @@ -0,0 +1,9 @@ +mpas-ocean: + streams: + coeffs_reconstruct: + type: output + filename_template: coeffs_reconstruct.nc + output_interval: 0000_00:00:01 + clobber_mode: truncate + contents: + - coeffs_reconstruct diff --git a/polaris/ocean/model/ocean_model_step.py b/polaris/ocean/model/ocean_model_step.py index fac215dbe9..1032f507c0 100644 --- a/polaris/ocean/model/ocean_model_step.py +++ b/polaris/ocean/model/ocean_model_step.py @@ -221,6 +221,22 @@ def dynamic_model_config(self, at_setup: bool) -> None: if self.update_eos: self.update_namelist_eos() + if self.config.has_option('ocean', 'write_coeffs_reconstruct'): + self.write_coeffs_reconstruct = self.config.get( + 'ocean', 'write_coeffs_reconstruct' + ) + else: + self.write_coeffs_reconstruct = False + if self.write_coeffs_reconstruct: + model = self.config.get('ocean', 'model') + if not model == 'mpas-ocean': + raise ValueError( + 'Coefficients for vector reconstruction can only be ' + 'written for ocean model MPAS-Ocean' + ) + self.add_yaml_file( + 'polaris.ocean.config', 'coeffs_reconstruct.yaml' + ) def constrain_resources(self, available_cores: Dict[str, Any]) -> None: """ diff --git a/polaris/ocean/mpas_ocean.cfg b/polaris/ocean/mpas_ocean.cfg index 89b4fda20a..d677fac751 100644 --- a/polaris/ocean/mpas_ocean.cfg +++ b/polaris/ocean/mpas_ocean.cfg @@ -1,5 +1,11 @@ # This config file has default config options for MPAS-Ocean +[ocean] + +# Whether to write coeffs for reconstructing vector quantities during +# the forward run +write_coeffs_reconstruct = False + # The paths section points polaris to external paths [paths] From 771813c7339e4f37a6332a59003ed96620256a5f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 26 Jan 2026 18:55:56 -0600 Subject: [PATCH 02/27] Remove later: test write coeffs --- polaris/ocean/mpas_ocean.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/polaris/ocean/mpas_ocean.cfg b/polaris/ocean/mpas_ocean.cfg index d677fac751..1641747582 100644 --- a/polaris/ocean/mpas_ocean.cfg +++ b/polaris/ocean/mpas_ocean.cfg @@ -4,7 +4,7 @@ # Whether to write coeffs for reconstructing vector quantities during # the forward run -write_coeffs_reconstruct = False +write_coeffs_reconstruct = True # The paths section points polaris to external paths [paths] From f8e5d4a897f15959a097333e98b731fc1975a0ec Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 3 Apr 2026 13:22:57 -0500 Subject: [PATCH 03/27] If model is Omega, link coeffs file from database --- polaris/mesh/base/__init__.py | 7 +++++++ polaris/tasks/ocean/barotropic_gyre/forward.py | 5 +++++ polaris/tasks/ocean/cosine_bell/forward.py | 11 +++++++++++ polaris/tasks/ocean/manufactured_solution/forward.py | 7 +++++++ polaris/tasks/ocean/merry_go_round/forward.py | 6 ++++++ polaris/tasks/ocean/sphere_transport/forward.py | 11 +++++++++++ 6 files changed, 47 insertions(+) diff --git a/polaris/mesh/base/__init__.py b/polaris/mesh/base/__init__.py index c5a8965c19..66293c253e 100644 --- a/polaris/mesh/base/__init__.py +++ b/polaris/mesh/base/__init__.py @@ -41,3 +41,10 @@ def get_base_mesh_steps(): base_mesh_steps.append(base_mesh_step) return base_mesh_steps + + +def parse_mesh_filepath(mesh_path): + component = mesh_path.split('/')[0] + database = '/'.join(mesh_path.split('/')[1:-1]) + mesh_name = mesh_path.split('/')[-1] + return component, database, mesh_name diff --git a/polaris/tasks/ocean/barotropic_gyre/forward.py b/polaris/tasks/ocean/barotropic_gyre/forward.py index 9efe0dcb56..56d8172f0c 100644 --- a/polaris/tasks/ocean/barotropic_gyre/forward.py +++ b/polaris/tasks/ocean/barotropic_gyre/forward.py @@ -277,6 +277,11 @@ def setup(self): # TODO: remove as soon as Omega no longer hard-codes this file if model == 'omega': self.add_input_file(filename='OmegaMesh.nc', target='init.nc') + self.add_input_file( + target='coeffs.nc', + filename='coeffs.nc', + database='barotropic_gyre', + ) def compute_max_time_step(self, config): """ diff --git a/polaris/tasks/ocean/cosine_bell/forward.py b/polaris/tasks/ocean/cosine_bell/forward.py index 2a0b844326..daac167c99 100644 --- a/polaris/tasks/ocean/cosine_bell/forward.py +++ b/polaris/tasks/ocean/cosine_bell/forward.py @@ -1,3 +1,4 @@ +from polaris.mesh.base import parse_mesh_filepath from polaris.ocean.convergence.spherical import SphericalConvergenceForward @@ -72,6 +73,7 @@ def __init__( refinement=refinement, ) self.do_restart = do_restart + self.mesh_path = mesh.path def setup(self): """ @@ -82,3 +84,12 @@ def setup(self): # TODO: remove as soon as Omega no longer hard-codes this file if model == 'omega': self.add_input_file(filename='OmegaMesh.nc', target='init.nc') + component, database, mesh_name = parse_mesh_filepath( + self.mesh_path + ) + self.add_input_file( + filename='coeffs.nc', + target=f'{mesh_name}_coeffs.nc', + database_component=component, + database=database, + ) diff --git a/polaris/tasks/ocean/manufactured_solution/forward.py b/polaris/tasks/ocean/manufactured_solution/forward.py index e226d63c9f..86e0d6620c 100644 --- a/polaris/tasks/ocean/manufactured_solution/forward.py +++ b/polaris/tasks/ocean/manufactured_solution/forward.py @@ -89,6 +89,7 @@ def __init__( validate_vars=['layerThickness', 'normalVelocity'], check_properties=['mass conservation'], ) + self.init = init self.del2 = del2 self.del4 = del4 @@ -102,6 +103,12 @@ def setup(self): # TODO: remove as soon as Omega no longer hard-codes this file if model == 'omega': self.add_input_file(filename='OmegaMesh.nc', target='init.nc') + mesh_name = self.init.path.split('/')[-1] + self.add_input_file( + target=f'{mesh_name}_coeffs.nc', + filename='coeffs.nc', + database='manufactured_solution', + ) def compute_cell_count(self): """ diff --git a/polaris/tasks/ocean/merry_go_round/forward.py b/polaris/tasks/ocean/merry_go_round/forward.py index 47756ad854..b1d301ba41 100644 --- a/polaris/tasks/ocean/merry_go_round/forward.py +++ b/polaris/tasks/ocean/merry_go_round/forward.py @@ -75,6 +75,7 @@ def __init__( ) self.order = vert_adv_order self.limiter = limiter + self.mesh_name = init.path.split('/')[-1] def setup(self): """ @@ -86,6 +87,11 @@ def setup(self): # TODO: remove as soon as Omega no longer hard-codes this file if model == 'omega': self.add_input_file(filename='OmegaMesh.nc', target='init.nc') + self.add_input_file( + target=f'{self.mesh_name}_coeffs.nc', + filename='coeffs.nc', + database='merry_go_round', + ) def dynamic_model_config(self, at_setup): """ diff --git a/polaris/tasks/ocean/sphere_transport/forward.py b/polaris/tasks/ocean/sphere_transport/forward.py index 6d72c57aa9..74f852593c 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.py +++ b/polaris/tasks/ocean/sphere_transport/forward.py @@ -1,3 +1,4 @@ +from polaris.mesh.base import parse_mesh_filepath from polaris.ocean.convergence.spherical import SphericalConvergenceForward @@ -78,6 +79,7 @@ def __init__( refinement_factor=refinement_factor, refinement=refinement, ) + self.mesh_path = base_mesh.path def setup(self): """ @@ -88,3 +90,12 @@ def setup(self): # TODO: remove as soon as Omega no longer hard-codes this file if model == 'omega': self.add_input_file(filename='OmegaMesh.nc', target='init.nc') + component, database, mesh_name = parse_mesh_filepath( + self.mesh_path + ) + self.add_input_file( + filename='coeffs.nc', + target=f'{mesh_name}_coeffs.nc', + database_component=component, + database=database, + ) From 64da2c1c2e1efef417ff3fe1fa4e793d2cb14cb6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 3 Apr 2026 14:59:58 -0500 Subject: [PATCH 04/27] Clean-up omega_pr suite order --- polaris/suites/ocean/omega_pr.txt | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/polaris/suites/ocean/omega_pr.txt b/polaris/suites/ocean/omega_pr.txt index 57c518c8fb..4e433af8bc 100644 --- a/polaris/suites/ocean/omega_pr.txt +++ b/polaris/suites/ocean/omega_pr.txt @@ -1,10 +1,17 @@ -#ocean/planar/barotropic_channel/default # Supported but fails -ocean/planar/merry_go_round/default ocean/planar/barotropic_gyre/munk/free-slip ocean/planar/manufactured_solution/convergence_both/default ocean/planar/manufactured_solution/convergence_both/del2 ocean/planar/manufactured_solution/convergence_both/del4 +ocean/planar/merry_go_round/default ocean/spherical/qu/rotation_2d ocean/spherical/icos/rotation_2d ocean/spherical/icos/cosine_bell/decomp ocean/spherical/icos/cosine_bell/restart + +# Supported but fails +#ocean/planar/barotropic_channel/default + +# Supported but expected to fail until higher-order advection is available +#ocean/spherical/icos/nondivergent_2d +#ocean/spherical/icos/divergent_2d +#ocean/spherical/icos/correlated_tracers_2d From 2e84e77f34bea925b485594cedc87c1157706ac4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 3 Apr 2026 16:40:03 -0500 Subject: [PATCH 05/27] Use mpas_tools reconstruct utility in open_model_dataset and test in sphere_transport viz --- polaris/tasks/ocean/__init__.py | 44 ++++++++++++++++++- .../sphere_transport/sphere_transport.cfg | 27 ++++++++++++ polaris/tasks/ocean/sphere_transport/viz.py | 25 ++++++++++- 3 files changed, 93 insertions(+), 3 deletions(-) diff --git a/polaris/tasks/ocean/__init__.py b/polaris/tasks/ocean/__init__.py index 243aa09f18..d5e9461f1b 100644 --- a/polaris/tasks/ocean/__init__.py +++ b/polaris/tasks/ocean/__init__.py @@ -4,6 +4,7 @@ import xarray as xr from mpas_tools.io import write_netcdf +from mpas_tools.vector.reconstruct import reconstruct_variable from ruamel.yaml import YAML from polaris import Component @@ -233,7 +234,14 @@ def map_var_list_from_native_model(self, var_list): ] return renamed_vars - def open_model_dataset(self, filename, **kwargs): + def open_model_dataset( + self, + filename, + ds_mesh=None, + reconstruct_variables=None, + coeffs_reconstruct=None, + **kwargs, + ): """ Open the given dataset, mapping variable and dimension names from Omega to MPAS-Ocean names if appropriate @@ -243,6 +251,15 @@ def open_model_dataset(self, filename, **kwargs): filename : str The path for the NetCDF file to open + ds_mesh : xr.Dataset + The MPAS mesh dataset for the ocean model. + + reconstruct_variables : list of str + List of variable names to reconstruct in the dataset. + + coeffs_reconstruct : xarray.DataArray, optional + Coefficients used for reconstructing variables. + kwargs keyword arguments passed to `xarray.open_dataset()` @@ -253,6 +270,31 @@ def open_model_dataset(self, filename, **kwargs): """ ds = xr.open_dataset(filename, **kwargs) ds = self.map_from_native_model_vars(ds) + if ds_mesh is None: + ds_mesh = xr.Dataset() + if reconstruct_variables is None: + reconstruct_variables = [] + if coeffs_reconstruct is None: + coeffs_reconstruct = xr.DataArray() + for variable in reconstruct_variables: + if variable in ds.keys(): + if 'normal' in variable: + out_var_name = variable.replace('normal', '').lower() + else: + out_var_name = variable + reconstruct_variable( + out_var_name, + ds[variable], + ds_mesh, + coeffs_reconstruct, + ds, + quiet=True, + ) + else: + raise ValueError( + f'User requested vector reconstruction for {variable} but ' + f"it isn't present in the dataset." + ) return ds def _has_ocean_io_model_steps(self, tasks) -> Tuple[bool, bool]: diff --git a/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg b/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg index b5a4da2a91..4f0432ecf2 100644 --- a/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg +++ b/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg @@ -138,3 +138,30 @@ norm_type = linear # A dictionary with keywords for the norm norm_args = {'vmin': -0.25, 'vmax': 0.25} + +# options for velocity visualization for the sphere transport test case +[sphere_transport_viz_velocity] + +# colormap options +# colormap +colormap_name = cmo.balance + +# the type of norm used in the colormap +norm_type = linear + +# A dictionary with keywords for the norm +norm_args = {'vmin': -50., 'vmax': 50.} + + +# options for plotting tracer differences from sphere transport tests +[sphere_transport_viz_velocity_diff] + +# colormap options +# colormap +colormap_name = cmo.balance + +# the type of norm used in the colormap +norm_type = linear + +# A dictionary with keywords for the norm +norm_args = {'vmin': -0.1, 'vmax': 0.1} diff --git a/polaris/tasks/ocean/sphere_transport/viz.py b/polaris/tasks/ocean/sphere_transport/viz.py index ef85d45b63..f0a5ce231d 100644 --- a/polaris/tasks/ocean/sphere_transport/viz.py +++ b/polaris/tasks/ocean/sphere_transport/viz.py @@ -1,5 +1,6 @@ import cmocean # noqa: F401 import numpy as np +import xarray as xr from polaris.mpas import time_since_start from polaris.ocean.model import OceanIOStep @@ -57,6 +58,9 @@ def __init__( self.add_input_file( filename='output.nc', work_dir_target=f'{forward.path}/output.nc' ) + self.add_input_file( + filename='coeffs.nc', work_dir_target=f'{forward.path}/coeffs.nc' + ) self.mesh_name = mesh_name variables_to_plot = dict( { @@ -64,6 +68,8 @@ def __init__( 'tracer2': 'tracer', 'tracer3': 'tracer', 'layerThickness': 'h', + 'velocityZonal': 'velocity', + 'velocityMeridional': 'velocity', } ) self.variables_to_plot = variables_to_plot @@ -82,10 +88,25 @@ def run(self): run_duration = config.getfloat('convergence_forward', 'run_duration') variables_to_plot = self.variables_to_plot - ds_init = self.open_model_dataset('initial_state.nc') + ds_mesh = xr.open_dataset('mesh.nc') + ds_coeff = xr.open_dataset('coeffs.nc') + coeffs_reconstruct = ds_coeff.coeffs_reconstruct + ds_init = self.open_model_dataset( + 'initial_state.nc', + decode_times=False, + ds_mesh=ds_mesh, + reconstruct_variables=['normalVelocity'], + coeffs_reconstruct=coeffs_reconstruct, + ) ds_init = ds_init[variables_to_plot.keys()].isel(Time=0, nVertLevels=0) - ds_out = self.open_model_dataset('output.nc', decode_times=False) + ds_out = self.open_model_dataset( + 'output.nc', + decode_times=False, + ds_mesh=ds_mesh, + reconstruct_variables=['normalVelocity'], + coeffs_reconstruct=coeffs_reconstruct, + ) s_per_hour = 3600.0 # Visualization at halfway around the globe (provided run duration is From 9fee43c29d8c7fe4a5eab62343e7059d3df23ba6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 3 Apr 2026 21:21:21 -0500 Subject: [PATCH 06/27] Ensure coeff reconstruction is not used when variables present --- polaris/tasks/ocean/__init__.py | 26 +++++++++++++++----- polaris/tasks/ocean/sphere_transport/init.py | 2 ++ polaris/tasks/ocean/sphere_transport/viz.py | 25 ++++++++++--------- 3 files changed, 36 insertions(+), 17 deletions(-) diff --git a/polaris/tasks/ocean/__init__.py b/polaris/tasks/ocean/__init__.py index d5e9461f1b..4f5dd2d86a 100644 --- a/polaris/tasks/ocean/__init__.py +++ b/polaris/tasks/ocean/__init__.py @@ -237,9 +237,9 @@ def map_var_list_from_native_model(self, var_list): def open_model_dataset( self, filename, - ds_mesh=None, + mesh_filename=None, reconstruct_variables=None, - coeffs_reconstruct=None, + coeffs_filename=None, **kwargs, ): """ @@ -270,18 +270,32 @@ def open_model_dataset( """ ds = xr.open_dataset(filename, **kwargs) ds = self.map_from_native_model_vars(ds) - if ds_mesh is None: - ds_mesh = xr.Dataset() if reconstruct_variables is None: reconstruct_variables = [] - if coeffs_reconstruct is None: - coeffs_reconstruct = xr.DataArray() for variable in reconstruct_variables: if variable in ds.keys(): if 'normal' in variable: out_var_name = variable.replace('normal', '').lower() else: out_var_name = variable + if ( + f'{out_var_name}Zonal' in ds.keys() + and f'{out_var_name}Meridional' in ds.keys() + ): + return ds + if mesh_filename is None: + raise ValueError( + 'mesh_filename must be provided to ' + f'open_model_dataset for reconstruction of {variable}' + ) + ds_mesh = xr.open_dataset(mesh_filename) + if coeffs_filename is None: + raise ValueError( + 'coeffs_filename must be provided to ' + f'open_model_dataset for reconstruction of {variable}' + ) + ds_coeff = xr.open_dataset(coeffs_filename) + coeffs_reconstruct = ds_coeff.coeffs_reconstruct reconstruct_variable( out_var_name, ds[variable], diff --git a/polaris/tasks/ocean/sphere_transport/init.py b/polaris/tasks/ocean/sphere_transport/init.py index 576daf7dd3..767e5eef4f 100644 --- a/polaris/tasks/ocean/sphere_transport/init.py +++ b/polaris/tasks/ocean/sphere_transport/init.py @@ -192,6 +192,8 @@ def run(self): ) normalVelocity, _ = xr.broadcast(normalVelocity, ds.refZMid) ds['normalVelocity'] = normalVelocity.expand_dims(dim='Time', axis=0) + ds['velocityZonal'] = np.nan * xr.zeros_like(ds.temperature) + ds['velocityMeridional'] = np.nan * xr.zeros_like(ds.temperature) ds['fCell'] = xr.zeros_like(ds_mesh.xCell) ds['fEdge'] = xr.zeros_like(ds_mesh.xEdge) diff --git a/polaris/tasks/ocean/sphere_transport/viz.py b/polaris/tasks/ocean/sphere_transport/viz.py index f0a5ce231d..8adf8de2c0 100644 --- a/polaris/tasks/ocean/sphere_transport/viz.py +++ b/polaris/tasks/ocean/sphere_transport/viz.py @@ -1,6 +1,5 @@ import cmocean # noqa: F401 import numpy as np -import xarray as xr from polaris.mpas import time_since_start from polaris.ocean.model import OceanIOStep @@ -58,9 +57,7 @@ def __init__( self.add_input_file( filename='output.nc', work_dir_target=f'{forward.path}/output.nc' ) - self.add_input_file( - filename='coeffs.nc', work_dir_target=f'{forward.path}/coeffs.nc' - ) + self.forward = forward self.mesh_name = mesh_name variables_to_plot = dict( { @@ -78,6 +75,15 @@ def __init__( self.add_output_file(f'{var}_final.png') self.add_output_file(f'{var}_diff.png') + def setup(self): + model = self.config.get('ocean', 'model') + # TODO: remove as soon as Omega no longer needs this file + if model == 'omega': + self.add_input_file( + filename='coeffs.nc', + work_dir_target=f'{self.forward.path}/coeffs.nc', + ) + def run(self): """ Run this step of the test case @@ -88,24 +94,21 @@ def run(self): run_duration = config.getfloat('convergence_forward', 'run_duration') variables_to_plot = self.variables_to_plot - ds_mesh = xr.open_dataset('mesh.nc') - ds_coeff = xr.open_dataset('coeffs.nc') - coeffs_reconstruct = ds_coeff.coeffs_reconstruct ds_init = self.open_model_dataset( 'initial_state.nc', decode_times=False, - ds_mesh=ds_mesh, + mesh_filename='mesh.nc', reconstruct_variables=['normalVelocity'], - coeffs_reconstruct=coeffs_reconstruct, + coeffs_filename='coeffs.nc', ) ds_init = ds_init[variables_to_plot.keys()].isel(Time=0, nVertLevels=0) ds_out = self.open_model_dataset( 'output.nc', decode_times=False, - ds_mesh=ds_mesh, + mesh_filename='mesh.nc', reconstruct_variables=['normalVelocity'], - coeffs_reconstruct=coeffs_reconstruct, + coeffs_filename='coeffs.nc', ) s_per_hour = 3600.0 From 9bd7e1ebd81b08e5681f0fa642b214aab330b27e Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Apr 2026 10:35:42 +0200 Subject: [PATCH 07/27] Update to v1.0.0-alpha.3 --- polaris/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/polaris/version.py b/polaris/version.py index e48e451469..8646659691 100644 --- a/polaris/version.py +++ b/polaris/version.py @@ -1 +1 @@ -__version__ = '1.0.0-alpha.2' +__version__ = '1.0.0-alpha.3' From 2eb1793f1510bdda6cb0536003ac461ee1d79c21 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Apr 2026 12:28:25 +0200 Subject: [PATCH 08/27] Update Omega submodule This update brings in: * https://github.com/E3SM-Project/Omega/pull/362 * https://github.com/E3SM-Project/Omega/pull/359 * https://github.com/E3SM-Project/Omega/pull/374 --- e3sm_submodules/Omega | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/e3sm_submodules/Omega b/e3sm_submodules/Omega index c409e9d53d..c2542a22e0 160000 --- a/e3sm_submodules/Omega +++ b/e3sm_submodules/Omega @@ -1 +1 @@ -Subproject commit c409e9d53dc116c7baefcf1bfcb740a5149ae786 +Subproject commit c2542a22e05b910b8688e10ed5df0521d9ddc4d6 From d6be8208c7a30f5ce65f82a3aef1075691f5b95d Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Apr 2026 15:14:38 +0200 Subject: [PATCH 09/27] Update to mache 3.3.0 --- deploy/cli_spec.json | 8 ++++---- deploy/pins.cfg | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/deploy/cli_spec.json b/deploy/cli_spec.json index 8bff7a8d52..7dde588c26 100644 --- a/deploy/cli_spec.json +++ b/deploy/cli_spec.json @@ -1,7 +1,7 @@ { "meta": { "software": "polaris", - "mache_version": "3.2.0", + "mache_version": "3.3.0", "description": "Deploy polaris environment" }, "arguments": [ @@ -18,9 +18,9 @@ "route": ["deploy", "bootstrap", "run"] }, { - "flags": ["--prefix"], - "dest": "prefix", - "help": "Install the environment into this prefix (directory). Overrides deploy/config.yaml.j2.", + "flags": ["--pixi-path", "--prefix"], + "dest": "pixi_path", + "help": "Install the pixi environment at this path (directory). Overrides deploy/config.yaml.j2. `--prefix` is deprecated.", "route": ["deploy", "bootstrap", "run"] }, { diff --git a/deploy/pins.cfg b/deploy/pins.cfg index a023be6966..dfaeafad54 100644 --- a/deploy/pins.cfg +++ b/deploy/pins.cfg @@ -3,7 +3,7 @@ bootstrap_python = 3.14 python = 3.14 geometric_features = 1.6.1 -mache = 3.2.0 +mache = 3.3.0 mpas_tools = 1.4.0 otps = 2021.10 parallelio = 2.6.9 From 136de47e879c70414a8b35383c71ba6f989ef60f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 4 Apr 2026 12:49:19 +0200 Subject: [PATCH 10/27] Add design doc for vector reconstruction --- docs/design_docs/index.md | 1 + docs/design_docs/vector_reconstruction.md | 666 ++++++++++++++++++++++ 2 files changed, 667 insertions(+) create mode 100644 docs/design_docs/vector_reconstruction.md diff --git a/docs/design_docs/index.md b/docs/design_docs/index.md index 64469c646a..525594c31d 100644 --- a/docs/design_docs/index.md +++ b/docs/design_docs/index.md @@ -6,5 +6,6 @@ :titlesonly: true shared_steps +vector_reconstruction template ``` diff --git a/docs/design_docs/vector_reconstruction.md b/docs/design_docs/vector_reconstruction.md new file mode 100644 index 0000000000..2a34710a6a --- /dev/null +++ b/docs/design_docs/vector_reconstruction.md @@ -0,0 +1,666 @@ +# Vector Reconstruction at MPAS Cell Centers + +date: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +## Summary + +This design proposes a new workflow for generating cell-center vector +reconstruction coefficients for MPAS meshes in Python within `polaris`. The +primary goal is for `polaris` to compute mesh-based reconstruction metadata and +coefficients, store them on an MPAS mesh, and enable Omega or MPAS-Ocean to +reconstruct vector fields at runtime in 3D Cartesian or zonal/meridional form. +Python is the intended home for coefficient generation because it provides the +most flexibility for experimenting with reconstruction methods, stencil choices, +and mesh-field layouts. + +The initial target field is `normalVelocity`, but the coefficient-generation +workflow should be general enough to apply to any field represented as a scalar +component normal to MPAS edges. The design follows the best-performing +practical methods identified in Peixoto and Barros (2014). The first +implementation should use linear least-squares (LSQ) reconstruction on all cells, +because offline coefficient generation in Python makes the setup cost less +important than simplicity, uniformity, and maintainability. A hybrid method +that combines Perot's inexpensive cell-centered reconstruction on well-aligned +cells with a linear least-squares reconstruction on poorly aligned cells should +be retained as a backup optimization path if runtime performance of the all-LSQ +approach proves insufficient. + +The design treats the reconstructed field as the sum of: + +$$ +\mathbf{u} = \mathbf{u}_t + w \hat{\mathbf{r}}, +$$ + +where $\mathbf{u}_t$ is tangent to the sphere and reconstructed from edge-normal +components, and $w \hat{\mathbf{r}}$ is an optional radial contribution aligned +with the local vertical unit vector. For velocity reconstruction, +`vertVelocityTop` may be supplied on cell interfaces and converted to layer +midpoints; if it is omitted, the radial contribution is assumed to be zero. + +Although velocity is the motivating use case, the core reconstruction should be +field-agnostic. The implementation will therefore separate: + +- a generic reconstruction of tangential vectors from edge-normal components, +- basis transformations between Cartesian and zonal/meridional coordinates, and +- generic 3D composition from a tangential component and an optional + interface-centered radial component. + +Success means that `polaris` gains a reusable, documented, and testable +coefficient-generation workflow in `polaris/mpas` that: + +- computes least-squares reconstruction coefficients and stencil metadata that + can be stored on an MPAS mesh and consumed at runtime by Omega or MPAS-Ocean, +- treats Python in `polaris` as the authoritative implementation for generating + these coefficients, +- reconstructs tangential vectors at MPAS cell centers with second-order + behavior on spherical MPAS meshes, +- supports outputs in 3D Cartesian and zonal/meridional form at cell centers + and layer midpoints, and +- is structured so later extension to arbitrary-point reconstruction can reuse + the same geometry and reconstruction cache. + +## Requirements + +### Requirement: Generic Tangential Reconstruction at Cell Centers + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +`polaris` shall support offline generation of reconstruction coefficients for a +tangential vector field at MPAS cell centers from scalar components given +normal to MPAS edges on spherical meshes. The capability must be general enough +to apply to fields other than velocity. + +The initial implementation shall focus on reconstruction at cell centers rather +than arbitrary points. The generated coefficients shall use only local mesh +geometry and local edge data, and they shall avoid methods that are known to be +less robust or more expensive than needed for second-order accuracy on MPAS +Voronoi meshes. + +The primary output of the implementation shall be coefficient and stencil data +that can be added to an MPAS mesh and used by Omega or MPAS-Ocean at runtime. + +### Requirement: 3D Vector Outputs at Layer Midpoints + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +`polaris` shall support a generic interface that returns reconstructed vector +fields at cell centers and layer midpoints in: + +- 3D Cartesian components, and +- zonal and meridional components. + +If a radial interface field is supplied on `nVertLevelsP1`, the reconstructed +vector shall include a radial contribution consistent with the MPAS local +vertical direction. If no radial interface field is supplied, the radial +contribution shall default to zero. + +For velocity, the corresponding radial interface field is `vertVelocityTop`. + +### Requirement: Reusable Precomputation and Extensibility + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +The implementation shall support precomputation of mesh-dependent reconstruction +metadata and weights so repeated application to many time slices or vertical +levels is efficient. The implementation shall be organized so later support for +reconstruction to arbitrary points can reuse the same geometry utilities, +stencil selection, and basis transforms without redesigning the core interface. + +The implementation shall define a mesh-storage format for reconstruction +coefficients and stencils that can be read directly at runtime by Omega or +MPAS-Ocean. + +The implementation shall not rely on online coefficient generation in +MPAS-Ocean or Omega. New reconstruction coefficients shall be generated in +Python in +`polaris` and then stored on the mesh for runtime use. + +## Algorithm Design + +### Algorithm Design: Generic Tangential Reconstruction at Cell Centers + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +Let $u_e$ be the scalar component of a tangential vector field normal to edge +$e$, evaluated at the edge midpoint. The goal is to reconstruct a tangential +vector $\mathbf{u}_t$ at a cell center $\mathbf{x}_c$. + +Following Peixoto and Barros (2014), the baseline method for the first +implementation should be linear least-squares reconstruction on all cells using +a Voronoi-centered two-ring stencil. This choice is motivated by the paper's +results for cell-centered and Voronoi-centered reconstruction: + +- least-squares reconstruction with a two-ring stencil is second-order and far + less sensitive to grid imprinting, +- it avoids cell-by-cell method switching and alignment-threshold tuning in the + primary implementation, and +- because coefficients are generated offline in Python, setup cost is less + important than a simple and uniform formulation. + +The same paper also supports a hybrid Perot/least-squares method as an +effective optimization. That hybrid should be retained as a fallback option if +runtime application of the all-LSQ coefficients proves too expensive. + +For MPAS cells, the reconstructed point is the cell center, interpreted as the +local reconstruction point. The implementation should use the local tangent +plane at the cell center, because the source data are tangent to the sphere and +the paper's spherical extensions are based on local projection to the tangent +plane. + +For a cell with local tangent coordinates $\boldsymbol{\xi}$ centered at the +cell center: + +1. Project nearby edge midpoints and edge-normal vectors into the tangent plane. +2. Apply least-squares reconstruction in that plane. +3. Map the reconstructed tangent vector back to 3D Cartesian coordinates. + +In the least-squares approach, a linear tangential vector field is fit in +tangent coordinates: + +$$ +\mathbf{u}_t(\boldsymbol{\xi}) = \mathbf{a}_0 + A \boldsymbol{\xi}, +$$ + +subject to the observed edge-normal constraints over a two-ring stencil. The +value at the cell center is then $\mathbf{a}_0$. As in the paper, the least +squares system should be solved in a weighted form to improve conditioning. + +The least-squares stencil should be Voronoi-centered and use two levels of +neighbors. In MPAS connectivity terms, this should be interpreted as the union +of `edgesOnVertex` over the `verticesOnCell` of the target cell, with duplicate +edges removed. Equivalently, this is the set of edges incident on the vertices +of the target Voronoi cell, not the much larger set obtained from +`edgesOnCell` of `cellsOnCell`. + +For a hexagon, each of the 6 vertices contributes one edge on the cell boundary +and one additional outer edge after duplicates are removed, giving 12 unique +edges. For a pentagon, the same construction gives 10 unique edges. This +matches the LS12 method described in the paper for cell-centered or +Voronoi-centered reconstruction. + +If runtime benchmarks later show that the all-LSQ approach is too expensive, a +secondary implementation path can adopt the paper's hybrid strategy. In that +case, Perot's method would be applied on well-aligned cells and least-squares +would be retained on the remaining cells. The Perot formula is: + +$$ +\mathbf{u}_{t,c}^{P} \approx +\frac{1}{A_c}\sum_{e \in c} l_e \, \mathbf{r}_e \, u_e, +$$ + +where $A_c$ is the cell area in the local projected plane, $l_e$ is edge +length, and $\mathbf{r}_e$ is the projected vector from the reconstruction +point to the midpoint of edge $e$. If this optimization path is pursued, the +alignment-index criterion from Peixoto and Barros should be revisited as part +of the performance study rather than being fixed in the first implementation. + +### Algorithm Design: 3D Vector Outputs at Layer Midpoints + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +The generic reconstruction produces only the tangential part of the vector. A +full 3D vector at layer midpoint $k$ should be constructed as: + +$$ +\mathbf{u}_{c,k} = +\mathbf{u}_{t,c,k} + w_{c,k} \hat{\mathbf{r}}_c, +$$ + +where $\hat{\mathbf{r}}_c$ is the local vertical unit vector at the cell center. + +If a radial interface field is provided on `nVertLevelsP1`, the midpoint radial +component should be formed by averaging adjacent interface values: + +$$ +w_{c,k} = \frac{1}{2}\left(w^\mathrm{top}_{c,k} + +w^\mathrm{top}_{c,k+1}\right). +$$ + +This yields a 3D vector defined at cell centers and layer midpoints. For +velocity, this is the expected target location for reconstructed 3D velocity in +MPAS-O analysis workflows. + +If no radial interface field is supplied, then: + +$$ +w_{c,k} = 0. +$$ + +The output should be available in both Cartesian and local geographic bases. +Using longitude $\lambda$ and latitude $\theta$ at the cell center, the +zonal-meridional decomposition should follow the same relationships used in the +MPAS framework: + +$$ +u_\mathrm{zonal} = -u_x \sin\lambda + u_y \cos\lambda, +$$ + +$$ +u_\mathrm{meridional} = +-(u_x \cos\lambda + u_y \sin\lambda)\sin\theta + u_z \cos\theta. +$$ + +These formulas are appropriate both for tangential-only reconstruction and for +the full vector after radial composition. + +### Algorithm Design: Reusable Precomputation and Extensibility + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +The reconstruction should be divided into a mesh-dependent setup phase and a +field-dependent apply phase. The setup phase is the primary deliverable for +`polaris`, because it produces the coefficients and stencil metadata to be +stored on the mesh. The setup phase should be implemented in Python and should +be considered the authoritative implementation for the new reconstruction +method. The apply phase is still useful in Python for testing, analysis, and workflow +development, while Omega or MPAS-Ocean will consume the stored coefficients at +runtime rather than recomputing them. + +The setup phase should precompute, for each cell: + +- the edge stencil, +- projected edge geometry in the local tangent plane, +- a weighted least-squares pseudo-inverse or equivalent coefficient matrix for + the cell-center value. + +The apply phase should then reduce to sparse local linear algebra. For a given +edge-normal field and vertical level, the tangential reconstruction at each cell +should be a weighted sum over the local stencil. + +This separation is important because: + +- the mesh geometry does not change between time slices, +- the same reconstruction is often applied over many vertical levels, and +- later arbitrary-point reconstruction can reuse the same tangent-plane + machinery, stencil construction, and basis transforms. + +The initial design intentionally excludes arbitrary-point interpolation with +spherical Wachspress coordinates. However, the reconstruction cache should be +built so this can be added later by introducing a second stage that evaluates +the field at vertices or other fixed locations and then interpolates. + +## Implementation + +### Implementation: Generic Tangential Reconstruction at Cell Centers + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +The proposed implementation adds a new coefficient-generation module to +`polaris/mpas`: + +- `polaris/mpas/reconstruct.py` + +and a companion module for basis transforms and 3D composition: + +- `polaris/mpas/vector.py` + +The public API should be field-agnostic at its core. A likely structure is: + +```python +def build_reconstruction_cache(ds_mesh, method="lsq"): + ... + + +def build_reconstruction_mesh_fields(ds_mesh, method="lsq"): + ... + + +def reconstruct_tangential_cell_center(edge_normal_field, ds_mesh, + cache=None): + ... + + +def tangential_cartesian_to_zonal_meridional(u_x, u_y, u_z, ds_mesh): + ... + + +def reconstruct_3d_cell_center(edge_normal_field, ds_mesh, + radial_interface_field=None, cache=None): + ... +``` + +`build_reconstruction_cache()` should: + +- validate that the mesh is spherical, +- determine the local tangent basis and local vertical direction at each cell, +- build the local stencil and coefficient arrays, and +- return a lightweight cache object, likely as an `xarray.Dataset` or a small + dataclass containing `xarray.DataArray` members. + +The default method should be all-cell least-squares. An optional `method` +argument may later support `hybrid`, but that should not be the baseline path. + +`build_reconstruction_mesh_fields()` should convert this cache into the mesh +variables needed by MPAS-Ocean or Omega at runtime. + +`reconstruct_tangential_cell_center()` should accept any edge-normal field with +an `nEdges` dimension, preserving all non-horizontal dimensions. This makes the +core implementation usable for velocity, pressure-gradient-like fields, or any +other tangential vector quantity represented through edge-normal components. + +`reconstruct_3d_cell_center()` should build on the tangential reconstruction by +adding an optional radial interface field and returning Cartesian and +zonal/meridional components. This function should be field-agnostic. + +For xarray compatibility, the Python reconstruction implementation should +vectorize over all +non-horizontal dimensions and apply the precomputed weights only along `nEdges`. +This will make the same code work naturally for: + +- `nEdges`, +- `nEdges, nVertLevels`, +- `Time, nEdges, nVertLevels`, + +and similar combinations. + +The least-squares path should avoid repeated dense solves at runtime. Instead, +the setup phase should precompute the operator that maps the stencil's edge +samples directly to the reconstructed Cartesian vector at the cell center. + +Pseudo code for setup: + +```python +for cell in cells: + basis = build_local_tangent_basis(cell) + stencil = build_two_ring_edge_stencil(cell) + coeffs[cell] = build_lsq_center_coeffs(cell, stencil, basis) +``` + +Pseudo code for apply: + +```python +for cell in cells: + tangential_xyz[..., cell] = sum( + coeffs[cell, i, :] * edge_normal[..., stencil[cell, i]] + for i in range(n_stencil[cell]) + ) +``` + +Because existing `polaris/mpas` modules are lightweight utilities, this design +fits the package style well. `polaris/mpas/__init__.py` should export the main +public functions. + +This design does not propose porting the new coefficient-generation +algorithm to Fortran for online execution in MPAS-Ocean, and it does not +propose computing these coefficients directly in Omega. Instead, the Python +implementation in `polaris` should generate the coefficients once for a given +mesh, after which runtime codes should only consume the stored mesh fields. + +The mesh-storage design needs special attention. Existing MPAS meshes already +contain: + +- `coeffs_reconstruct(R3, maxEdges, nCells)`, + +and the current MPAS runtime reconstruction uses `edgesOnCell` as the implicit +stencil. This is sufficient for the existing one-ring reconstruction but not +for the all-cell least-squares design described here, which requires a two-ring +stencil with up to 12 edges on hexagons and 10 on pentagons. + +Therefore, the least-squares design requires runtime-readable mesh fields +that explicitly encode both: + +- the reconstruction stencil for each cell, and +- the coefficients associated with that stencil. + +The recommended target representation is: + +- `reconstructEdgeStencil(maxEdges2, nCells)`, +- `nReconstructEdges(nCells)`, +- `coeffs_reconstruct(R3, maxEdges2, nCells)`. + +Here, `maxEdges2` is already an MPAS mesh dimension and is 12 on standard +hexagon-pentagon meshes, which matches the LS12 stencil width described in the +paper. This makes `maxEdges2` a natural candidate for storing the two-ring +least-squares stencil. + +This design explicitly recommends enlarging `coeffs_reconstruct` from +`maxEdges` to `maxEdges2` rather than introducing a second coefficient +variable. This keeps the runtime interface clearer by preserving the existing +coefficient name while expanding its stencil capacity to support both the +current one-ring case and the new two-ring least-squares case. + +### Implementation: 3D Vector Outputs at Layer Midpoints + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +`polaris/mpas/vector.py` should contain the generic 3D-composition logic that +sits on top of the tangential reconstruction and the coefficient-generation +workflow. + +The generic 3D wrapper should: + +1. reconstruct the tangential Cartesian vector from an edge-normal field, +2. broadcast the tangential result to all non-horizontal dimensions, +3. derive midpoint radial values from an optional radial interface field, +4. add the radial contribution in the local vertical direction, and +5. derive zonal and meridional components from the resulting Cartesian vector. + +For the optional radial interface field, the implementation should: + +- require an `nCells` dimension, +- require `nVertLevelsP1` if present, +- form midpoint values on `nVertLevels`, +- preserve any leading dimensions such as `Time`. + +If the radial interface field is omitted, the wrapper should behave exactly +like a horizontal-only reconstruction, with zero radial contribution. + +The generic wrapper should return an `xarray.Dataset` with clearly named output +fields. A likely default is: + +- `vectorX`, +- `vectorY`, +- `vectorZ`, +- `vectorZonal`, +- `vectorMeridional`. + +To align with MPAS naming and user expectations, the design should also discuss +how this generic output maps onto velocity-oriented runtime fields such as: + +- `uReconstructX`, +- `uReconstructY`, +- `uReconstructZ`, +- `uReconstructZonal`, +- `uReconstructMeridional`. + +The core implementation should remain field-agnostic even if some downstream +uses map the outputs onto velocity-specific names. + +At runtime, Omega or MPAS-Ocean should reconstruct: + +1. the tangential Cartesian vector from the stored stencil and coefficients, +2. the zonal and meridional components from the reconstructed Cartesian vector, +3. the full 3D vector by adding the radial contribution derived from an + optional radial interface field when needed. + +This keeps the coefficient storage generic while allowing field-specific +post-processing in the runtime model. Runtime codes should consume the stored +coefficients only; they should not be responsible for recomputing reconstruction +coefficients. + +### Implementation: Reusable Precomputation and Extensibility + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +The reconstruction cache should be serializable or at least representable as an +`xarray.Dataset`, so it can be inspected during debugging and written to a mesh +file. + +Recommended cache contents include: + +- `nStencilEdges` on `nCells`, +- `stencilEdges` on `nCells, maxStencilEdges`, +- `reconstructCoeffs` on `nCells, maxStencilEdges, R3`, +- local basis vectors on `nCells, R3`. + +This design intentionally precomputes coefficients for cell-center values only. +That keeps the initial implementation compact and directly aligned with the user +need. A later arbitrary-point design can either: + +- add a second cache for fixed target points such as vertices, or +- add target-specific coefficient builders that reuse the same tangent-plane and + stencil code. + +Project-management-wise, the work can proceed in four steps: + +1. implement the offline least-squares coefficient and stencil generator in + `polaris`, +2. define the MPAS/Omega mesh-field representation for those coefficients and + write them to mesh files, +3. implement basis transforms and the Python 3D reconstruction wrapper with an + optional radial interface field, +4. document the API and add examples for coefficient generation and + reconstruction from stored coefficients. + +If performance studies later indicate that runtime application is too costly, a +follow-on effort can evaluate the hybrid Perot/least-squares optimization. + +## Testing + +### Testing and Validation: Generic Tangential Reconstruction at Cell Centers + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +Testing should include unit tests and convergence tests. + +Unit tests should verify: + +- constant tangential vector fields reconstruct with small error, +- the reconstructed vector is tangent to the sphere when no radial component is + supplied, +- the cache and direct-setup paths give identical results, +- mesh-output fields contain the expected stencil width and coefficient layout. + +Convergence tests should use analytic spherical flows already available in +Polaris test cases, such as the sphere transport or cosine-bell style flows. +The test procedure should: + +1. evaluate analytic zonal/meridional velocity on edges, +2. convert to edge-normal components, +3. reconstruct at cell centers, +4. compare with the analytic vector at cell centers. + +The expected result is second-order convergence in RMS error for the +least-squares scheme on spherical MPAS meshes, with substantially reduced grid +imprinting relative to Perot alone. + +For the coefficient-generation path, tests should also confirm that applying the +stored mesh coefficients reproduces the same result as the direct Python +reconstruction. + +Performance tests should compare: + +- coefficient-generation cost in Python for representative mesh sizes, +- runtime application cost of the stored least-squares coefficients, +- mesh storage footprint for the widened stencil representation. + +If these tests show that runtime cost is too high, they should motivate +evaluation of the hybrid optimization path. + +### Testing and Validation: 3D Vector Outputs at Layer Midpoints + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +3D-vector tests should verify: + +- Cartesian and zonal/meridional outputs are mutually consistent, +- omission of the radial interface field yields zero radial contribution, +- midpoint radial velocity is computed correctly from interface values, +- shape and dimension handling are correct for `nVertLevels`, + `nVertLevelsP1`, and optional `Time`. + +A simple analytic test can define: + +- a known tangential horizontal velocity, +- a known radial profile on interfaces, +- an expected full Cartesian vector at layer midpoints. + +The generic wrapper should reproduce the known Cartesian and zonal/meridional +fields to within expected truncation error. A velocity-specific test case using +`normalVelocity` and `vertVelocityTop` should remain part of the test suite as +the motivating example. + +### Testing and Validation: Reusable Precomputation and Extensibility + +Date last modified: 2026/04/04 + +Contributors: + +- Xylar Asay-Davis +- Codex + +For the primary workflow, tests should also verify that mesh-writing preserves +the coefficient fields exactly and that a runtime-style application of those +fields gives identical answers before and after I/O. + +Documentation tests or example scripts should also be added to demonstrate: + +- generic tangential reconstruction from an arbitrary edge-normal field, +- velocity reconstruction from `normalVelocity`, +- velocity reconstruction from `normalVelocity` and `vertVelocityTop`. + +The initial implementation does not need regression-suite integration beyond +unit and convergence tests, but if the functionality becomes part of standard +analysis or initialization workflows, a lightweight analysis regression test +should be added later to guard against changes in reconstruction coefficients or +basis transforms. From 024b5b7b9efc351a2ee0648307f811a078728f23 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 4 Apr 2026 13:00:07 +0200 Subject: [PATCH 11/27] Add reference --- docs/design_docs/vector_reconstruction.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/design_docs/vector_reconstruction.md b/docs/design_docs/vector_reconstruction.md index 2a34710a6a..c6b17802e5 100644 --- a/docs/design_docs/vector_reconstruction.md +++ b/docs/design_docs/vector_reconstruction.md @@ -664,3 +664,9 @@ unit and convergence tests, but if the functionality becomes part of standard analysis or initialization workflows, a lightweight analysis regression test should be added later to guard against changes in reconstruction coefficients or basis transforms. + +## References + +Peixoto, P. S., and S. R. M. Barros (2014), On vector field reconstructions for +semi-Lagrangian transport methods on geodesic staggered grids, *Journal of +Computational Physics*, 273, 185-211, DOI:[10.1016/j.jcp.2014.04.043](https://doi.org/10.1016/j.jcp.2014.04.043). From 3e6a0e3740fb0582ef814a1a671a5075d4fe91ea Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 5 Apr 2026 09:50:02 +0200 Subject: [PATCH 12/27] Clarify vector reconstruction design scope and API Revise the vector reconstruction design doc to address review feedback. - move the proposed implementation from polaris/mpas to polaris/mesh - clarify that this PR is limited to Polaris design and mesh-field writing - state that MPAS-Ocean and Omega support will follow in later PRs - clarify that zonal/meridional outputs are tangential and exclude the radial part - explain that tangential reconstruction is produced in Cartesian form first - clarify why the in-memory cache and mesh-field export are separate steps - replace a few software-oriented phrases with plainer wording - tighten testing plans for culled meshes and nonzero radial contributions - narrow contributor/date updates to sections with substantive review-driven edits --- docs/design_docs/vector_reconstruction.md | 215 +++++++++++++--------- 1 file changed, 128 insertions(+), 87 deletions(-) diff --git a/docs/design_docs/vector_reconstruction.md b/docs/design_docs/vector_reconstruction.md index c6b17802e5..3df2fc3ba0 100644 --- a/docs/design_docs/vector_reconstruction.md +++ b/docs/design_docs/vector_reconstruction.md @@ -1,10 +1,11 @@ # Vector Reconstruction at MPAS Cell Centers -date: 2026/04/04 +date: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex ## Summary @@ -12,11 +13,14 @@ Contributors: This design proposes a new workflow for generating cell-center vector reconstruction coefficients for MPAS meshes in Python within `polaris`. The primary goal is for `polaris` to compute mesh-based reconstruction metadata and -coefficients, store them on an MPAS mesh, and enable Omega or MPAS-Ocean to -reconstruct vector fields at runtime in 3D Cartesian or zonal/meridional form. -Python is the intended home for coefficient generation because it provides the -most flexibility for experimenting with reconstruction methods, stencil choices, -and mesh-field layouts. +coefficients and write them to an MPAS mesh. A later MPAS-Ocean or Omega +implementation could read those mesh fields, but changes in those codes are out +of scope for this PR. Python is the intended home for coefficient generation +because it provides the most flexibility for experimenting with reconstruction +methods, stencil choices, and mesh-field layouts. The intent is that this +Python workflow will eventually become the authoritative way to generate these +coefficients, with follow-on PRs adding or updating support in MPAS-Ocean and +Omega to consume the stored mesh fields. The initial target field is `normalVelocity`, but the coefficient-generation workflow should be general enough to apply to any field represented as a scalar @@ -51,16 +55,16 @@ field-agnostic. The implementation will therefore separate: interface-centered radial component. Success means that `polaris` gains a reusable, documented, and testable -coefficient-generation workflow in `polaris/mpas` that: +coefficient-generation workflow in `polaris/mesh` that: - computes least-squares reconstruction coefficients and stencil metadata that can be stored on an MPAS mesh and consumed at runtime by Omega or MPAS-Ocean, -- treats Python in `polaris` as the authoritative implementation for generating - these coefficients, +- establishes Python in `polaris` as the implementation that is expected to + become the authoritative source for generating these coefficients, - reconstructs tangential vectors at MPAS cell centers with second-order behavior on spherical MPAS meshes, -- supports outputs in 3D Cartesian and zonal/meridional form at cell centers - and layer midpoints, and +- supports 3D Cartesian outputs and tangential zonal/meridional outputs at + cell centers and layer midpoints, and - is structured so later extension to arbitrary-point reconstruction can reuse the same geometry and reconstruction cache. @@ -91,18 +95,19 @@ that can be added to an MPAS mesh and used by Omega or MPAS-Ocean at runtime. ### Requirement: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex `polaris` shall support a generic interface that returns reconstructed vector -fields at cell centers and layer midpoints in: +fields at cell centers and layer midpoints as: - 3D Cartesian components, and -- zonal and meridional components. +- zonal and meridional components of the tangential part of the vector. If a radial interface field is supplied on `nVertLevelsP1`, the reconstructed vector shall include a radial contribution consistent with the MPAS local @@ -113,11 +118,12 @@ For velocity, the corresponding radial interface field is `vertVelocityTop`. ### Requirement: Reusable Precomputation and Extensibility -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex The implementation shall support precomputation of mesh-dependent reconstruction @@ -132,8 +138,10 @@ MPAS-Ocean. The implementation shall not rely on online coefficient generation in MPAS-Ocean or Omega. New reconstruction coefficients shall be generated in -Python in -`polaris` and then stored on the mesh for runtime use. +Python in `polaris` and then stored on the mesh for runtime use. Defining those +mesh fields is in scope for this design. Implementing MPAS-Ocean or Omega code +changes to read them is not, but follow-on PRs are expected to add that +support. ## Algorithm Design @@ -221,11 +229,12 @@ of the performance study rather than being fixed in the first implementation. ### Algorithm Design: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex The generic reconstruction produces only the tangential part of the vector. A @@ -256,10 +265,15 @@ $$ w_{c,k} = 0. $$ -The output should be available in both Cartesian and local geographic bases. -Using longitude $\lambda$ and latitude $\theta$ at the cell center, the -zonal-meridional decomposition should follow the same relationships used in the -MPAS framework: +The tangential reconstruction is first returned in 3D Cartesian components +$(u_x, u_y, u_z)$ even when the vector lies entirely in the tangent plane. If a +radial component is supplied, it is added in Cartesian form before any basis +change is applied. + +The output should also be available in the local geographic basis. Using +longitude $\lambda$ and latitude $\theta$ at the cell center, the zonal and +meridional components should be computed by projecting the Cartesian vector onto +the local eastward and northward directions used in the MPAS framework: $$ u_\mathrm{zonal} = -u_x \sin\lambda + u_y \cos\lambda, @@ -271,36 +285,39 @@ u_\mathrm{meridional} = $$ These formulas are appropriate both for tangential-only reconstruction and for -the full vector after radial composition. +the full vector after radial composition. Because the zonal and meridional +directions are tangent to the sphere, the radial contribution does not appear in +the zonal or meridional results. ### Algorithm Design: Reusable Precomputation and Extensibility -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis - Codex -The reconstruction should be divided into a mesh-dependent setup phase and a -field-dependent apply phase. The setup phase is the primary deliverable for -`polaris`, because it produces the coefficients and stencil metadata to be -stored on the mesh. The setup phase should be implemented in Python and should -be considered the authoritative implementation for the new reconstruction -method. The apply phase is still useful in Python for testing, analysis, and workflow -development, while Omega or MPAS-Ocean will consume the stored coefficients at -runtime rather than recomputing them. +The reconstruction should be divided into a mesh-dependent +coefficient-generation step and a field-dependent reconstruction step that uses +those precomputed coefficients. The coefficient-generation step is the primary +deliverable for `polaris` because it produces the coefficients and stencil +metadata to be stored on the mesh. This step should be implemented in Python +and is intended to become the authoritative implementation for the new +reconstruction method. The reconstruction step is also useful in Python for +testing and analysis. Changes to MPAS-Ocean or Omega to read the stored fields +would be follow-on work, not part of this PR. -The setup phase should precompute, for each cell: +The coefficient-generation step should precompute, for each cell: - the edge stencil, - projected edge geometry in the local tangent plane, - a weighted least-squares pseudo-inverse or equivalent coefficient matrix for the cell-center value. -The apply phase should then reduce to sparse local linear algebra. For a given -edge-normal field and vertical level, the tangential reconstruction at each cell -should be a weighted sum over the local stencil. +The reconstruction step should then reduce to weighted sums over each cell's +local stencil. For a given edge-normal field and vertical level, the tangential +reconstruction at each cell should be a weighted sum over the local stencil. This separation is important because: @@ -318,30 +335,32 @@ the field at vertices or other fixed locations and then interpolates. ### Implementation: Generic Tangential Reconstruction at Cell Centers -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex The proposed implementation adds a new coefficient-generation module to -`polaris/mpas`: +`polaris/mesh`: -- `polaris/mpas/reconstruct.py` +- `polaris/mesh/reconstruct.py` and a companion module for basis transforms and 3D composition: -- `polaris/mpas/vector.py` +- `polaris/mesh/vector.py` -The public API should be field-agnostic at its core. A likely structure is: +The public API should be generic enough to support velocity and other +edge-normal quantities. A likely structure is: ```python def build_reconstruction_cache(ds_mesh, method="lsq"): ... -def build_reconstruction_mesh_fields(ds_mesh, method="lsq"): +def build_reconstruction_mesh_fields(cache): ... @@ -364,27 +383,33 @@ def reconstruct_3d_cell_center(edge_normal_field, ds_mesh, - validate that the mesh is spherical, - determine the local tangent basis and local vertical direction at each cell, - build the local stencil and coefficient arrays, and -- return a lightweight cache object, likely as an `xarray.Dataset` or a small - dataclass containing `xarray.DataArray` members. +- return a compact collection of precomputed arrays, likely as an + `xarray.Dataset` or a small dataclass containing `xarray.DataArray` members. The default method should be all-cell least-squares. An optional `method` argument may later support `hybrid`, but that should not be the baseline path. -`build_reconstruction_mesh_fields()` should convert this cache into the mesh -variables needed by MPAS-Ocean or Omega at runtime. +`build_reconstruction_mesh_fields()` should convert that full in-memory cache +into the smaller set of MPAS mesh variables that belong in a mesh file. Keeping +these steps separate is useful because the cache may include projected +coordinates, basis vectors, and other helper arrays that are convenient in +Python but do not need to be written to the mesh. `reconstruct_tangential_cell_center()` should accept any edge-normal field with an `nEdges` dimension, preserving all non-horizontal dimensions. This makes the core implementation usable for velocity, pressure-gradient-like fields, or any -other tangential vector quantity represented through edge-normal components. +other tangential vector quantity represented through edge-normal components. It +should return a 3-component Cartesian vector that is tangent to the sphere at +each cell center, so `tangential_cartesian_to_zonal_meridional()` can be +applied directly to its output. `reconstruct_3d_cell_center()` should build on the tangential reconstruction by adding an optional radial interface field and returning Cartesian and -zonal/meridional components. This function should be field-agnostic. +zonal/meridional components. This function should remain usable for velocity and +other edge-normal quantities. -For xarray compatibility, the Python reconstruction implementation should -vectorize over all -non-horizontal dimensions and apply the precomputed weights only along `nEdges`. +For xarray compatibility, the Python reconstruction implementation should work +with extra dimensions and apply the precomputed weights only along `nEdges`. This will make the same code work naturally for: - `nEdges`, @@ -393,11 +418,12 @@ This will make the same code work naturally for: and similar combinations. -The least-squares path should avoid repeated dense solves at runtime. Instead, -the setup phase should precompute the operator that maps the stencil's edge -samples directly to the reconstructed Cartesian vector at the cell center. +The least-squares path should avoid repeated dense solves when reconstructing a +field. Instead, the coefficient-generation step should precompute the matrix +that maps the stencil's edge-normal values directly to the reconstructed +Cartesian vector at the cell center. -Pseudo code for setup: +Pseudo code for generating coeffs: ```python for cell in cells: @@ -406,7 +432,7 @@ for cell in cells: coeffs[cell] = build_lsq_center_coeffs(cell, stencil, basis) ``` -Pseudo code for apply: +Pseudo code for applying coeffs: ```python for cell in cells: @@ -416,15 +442,20 @@ for cell in cells: ) ``` -Because existing `polaris/mpas` modules are lightweight utilities, this design -fits the package style well. `polaris/mpas/__init__.py` should export the main -public functions. +`polaris/mesh/__init__.py` should export the main public functions. This design does not propose porting the new coefficient-generation algorithm to Fortran for online execution in MPAS-Ocean, and it does not propose computing these coefficients directly in Omega. Instead, the Python implementation in `polaris` should generate the coefficients once for a given mesh, after which runtime codes should only consume the stored mesh fields. +Follow-on PRs should update MPAS-Ocean and Omega to use those fields. + +#### Mesh Fields Written by Polaris + +This subsection defines the mesh fields that `polaris` will write. It includes +these details so the mesh-file interface is clear, but it does not propose +changes to MPAS-Ocean or Omega in this PR. The mesh-storage design needs special attention. Existing MPAS meshes already contain: @@ -461,14 +492,15 @@ current one-ring case and the new two-ring least-squares case. ### Implementation: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex -`polaris/mpas/vector.py` should contain the generic 3D-composition logic that +`polaris/mesh/vector.py` should contain the generic 3D-composition logic that sits on top of the tangential reconstruction and the coefficient-generation workflow. @@ -499,8 +531,8 @@ fields. A likely default is: - `vectorZonal`, - `vectorMeridional`. -To align with MPAS naming and user expectations, the design should also discuss -how this generic output maps onto velocity-oriented runtime fields such as: +If later work wants velocity-oriented names, the design should also note how +this generic output could map onto fields such as: - `uReconstructX`, - `uReconstructY`, @@ -508,28 +540,19 @@ how this generic output maps onto velocity-oriented runtime fields such as: - `uReconstructZonal`, - `uReconstructMeridional`. -The core implementation should remain field-agnostic even if some downstream -uses map the outputs onto velocity-specific names. - -At runtime, Omega or MPAS-Ocean should reconstruct: - -1. the tangential Cartesian vector from the stored stencil and coefficients, -2. the zonal and meridional components from the reconstructed Cartesian vector, -3. the full 3D vector by adding the radial contribution derived from an - optional radial interface field when needed. - -This keeps the coefficient storage generic while allowing field-specific -post-processing in the runtime model. Runtime codes should consume the stored -coefficients only; they should not be responsible for recomputing reconstruction -coefficients. +The core implementation should remain usable for velocity and for other +edge-normal quantities even if some later workflow maps the outputs onto +velocity-specific names. Choosing and implementing those downstream names is out +of scope for this PR. ### Implementation: Reusable Precomputation and Extensibility -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex The reconstruction cache should be serializable or at least representable as an @@ -551,17 +574,21 @@ need. A later arbitrary-point design can either: - add target-specific coefficient builders that reuse the same tangent-plane and stencil code. -Project-management-wise, the work can proceed in four steps: +The Polaris work described here can proceed in four steps: 1. implement the offline least-squares coefficient and stencil generator in `polaris`, -2. define the MPAS/Omega mesh-field representation for those coefficients and - write them to mesh files, +2. define the mesh-field representation written by `polaris` and write those + fields to mesh files, 3. implement basis transforms and the Python 3D reconstruction wrapper with an optional radial interface field, 4. document the API and add examples for coefficient generation and reconstruction from stored coefficients. +Any MPAS-Ocean or Omega code changes to read the new mesh fields would be +separate follow-on work after this Polaris design and implementation are in +place. + If performance studies later indicate that runtime application is too costly, a follow-on effort can evaluate the hybrid Perot/least-squares optimization. @@ -569,11 +596,12 @@ follow-on effort can evaluate the hybrid Perot/least-squares optimization. ### Testing and Validation: Generic Tangential Reconstruction at Cell Centers -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex Testing should include unit tests and convergence tests. @@ -583,7 +611,6 @@ Unit tests should verify: - constant tangential vector fields reconstruct with small error, - the reconstructed vector is tangent to the sphere when no radial component is supplied, -- the cache and direct-setup paths give identical results, - mesh-output fields contain the expected stencil width and coefficient layout. Convergence tests should use analytic spherical flows already available in @@ -603,6 +630,10 @@ For the coefficient-generation path, tests should also confirm that applying the stored mesh coefficients reproduces the same result as the direct Python reconstruction. +At least one test case should use a culled ocean mesh to confirm that stencil +construction and reconstruction behave sensibly near topography and other +non-global boundaries. + Performance tests should compare: - coefficient-generation cost in Python for representative mesh sizes, @@ -614,11 +645,12 @@ evaluation of the hybrid optimization path. ### Testing and Validation: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: - Xylar Asay-Davis +- Carolyn Begeman - Codex 3D-vector tests should verify: @@ -626,6 +658,10 @@ Contributors: - Cartesian and zonal/meridional outputs are mutually consistent, - omission of the radial interface field yields zero radial contribution, - midpoint radial velocity is computed correctly from interface values, +- a case with nonzero radial contribution is reconstructed correctly in + Cartesian form, and +- the zonal and meridional outputs are unchanged by the added radial + contribution, - shape and dimension handling are correct for `nVertLevels`, `nVertLevelsP1`, and optional `Time`. @@ -642,7 +678,7 @@ the motivating example. ### Testing and Validation: Reusable Precomputation and Extensibility -Date last modified: 2026/04/04 +Date last modified: 2026/04/05 Contributors: @@ -650,8 +686,13 @@ Contributors: - Codex For the primary workflow, tests should also verify that mesh-writing preserves -the coefficient fields exactly and that a runtime-style application of those -fields gives identical answers before and after I/O. +the coefficient fields exactly and that applying those stored fields gives +identical answers before and after I/O. + +Tests in this section should also verify that converting the in-memory cache to +mesh fields preserves the coefficients needed for reconstruction, and that +applying those stored fields reproduces the same result as applying the +in-memory cache directly. Documentation tests or example scripts should also be added to demonstrate: From a65f7eb0071adf5879406ab8e6ebe8394ad5cd66 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 6 Apr 2026 17:23:20 +0200 Subject: [PATCH 13/27] Add radial component to the local geographic reconstruction --- docs/design_docs/vector_reconstruction.md | 76 ++++++++++++++--------- 1 file changed, 47 insertions(+), 29 deletions(-) diff --git a/docs/design_docs/vector_reconstruction.md b/docs/design_docs/vector_reconstruction.md index 3df2fc3ba0..4d868f7e54 100644 --- a/docs/design_docs/vector_reconstruction.md +++ b/docs/design_docs/vector_reconstruction.md @@ -1,6 +1,6 @@ # Vector Reconstruction at MPAS Cell Centers -date: 2026/04/05 +date: 2026/04/06 Contributors: @@ -50,7 +50,7 @@ Although velocity is the motivating use case, the core reconstruction should be field-agnostic. The implementation will therefore separate: - a generic reconstruction of tangential vectors from edge-normal components, -- basis transformations between Cartesian and zonal/meridional coordinates, and +- basis transformations between Cartesian and local geographic coordinates, and - generic 3D composition from a tangential component and an optional interface-centered radial component. @@ -63,7 +63,7 @@ coefficient-generation workflow in `polaris/mesh` that: become the authoritative source for generating these coefficients, - reconstructs tangential vectors at MPAS cell centers with second-order behavior on spherical MPAS meshes, -- supports 3D Cartesian outputs and tangential zonal/meridional outputs at +- supports 3D Cartesian outputs and local zonal/meridional/radial outputs at cell centers and layer midpoints, and - is structured so later extension to arbitrary-point reconstruction can reuse the same geometry and reconstruction cache. @@ -95,7 +95,7 @@ that can be added to an MPAS mesh and used by Omega or MPAS-Ocean at runtime. ### Requirement: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/05 +Date last modified: 2026/04/06 Contributors: @@ -107,13 +107,16 @@ Contributors: fields at cell centers and layer midpoints as: - 3D Cartesian components, and -- zonal and meridional components of the tangential part of the vector. +- zonal, meridional, and radial components in the local geographic basis. If a radial interface field is supplied on `nVertLevelsP1`, the reconstructed vector shall include a radial contribution consistent with the MPAS local vertical direction. If no radial interface field is supplied, the radial contribution shall default to zero. +For a horizontal-only reconstruction, the local geographic radial output shall +therefore be zero and may be ignored by downstream users. + For velocity, the corresponding radial interface field is `vertVelocityTop`. ### Requirement: Reusable Precomputation and Extensibility @@ -229,7 +232,7 @@ of the performance study rather than being fixed in the first implementation. ### Algorithm Design: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/05 +Date last modified: 2026/04/06 Contributors: @@ -270,10 +273,11 @@ $(u_x, u_y, u_z)$ even when the vector lies entirely in the tangent plane. If a radial component is supplied, it is added in Cartesian form before any basis change is applied. -The output should also be available in the local geographic basis. Using -longitude $\lambda$ and latitude $\theta$ at the cell center, the zonal and -meridional components should be computed by projecting the Cartesian vector onto -the local eastward and northward directions used in the MPAS framework: +The output should also be available in the local geographic basis as zonal, +meridional, and radial components. Using longitude $\lambda$ and latitude +$\theta$ at the cell center, these components should be computed by projecting +the Cartesian vector onto the local eastward, northward, and vertical +directions used in the MPAS framework: $$ u_\mathrm{zonal} = -u_x \sin\lambda + u_y \cos\lambda, @@ -284,10 +288,16 @@ u_\mathrm{meridional} = -(u_x \cos\lambda + u_y \sin\lambda)\sin\theta + u_z \cos\theta. $$ +$$ +u_\mathrm{radial} = +(u_x \cos\lambda + u_y \sin\lambda)\cos\theta + u_z \sin\theta. +$$ + These formulas are appropriate both for tangential-only reconstruction and for -the full vector after radial composition. Because the zonal and meridional -directions are tangent to the sphere, the radial contribution does not appear in -the zonal or meridional results. +the full vector after radial composition. For a tangential-only reconstruction, +$u_\mathrm{radial}$ should be zero up to truncation error. When a radial +contribution is included, it should appear explicitly in +$u_\mathrm{radial}$. ### Algorithm Design: Reusable Precomputation and Extensibility @@ -335,7 +345,7 @@ the field at vertices or other fixed locations and then interpolates. ### Implementation: Generic Tangential Reconstruction at Cell Centers -Date last modified: 2026/04/05 +Date last modified: 2026/04/06 Contributors: @@ -369,7 +379,7 @@ def reconstruct_tangential_cell_center(edge_normal_field, ds_mesh, ... -def tangential_cartesian_to_zonal_meridional(u_x, u_y, u_z, ds_mesh): +def cartesian_to_local_geographic(u_x, u_y, u_z, ds_mesh): ... @@ -400,12 +410,12 @@ an `nEdges` dimension, preserving all non-horizontal dimensions. This makes the core implementation usable for velocity, pressure-gradient-like fields, or any other tangential vector quantity represented through edge-normal components. It should return a 3-component Cartesian vector that is tangent to the sphere at -each cell center, so `tangential_cartesian_to_zonal_meridional()` can be -applied directly to its output. +each cell center, so `cartesian_to_local_geographic()` can be applied directly +to its output. `reconstruct_3d_cell_center()` should build on the tangential reconstruction by adding an optional radial interface field and returning Cartesian and -zonal/meridional components. This function should remain usable for velocity and +local geographic components. This function should remain usable for velocity and other edge-normal quantities. For xarray compatibility, the Python reconstruction implementation should work @@ -492,7 +502,7 @@ current one-ring case and the new two-ring least-squares case. ### Implementation: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/05 +Date last modified: 2026/04/06 Contributors: @@ -510,7 +520,8 @@ The generic 3D wrapper should: 2. broadcast the tangential result to all non-horizontal dimensions, 3. derive midpoint radial values from an optional radial interface field, 4. add the radial contribution in the local vertical direction, and -5. derive zonal and meridional components from the resulting Cartesian vector. +5. derive zonal, meridional, and radial components from the resulting Cartesian + vector. For the optional radial interface field, the implementation should: @@ -529,7 +540,8 @@ fields. A likely default is: - `vectorY`, - `vectorZ`, - `vectorZonal`, -- `vectorMeridional`. +- `vectorMeridional`, +- `vectorRadial`. If later work wants velocity-oriented names, the design should also note how this generic output could map onto fields such as: @@ -538,7 +550,8 @@ this generic output could map onto fields such as: - `uReconstructY`, - `uReconstructZ`, - `uReconstructZonal`, -- `uReconstructMeridional`. +- `uReconstructMeridional`, +- `uReconstructRadial`. The core implementation should remain usable for velocity and for other edge-normal quantities even if some later workflow maps the outputs onto @@ -645,7 +658,7 @@ evaluation of the hybrid optimization path. ### Testing and Validation: 3D Vector Outputs at Layer Midpoints -Date last modified: 2026/04/05 +Date last modified: 2026/04/06 Contributors: @@ -655,13 +668,17 @@ Contributors: 3D-vector tests should verify: -- Cartesian and zonal/meridional outputs are mutually consistent, +- Cartesian and local geographic outputs are mutually consistent, - omission of the radial interface field yields zero radial contribution, - midpoint radial velocity is computed correctly from interface values, +- the reconstructed local geographic radial component equals the simple mean of + the adjacent top and bottom interface values, - a case with nonzero radial contribution is reconstructed correctly in Cartesian form, and -- the zonal and meridional outputs are unchanged by the added radial - contribution, +- the local geographic radial output matches the expected midpoint radial + velocity, and +- the zonal and meridional outputs remain consistent with the horizontal part + of the reconstructed vector, - shape and dimension handling are correct for `nVertLevels`, `nVertLevelsP1`, and optional `Time`. @@ -672,9 +689,10 @@ A simple analytic test can define: - an expected full Cartesian vector at layer midpoints. The generic wrapper should reproduce the known Cartesian and zonal/meridional -fields to within expected truncation error. A velocity-specific test case using -`normalVelocity` and `vertVelocityTop` should remain part of the test suite as -the motivating example. +fields, together with the expected radial field, to within expected truncation +error. A velocity-specific test case using `normalVelocity` and +`vertVelocityTop` should remain part of the test suite as the motivating +example. ### Testing and Validation: Reusable Precomputation and Extensibility From bcdb96f3406c72f2ac9dd96749ec2542cf3c6480 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 5 Apr 2026 15:46:14 +0200 Subject: [PATCH 14/27] Add AI (Codex and Copilot) instructions These are designed to keep these tools from straying too far from the Polaris coding style and to make the easier for developers to use without later fighting with pre-commit. --- .github/copilot-instructions.md | 10 ++++++ .github/instructions/python.instructions.md | 15 ++++++++ AGENTS.md | 32 +++++++++++++++++ docs/developers_guide/quick_start.md | 38 +++++++++++++++++++++ 4 files changed, 95 insertions(+) create mode 100644 .github/copilot-instructions.md create mode 100644 .github/instructions/python.instructions.md create mode 100644 AGENTS.md diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md new file mode 100644 index 0000000000..c11b3a833d --- /dev/null +++ b/.github/copilot-instructions.md @@ -0,0 +1,10 @@ +# Polaris Copilot Instructions + +Follow the repository's automated style configuration in +`pyproject.toml` and `.pre-commit-config.yaml`. + +- Keep changes consistent with existing Polaris patterns. +- For Python, follow the path-specific instructions in + `.github/instructions/python.instructions.md`. +- Prefer changes that pass the configured pre-commit hooks without + adding ignores or suppressions. diff --git a/.github/instructions/python.instructions.md b/.github/instructions/python.instructions.md new file mode 100644 index 0000000000..43450878cf --- /dev/null +++ b/.github/instructions/python.instructions.md @@ -0,0 +1,15 @@ +--- +applyTo: "**/*.py" +--- + +# Python Instructions + +- Keep lines at 79 characters or fewer whenever possible. +- Adhere to `ruff format` formatting. +- Keep imports at module scope whenever possible. Avoid local imports + unless they are needed for circular-import avoidance, lazy loading, or + optional dependencies. +- Avoid nested functions whenever possible. +- Prefer public functions before private helper functions whenever + practical. +- Prefer private module-level helper functions over nested helpers. diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000000..967ed1676c --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,32 @@ +# Polaris Agent Instructions + +These instructions apply to the whole repository unless a deeper +`AGENTS.md` overrides them. + +## Source of truth + +- Follow the repo's automated style and lint configuration in + `pyproject.toml` and `.pre-commit-config.yaml`. +- If an instruction here conflicts with automated tooling, follow the + automated tooling. + +## Python style + +- Keep Python lines at 79 characters or fewer whenever possible. +- Use `ruff format` style. Do not preserve manual formatting that Ruff + would rewrite. +- Keep imports at module scope whenever possible. Avoid local imports + unless they are needed to prevent circular imports, defer expensive + dependencies, or avoid optional dependency failures. +- Avoid nested functions whenever possible. Prefer private module-level + helpers instead. +- Put public functions before private helper functions whenever + practical. +- Name private helper functions with a leading underscore when that fits + existing repo conventions. + +## Validation + +- Run relevant pre-commit hooks on changed files before finishing when + practical. +- Prefer fixing lint and formatting issues rather than suppressing them. diff --git a/docs/developers_guide/quick_start.md b/docs/developers_guide/quick_start.md index 28055a9961..84077ca1d0 100644 --- a/docs/developers_guide/quick_start.md +++ b/docs/developers_guide/quick_start.md @@ -575,6 +575,44 @@ for some tips on checking code style in VS Code. Once you open a pull request for your feature, there is an additional PEP8 style check at this stage (again using pre-commit). +(dev-polaris-ai-instructions)= + +## AI coding assistant instructions + +Polaris also includes repository-level instructions for AI coding assistants: + +- `AGENTS.md` provides instructions for Codex and other agent-style tools that + support repository guidance. +- `.github/copilot-instructions.md` provides repository-wide instructions for + GitHub Copilot. +- `.github/instructions/python.instructions.md` provides Python-specific + Copilot instructions. + +These files are intended to capture style preferences that are useful for +human contributors and AI tools but are not always enforceable automatically. +Examples include preferring module-level helper functions over nested +functions, putting public functions before private helpers when practical, and +avoiding local imports unless they are needed. + +To keep maintenance manageable: + +- Put anything that can be enforced automatically in `pyproject.toml` and + `.pre-commit-config.yaml` first, then keep the AI instruction files focused + on higher-level preferences. +- Keep the shared guidance short and stable so it is practical to mirror in + both Codex and Copilot instruction files. +- Treat `AGENTS.md` as the easiest place to edit shared prose first, then + update the two Copilot instruction files to match when shared guidance + changes. +- Use `.github/instructions/*.instructions.md` only for path-specific guidance + such as Python conventions, rather than repeating general repository advice + in many places. +- If an AI instruction conflicts with automated tooling, the automated tooling + should win. + +When you add or change a style preference, update these three files in the same +pull request so they stay aligned. + (dev-polaris-repo-advanced)= ## Set up a polaris repository with worktrees: for advanced users From 8f5d6c897c4ad5960e93f61911837a285d0d0bb8 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 6 Apr 2026 17:18:41 +0200 Subject: [PATCH 15/27] Add AI instructions to follow templates for tasks added to docs. --- .github/copilot-instructions.md | 2 ++ .github/instructions/docs.instructions.md | 13 +++++++++++++ AGENTS.md | 8 ++++++++ .../dev_add_category_of_tasks/documenting.md | 5 +++-- 4 files changed, 26 insertions(+), 2 deletions(-) create mode 100644 .github/instructions/docs.instructions.md diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index c11b3a833d..3efc72d3ca 100644 --- a/.github/copilot-instructions.md +++ b/.github/copilot-instructions.md @@ -6,5 +6,7 @@ Follow the repository's automated style configuration in - Keep changes consistent with existing Polaris patterns. - For Python, follow the path-specific instructions in `.github/instructions/python.instructions.md`. +- For documentation in `docs/`, follow the path-specific instructions in + `.github/instructions/docs.instructions.md`. - Prefer changes that pass the configured pre-commit hooks without adding ignores or suppressions. diff --git a/.github/instructions/docs.instructions.md b/.github/instructions/docs.instructions.md new file mode 100644 index 0000000000..676b95a142 --- /dev/null +++ b/.github/instructions/docs.instructions.md @@ -0,0 +1,13 @@ +--- +applyTo: "docs/**/*.md" +--- + +# Documentation Instructions + +- When writing documentation for component tasks, follow the relevant + `template.md` format and inline instructions whenever a component + task template is available. +- Prefer starting from the existing template instead of writing + component task documentation pages from scratch. +- Keep new component task documentation consistent with neighboring + pages in the same component. diff --git a/AGENTS.md b/AGENTS.md index 967ed1676c..8cf300d3e0 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -25,6 +25,14 @@ These instructions apply to the whole repository unless a deeper - Name private helper functions with a leading underscore when that fits existing repo conventions. +## Documentation + +- When writing documentation for component tasks, follow the relevant + `template.md` format and its inline instructions whenever a component + task template is available. +- Prefer starting from the existing template instead of creating task + documentation pages from scratch. + ## Validation - Run relevant pre-commit hooks on changed files before finishing when diff --git a/docs/tutorials/dev_add_category_of_tasks/documenting.md b/docs/tutorials/dev_add_category_of_tasks/documenting.md index 7da444bafd..798396fee9 100644 --- a/docs/tutorials/dev_add_category_of_tasks/documenting.md +++ b/docs/tutorials/dev_add_category_of_tasks/documenting.md @@ -16,8 +16,9 @@ describing the category of tasks and its tasks and steps. For the user's guide, make a copy of `docs/users_guide//tasks/template.md` called `docs/users_guide//tasks/.md`. In that file, you -should describe the category of tasks and its tasks in a way that would be -relevant for a user wanting to run the task and look at the output. +should follow the template's format and inline instructions, then describe the +category of tasks and its tasks in a way that would be relevant for a user +wanting to run the task and look at the output. This file should describe all of the config options relevant the tasks collectively and each task (if it has its own config options), including what they are used for and whether it is a good idea to modify them. Add From 15672b7eb0cacdd01e1a6f9362063d527ce5edf4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 14:19:03 -0700 Subject: [PATCH 16/27] Add replacements argument for convergence forward --- polaris/ocean/convergence/forward.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py index 1f7643c6a1..16d8cbbc8d 100644 --- a/polaris/ocean/convergence/forward.py +++ b/polaris/ocean/convergence/forward.py @@ -36,6 +36,7 @@ def __init__( yaml_filename='forward.yaml', mesh_input_filename='mesh.nc', options=None, + replacements=None, graph_target=None, output_filename='output.nc', validate_vars=None, @@ -98,6 +99,7 @@ def __init__( graph_target=graph_target, ) + self.replacements = replacements self.refinement = refinement self.refinement_factor = refinement_factor self.package = package @@ -202,6 +204,8 @@ def dynamic_model_config(self, at_setup): output_interval=output_interval_str, output_freq=f'{output_freq}', ) + if self.replacements is not None: + replacements.update(self.replacements) self.add_yaml_file( self.package, From 22c25784ddf347007a6fc538fdeeb1573a3ea0f9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 15:00:46 -0700 Subject: [PATCH 17/27] Set flow type/id using yaml replacements; set PrescribeVelocityType to init for rotation_2d --- polaris/tasks/ocean/sphere_transport/forward.py | 15 ++++++++++----- polaris/tasks/ocean/sphere_transport/forward.yaml | 4 ++++ 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/forward.py b/polaris/tasks/ocean/sphere_transport/forward.py index 74f852593c..6fd241a27c 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.py +++ b/polaris/tasks/ocean/sphere_transport/forward.py @@ -56,10 +56,15 @@ def __init__( 'divergent_2d': 3, 'correlated_tracers_2d': 4, } - namelist_options = { - 'mpas-ocean': { - 'config_transport_tests_flow_id': flow_id[case_name] - } + flow_type = { + 'rotation_2d': 'Init', + 'nondivergent_2d': 'None', + 'divergent_2d': 'None', + 'correlated_tracers_2d': 'None', + } + replacements = { + 'flow_id': flow_id[case_name], + 'flow_type': flow_type[case_name], } validate_vars = ['normalVelocity', 'tracer1', 'tracer2', 'tracer3'] super().__init__( @@ -74,7 +79,7 @@ def __init__( output_filename='output.nc', validate_vars=validate_vars, check_properties=['mass conservation', 'tracer conservation'], - options=namelist_options, + replacements=replacements, graph_target=f'{base_mesh.path}/graph.info', refinement_factor=refinement_factor, refinement=refinement, diff --git a/polaris/tasks/ocean/sphere_transport/forward.yaml b/polaris/tasks/ocean/sphere_transport/forward.yaml index 7ef8c50369..a067ff5ace 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.yaml +++ b/polaris/tasks/ocean/sphere_transport/forward.yaml @@ -13,6 +13,8 @@ mpas-ocean: config_block_decomp_file_prefix: graph.info.part. advection: config_vert_coord_movement: impermeable_interfaces + transport_tests: + config_transport_tests_flow_id: {{ flow_id }} debug: config_disable_thick_sflux: true config_disable_vel_all_tend: true @@ -61,6 +63,8 @@ mpas-ocean: - relativeVorticityCell Omega: + State: + PrescribeVelocityType: {{ flow_type }} Tendencies: ThicknessFluxTendencyEnable : true PVTendencyEnable : false From fe6963282e8abe534f7000b0091f5068a617c43b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 16:47:16 -0700 Subject: [PATCH 18/27] Support non-divergent and divergent flow types --- polaris/tasks/ocean/sphere_transport/forward.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/forward.py b/polaris/tasks/ocean/sphere_transport/forward.py index 6fd241a27c..1af29ef37c 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.py +++ b/polaris/tasks/ocean/sphere_transport/forward.py @@ -58,9 +58,9 @@ def __init__( } flow_type = { 'rotation_2d': 'Init', - 'nondivergent_2d': 'None', - 'divergent_2d': 'None', - 'correlated_tracers_2d': 'None', + 'nondivergent_2d': 'NonDivergent', + 'divergent_2d': 'Divergent', + 'correlated_tracers_2d': 'NonDivergent', } replacements = { 'flow_id': flow_id[case_name], From 5c94f09c8cee5a541c9a10a23755d286de0fbdfa Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 17 Mar 2026 10:57:36 -0700 Subject: [PATCH 19/27] Use ds_mesh.sphere_radius for non-div, div --- polaris/tasks/ocean/sphere_transport/init.py | 4 ++-- .../sphere_transport/resources/flow_types.py | 18 ++++++++++-------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/init.py b/polaris/tasks/ocean/sphere_transport/init.py index 767e5eef4f..9751dcca67 100644 --- a/polaris/tasks/ocean/sphere_transport/init.py +++ b/polaris/tasks/ocean/sphere_transport/init.py @@ -173,7 +173,7 @@ def run(self): section = config[case_name] vel_amp = section.getfloat('vel_amp') u, v = flow_divergent( - 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour + 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour, sphere_radius ) elif ( case_name == 'nondivergent_2d' @@ -182,7 +182,7 @@ def run(self): section = config[case_name] vel_amp = section.getfloat('vel_amp') u, v = flow_nondivergent( - 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour + 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour, sphere_radius ) else: raise ValueError(f'Unexpected test case name {case_name}') diff --git a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py index bf973e69de..b1efaf1de3 100644 --- a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py +++ b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py @@ -8,7 +8,7 @@ ) -def flow_nondivergent(t, lon, lat, u_0, tau): +def flow_nondivergent(t, lon, lat, u_0, tau, sphere_radius): """ Compute a nondivergent velocity field @@ -40,14 +40,14 @@ def flow_nondivergent(t, lon, lat, u_0, tau): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (1 / tau) * ( - u_0 * (sin(lon_p) ** 2) * sin(2 * lat) * cost + 2.0 * pi * coslat + u = (u_0 * sphere_radius / tau) * ( + (sin(lon_p) ** 2) * sin(2 * lat) * cost + 2.0 * pi * coslat ) - v = (u_0 / tau) * sin(2 * lon_p) * coslat * cost + v = (u_0 * sphere_radius / tau) * sin(2 * lon_p) * coslat * cost return u, v -def flow_divergent(t, lon, lat, u_0, tau): +def flow_divergent(t, lon, lat, u_0, tau, sphere_radius): """ Compute a divergent velocity field @@ -79,11 +79,13 @@ def flow_divergent(t, lon, lat, u_0, tau): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (1 / tau) * ( - -u_0 * (sin(lon_p / 2) ** 2) * sin(2 * lat) * (coslat**2) * cost + u = (u_0 * sphere_radius / tau) * ( + -(sin(lon_p / 2) ** 2) * sin(2 * lat) * (coslat**2) * cost + 2.0 * pi * coslat ) - v = (u_0 / (2 * tau)) * sin(lon_p) * (coslat**3) * cost + v = (u_0 * sphere_radius / (2 * tau)) * ( + sin(lon_p) * (coslat**3) * cost + ) return u, v From 9d8b1f78b013eda1708c414e41b611a8197f58ef Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 17 Mar 2026 11:03:10 -0700 Subject: [PATCH 20/27] Use PCD constant earth_radius for sphere_transport --- polaris/tasks/ocean/sphere_transport/init.py | 18 +++++++------- .../sphere_transport/resources/flow_types.py | 24 +++++++++---------- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/init.py b/polaris/tasks/ocean/sphere_transport/init.py index 9751dcca67..16527553bc 100644 --- a/polaris/tasks/ocean/sphere_transport/init.py +++ b/polaris/tasks/ocean/sphere_transport/init.py @@ -1,6 +1,7 @@ import numpy as np import xarray as xr +from polaris.constants import get_constant from polaris.ocean.model import OceanIOStep from polaris.ocean.vertical import init_vertical_coord from polaris.tasks.ocean.sphere_transport.resources.flow_types import ( @@ -16,6 +17,8 @@ xyztrig, ) +earth_radius = get_constant('mean_radius') + class Init(OceanIOStep): """ @@ -80,7 +83,6 @@ def run(self): latEdge = ds_mesh.latEdge lonCell = ds_mesh.lonCell lonEdge = ds_mesh.lonEdge - sphere_radius = ds_mesh.sphere_radius ds = ds_mesh.copy() @@ -95,7 +97,7 @@ def run(self): ds['salinity'] = salinity * xr.ones_like(ds.temperature) # tracer1 - tracer1 = xyztrig(lonCell, latCell, sphere_radius) + tracer1 = xyztrig(lonCell, latCell, earth_radius) # tracer2 section = config['sphere_transport'] @@ -108,7 +110,7 @@ def run(self): radius, background_value, amplitude, - sphere_radius, + earth_radius, ) # tracer3 @@ -128,7 +130,7 @@ def run(self): radius, background_value, amplitude, - sphere_radius, + earth_radius, ) _, tracer1_array = np.meshgrid(ds.refZMid.values, tracer1) _, tracer2_array = np.meshgrid(ds.refZMid.values, tracer2) @@ -166,14 +168,12 @@ def run(self): case_name, 'rotation_vector', dtype=float ) vector = np.array(rotation_vector) - u, v = flow_rotation( - lonEdge, latEdge, vector, vel_pd * s_per_hour, sphere_radius - ) + u, v = flow_rotation(lonEdge, latEdge, vector, vel_pd * s_per_hour) elif case_name == 'divergent_2d': section = config[case_name] vel_amp = section.getfloat('vel_amp') u, v = flow_divergent( - 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour, sphere_radius + 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour ) elif ( case_name == 'nondivergent_2d' @@ -182,7 +182,7 @@ def run(self): section = config[case_name] vel_amp = section.getfloat('vel_amp') u, v = flow_nondivergent( - 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour, sphere_radius + 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour ) else: raise ValueError(f'Unexpected test case name {case_name}') diff --git a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py index b1efaf1de3..4be555fcaa 100644 --- a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py +++ b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py @@ -3,12 +3,15 @@ from mpas_tools.transects import lon_lat_to_cartesian from numpy import cos, pi, sin +from polaris.constants import get_constant from polaris.mesh.spherical import ( calc_vector_east_north, ) +earth_radius = get_constant('mean_radius') -def flow_nondivergent(t, lon, lat, u_0, tau, sphere_radius): + +def flow_nondivergent(t, lon, lat, u_0, tau): """ Compute a nondivergent velocity field @@ -40,14 +43,14 @@ def flow_nondivergent(t, lon, lat, u_0, tau, sphere_radius): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (u_0 * sphere_radius / tau) * ( + u = (u_0 * earth_radius / tau) * ( (sin(lon_p) ** 2) * sin(2 * lat) * cost + 2.0 * pi * coslat ) - v = (u_0 * sphere_radius / tau) * sin(2 * lon_p) * coslat * cost + v = (u_0 * earth_radius / tau) * sin(2 * lon_p) * coslat * cost return u, v -def flow_divergent(t, lon, lat, u_0, tau, sphere_radius): +def flow_divergent(t, lon, lat, u_0, tau): """ Compute a divergent velocity field @@ -79,17 +82,15 @@ def flow_divergent(t, lon, lat, u_0, tau, sphere_radius): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (u_0 * sphere_radius / tau) * ( + u = (u_0 * earth_radius / tau) * ( -(sin(lon_p / 2) ** 2) * sin(2 * lat) * (coslat**2) * cost + 2.0 * pi * coslat ) - v = (u_0 * sphere_radius / (2 * tau)) * ( - sin(lon_p) * (coslat**3) * cost - ) + v = (u_0 * earth_radius / (2 * tau)) * (sin(lon_p) * (coslat**3) * cost) return u, v -def flow_rotation(lon, lat, omega, tau, sphere_radius): +def flow_rotation(lon, lat, omega, tau): """ Compute a rotational velocity field @@ -108,9 +109,6 @@ def flow_rotation(lon, lat, omega, tau, sphere_radius): tau : float time in seconds for the flow to circumnavigate the sphere - sphere_radius : float - radius of the sphere - Returns ------- u : np.ndarray of type float @@ -120,7 +118,7 @@ def flow_rotation(lon, lat, omega, tau, sphere_radius): meridional velocity """ omega = (2.0 * pi / tau) * (omega / np.linalg.norm(omega)) - x, y, z = lon_lat_to_cartesian(lon, lat, sphere_radius, degrees=False) + x, y, z = lon_lat_to_cartesian(lon, lat, earth_radius, degrees=False) xyz = np.stack((x, y, z), axis=1) vel = np.cross(omega, np.transpose(xyz), axis=0) east, north = calc_vector_east_north(x, y, z) From 2265d1ae47a926dde2a23df4b26dfcde36325ed2 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 17 Mar 2026 11:22:19 -0700 Subject: [PATCH 21/27] Use config option mapping to ensure the same hadv order is used for sphere_transport tests --- polaris/ocean/model/mpaso_to_omega.yaml | 1 + polaris/tasks/ocean/sphere_transport/forward.yaml | 2 ++ 2 files changed, 3 insertions(+) diff --git a/polaris/ocean/model/mpaso_to_omega.yaml b/polaris/ocean/model/mpaso_to_omega.yaml index 74e75e4125..4b5f819b7e 100644 --- a/polaris/ocean/model/mpaso_to_omega.yaml +++ b/polaris/ocean/model/mpaso_to_omega.yaml @@ -76,6 +76,7 @@ config: advection: Advection options: config_thickness_flux_type: FluxThicknessType + config_horiz_tracer_adv_order: HorzTracerFluxOrder - section: advection: Tendencies diff --git a/polaris/tasks/ocean/sphere_transport/forward.yaml b/polaris/tasks/ocean/sphere_transport/forward.yaml index a067ff5ace..016a17afed 100644 --- a/polaris/tasks/ocean/sphere_transport/forward.yaml +++ b/polaris/tasks/ocean/sphere_transport/forward.yaml @@ -5,6 +5,8 @@ ocean: time_integration: config_dt: {{ dt }} config_time_integrator: {{ time_integrator }} + advection: + config_horiz_tracer_adv_order: 3 mpas-ocean: run_modes: From 1b8343c8e266c1b1e45ffa69dc2f6ce541bc0127 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 7 Apr 2026 06:49:01 -0700 Subject: [PATCH 22/27] Update E3SM submodule --- e3sm_submodules/Omega | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/e3sm_submodules/Omega b/e3sm_submodules/Omega index c2542a22e0..589689d3e9 160000 --- a/e3sm_submodules/Omega +++ b/e3sm_submodules/Omega @@ -1 +1 @@ -Subproject commit c2542a22e05b910b8688e10ed5df0521d9ddc4d6 +Subproject commit 589689d3e9dbff201e40fe7c11d7297aeafc63b1 From 730ac463d42501485e91701ae6b6f7e9561f81b5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 7 Apr 2026 14:54:37 -0700 Subject: [PATCH 23/27] Remove cpptrace as omega build dep --- polaris/build/build_omega.template | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/polaris/build/build_omega.template b/polaris/build/build_omega.template index bff15a5cdf..7e2023e514 100644 --- a/polaris/build/build_omega.template +++ b/polaris/build/build_omega.template @@ -26,8 +26,7 @@ git submodule update --init --recursive \ externals/YAKL \ externals/ekat \ externals/scorpio \ - externals/cpptrace \ - components/omega/external \ + components/omega/external/cpptrace \ cime cd ${cwd} From cadfcce064a0e5cf256fdf7fb5c794c646d710dd Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 7 Apr 2026 14:56:18 -0700 Subject: [PATCH 24/27] Include velocityZonal, Merid in initial condition for debugging --- polaris/tasks/ocean/sphere_transport/init.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/init.py b/polaris/tasks/ocean/sphere_transport/init.py index 16527553bc..1a07c7ebe8 100644 --- a/polaris/tasks/ocean/sphere_transport/init.py +++ b/polaris/tasks/ocean/sphere_transport/init.py @@ -169,12 +169,16 @@ def run(self): ) vector = np.array(rotation_vector) u, v = flow_rotation(lonEdge, latEdge, vector, vel_pd * s_per_hour) + u_cell, v_cell = flow_rotation(lonCell, latCell, vector, vel_pd * s_per_hour) elif case_name == 'divergent_2d': section = config[case_name] vel_amp = section.getfloat('vel_amp') u, v = flow_divergent( 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour ) + u_cell, v_cell = flow_divergent( + 0.0, lonCell, latCell, vel_amp, vel_pd * s_per_hour + ) elif ( case_name == 'nondivergent_2d' or case_name == 'correlated_tracers_2d' @@ -184,6 +188,9 @@ def run(self): u, v = flow_nondivergent( 0.0, lonEdge, latEdge, vel_amp, vel_pd * s_per_hour ) + u_cell, v_cell = flow_nondivergent( + 0.0, lonCell, latCell, vel_amp, vel_pd * s_per_hour + ) else: raise ValueError(f'Unexpected test case name {case_name}') @@ -192,8 +199,10 @@ def run(self): ) normalVelocity, _ = xr.broadcast(normalVelocity, ds.refZMid) ds['normalVelocity'] = normalVelocity.expand_dims(dim='Time', axis=0) - ds['velocityZonal'] = np.nan * xr.zeros_like(ds.temperature) - ds['velocityMeridional'] = np.nan * xr.zeros_like(ds.temperature) + u_cell_broadcast, _ = xr.broadcast(xr.DataArray(u_cell, dims='nCells'), ds.refZMid) + v_cell_broadcast, _ = xr.broadcast(xr.DataArray(v_cell, dims='nCells'), ds.refZMid) + ds['velocityZonal'] = u_cell_broadcast.expand_dims(dim='Time', axis=0) + ds['velocityMeridional'] = v_cell_broadcast.expand_dims(dim='Time', axis=0) ds['fCell'] = xr.zeros_like(ds_mesh.xCell) ds['fEdge'] = xr.zeros_like(ds_mesh.xEdge) From 2998766d7f3dd0302758b8790819b50a51b98e57 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 7 Apr 2026 14:57:01 -0700 Subject: [PATCH 25/27] Fix bug in div, nondiv flow type definitions --- .../ocean/sphere_transport/resources/flow_types.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py index 4be555fcaa..ec224e97c1 100644 --- a/polaris/tasks/ocean/sphere_transport/resources/flow_types.py +++ b/polaris/tasks/ocean/sphere_transport/resources/flow_types.py @@ -43,10 +43,10 @@ def flow_nondivergent(t, lon, lat, u_0, tau): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (u_0 * earth_radius / tau) * ( - (sin(lon_p) ** 2) * sin(2 * lat) * cost + 2.0 * pi * coslat + u = (earth_radius / tau) * ( + u_0 * (sin(lon_p) ** 2) * sin(2 * -lat) * cost + 2.0 * pi * coslat ) - v = (u_0 * earth_radius / tau) * sin(2 * lon_p) * coslat * cost + v = (u_0 * earth_radius / tau) * sin(2 * lon_p + pi) * coslat * cost return u, v @@ -82,8 +82,8 @@ def flow_divergent(t, lon, lat, u_0, tau): lon_p = lon - 2.0 * pi * t / tau coslat = cos(lat) cost = cos(pi * t / tau) - u = (u_0 * earth_radius / tau) * ( - -(sin(lon_p / 2) ** 2) * sin(2 * lat) * (coslat**2) * cost + u = (earth_radius / tau) * ( + -u_0 * (sin(lon_p / 2) ** 2) * sin(2 * lat) * (coslat**2) * cost + 2.0 * pi * coslat ) v = (u_0 * earth_radius / (2 * tau)) * (sin(lon_p) * (coslat**3) * cost) From 0ee6fe94d10109d568495982d006ac7da8cfa019 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Apr 2026 10:03:50 -0700 Subject: [PATCH 26/27] Reduce the timestep for all sphere_transport tests --- polaris/tasks/ocean/sphere_transport/sphere_transport.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg b/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg index 4f0432ecf2..bdb4026989 100644 --- a/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg +++ b/polaris/tasks/ocean/sphere_transport/sphere_transport.cfg @@ -40,7 +40,7 @@ error_type = l2 time_integrator = RK4 # RK4 time step per resolution (s/km), since dt is proportional to resolution -rk4_dt_per_km = 8.0 +rk4_dt_per_km = 1.0 # Run duration in hours run_duration = ${sphere_transport:vel_pd} From f2b323ba4f5c4e582b52f06ee60df98fcaa38d06 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Apr 2026 10:04:50 -0700 Subject: [PATCH 27/27] Update E3SM submodule --- e3sm_submodules/Omega | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/e3sm_submodules/Omega b/e3sm_submodules/Omega index 589689d3e9..852572a331 160000 --- a/e3sm_submodules/Omega +++ b/e3sm_submodules/Omega @@ -1 +1 @@ -Subproject commit 589689d3e9dbff201e40fe7c11d7297aeafc63b1 +Subproject commit 852572a3314a1f36b0b727c7772f752e1af2d7e6