diff --git a/.github/workflows/omega-buildnml-workflow.yml b/.github/workflows/omega-buildnml-workflow.yml new file mode 100644 index 000000000000..9f8eccc478c5 --- /dev/null +++ b/.github/workflows/omega-buildnml-workflow.yml @@ -0,0 +1,111 @@ +name: omega-buildnml + +on: + push: + branches: [develop] + paths: + - 'components/omega/cime_config/**' + - '!components/omega/cime_config/*.xml' + - 'components/omega/configs/Default.yml' + - 'components/omega/dev-conda.txt' + - '.github/workflows/omega-buildnml-workflow.yml' + + pull_request: + branches: [develop] + paths: + - 'components/omega/cime_config/**' + - '!components/omega/cime_config/*.xml' + - 'components/omega/configs/Default.yml' + - 'components/omega/dev-conda.txt' + - '.github/workflows/omega-buildnml-workflow.yml' + + workflow_dispatch: + +concurrency: + # Cancel in progress runs testing the same git ref + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +env: + PYTHON_VERSION: "3.13" + +jobs: + validate-config: + if: ${{ github.repository == 'E3SM-Project/E3SM' || + github.repository == 'E3SM-Project/Omega' }} + name: validate configuration files + runs-on: ubuntu-latest + timeout-minutes: 20 + defaults: + run: + shell: bash -l {0} + steps: + - name: Checkout Code Repository + uses: actions/checkout@v4 + with: + persist-credentials: false + + - name: Set up Conda Environment + uses: conda-incubator/setup-miniconda@v3 + with: + activate-environment: "omega_ci" + miniforge-version: latest + channels: conda-forge + channel-priority: strict + auto-update-conda: true + python-version: ${{ env.PYTHON_VERSION }} + + - name: Install dependencies + run: | + conda create -n omega_dev --file components/omega/dev-conda.txt \ + python=${{ env.PYTHON_VERSION }} + conda activate omega_dev + conda list + + - name: Validate Omega configuration files + run: | + conda activate omega_dev + cd components/omega/cime_config + ./validate_config.py + + unit-tests: + if: ${{ github.repository == 'E3SM-Project/E3SM' || + github.repository == 'E3SM-Project/Omega' }} + name: unit tests (python ${{ matrix.python-version }}) + runs-on: ubuntu-latest + timeout-minutes: 20 + strategy: + fail-fast: false + matrix: + python-version: ["3.9", "3.10", "3.11", "3.12", "3.13", "3.14"] + defaults: + run: + shell: bash -l {0} + steps: + - name: Checkout Code Repository + uses: actions/checkout@v4 + with: + persist-credentials: false + + - name: Set up Conda Environment + uses: conda-incubator/setup-miniconda@v3 + with: + activate-environment: "omega_ci" + miniforge-version: latest + channels: conda-forge + channel-priority: strict + auto-update-conda: true + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + conda create -n omega_dev --file components/omega/dev-conda.txt \ + python=${{ matrix.python-version }} + conda activate omega_dev + conda list + + - name: Run omega_buildnml unit tests + run: | + conda activate omega_dev + cd components/omega/cime_config + python -m pytest omega_buildnml/tests -v diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b9a45151ca7f..2a9bf1c75216 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -43,7 +43,7 @@ repos: - id: mypy args: ["--config=components/omega/mypy.cfg", "--show-error-codes"] verbose: true - additional_dependencies: ['types-requests'] + additional_dependencies: ['types-requests', 'types-PyYAML'] exclude: ^components/omega/cime_config/build.* # Can run individually with `pre-commit run fortitude --all-files` diff --git a/cime_config/allactive/config_compsets.xml b/cime_config/allactive/config_compsets.xml index c88a89a782f4..b8cdf3115a0e 100755 --- a/cime_config/allactive/config_compsets.xml +++ b/cime_config/allactive/config_compsets.xml @@ -85,6 +85,11 @@ 1850_EAM%CMIP6_ELM%CNPRDCTCBCTOP_MPASSI_MPASO_MOSART_SGLC_SWAV + + WCYCL1850NS-OMEGA + 1850_EAM%CMIP6_ELM%CNPRDCTCBCTOP_MPASSI_OMEGA_MOSART_SGLC_SWAV + + WCYCL2010NS diff --git a/components/cmake/build_omega.cmake b/components/cmake/build_omega.cmake index eff6c4061e71..fc2a9b871f17 100644 --- a/components/cmake/build_omega.cmake +++ b/components/cmake/build_omega.cmake @@ -1,8 +1,12 @@ function(build_omega) - # Set CIME source path relative to components - set(CIMESRC_PATH "../cime/src") + if (COMP_NAMES MATCHES ".*omega.*") - add_subdirectory("omega") + # Set CIME source path relative to components + set(CIMESRC_PATH "../cime/src") + + add_subdirectory("omega") + + endif() endfunction(build_omega) diff --git a/components/eam/bld/config_files/definition.xml b/components/eam/bld/config_files/definition.xml index 023903d11e57..39926b69a45c 100644 --- a/components/eam/bld/config_files/definition.xml +++ b/components/eam/bld/config_files/definition.xml @@ -110,7 +110,7 @@ Use clm, stub land or no lnd model in cam build: clm, slnd, none Use rtm, stub runoff or no runoff model in cam build: rtm, srof, none - + Use data ocean model (docn or dom), stub ocean (socn), or aqua planet ocean (aquaplanet) in cam build. When built from the CESM scripts the value of ocn may be set to pop, but this doesn't impact how CAM is built. diff --git a/components/eam/bld/configure b/components/eam/bld/configure index 304d595522c3..b01336a96388 100755 --- a/components/eam/bld/configure +++ b/components/eam/bld/configure @@ -212,7 +212,7 @@ OPTIONS Options for surface components used in standalone EAM mode: -ice Build EAM with sea ice model [cice | sice | none ]. Default: cice. - -ocn Build EAM with ocean model [docn | dom | socn | aquaplanet | pop | mpaso]. Default: docn. + -ocn Build EAM with ocean model [docn | dom | socn | aquaplanet | pop | mpaso | omega]. Default: docn. -lnd Build EAM with land model [clm | slnd | none]. Default: clm. -rof Build EAM with runoff model [rtm | srof | none]. Default: rtm. diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index b42f03b8564f..ded42bf7f582 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -180,6 +180,7 @@ 0 0 3 +3 diff --git a/components/omega/cime_config/buildnml b/components/omega/cime_config/buildnml index 9f86189bdb26..40ad253f93dc 100755 --- a/components/omega/cime_config/buildnml +++ b/components/omega/cime_config/buildnml @@ -7,29 +7,20 @@ Namelist creator for E3SM's omega component import sys from pathlib import Path -import yaml from CIME.buildnml import parse_input from CIME.case import Case from CIME.utils import SharedArea, expect, safe_copy - -nan_time_instant_str = "9999-12-31_00:00:00" - -omega_supported_grids = ["oQU240", "EC30to60E2r2"] - - -# YAML representer to preserve order of dicts when dumping to yaml -class OrderedSafeDumper(yaml.SafeDumper): - pass - - -def _dict_representer(dumper, data): - return dumper.represent_mapping( - yaml.resolver.BaseResolver.DEFAULT_MAPPING_TAG, data.items() - ) - - -OrderedSafeDumper.add_representer(dict, _dict_representer) -# end of YAML representer related code +from omega_buildnml import ( + build_omega_config, + build_runtime_overrides, + read_config_overrides, + read_default_config, + read_input_files_config, + read_user_overrides, + resolve_streams_files, + write_input_data_list, + write_yaml_mapping, +) def buildnml(case, caseroot, compname): @@ -52,20 +43,53 @@ def buildnml(case, caseroot, compname): run_dir = Path(case.get_value("RUNDIR")) case_root = Path(case.get_value("CASEROOT")) + case_build = Path(case.get_value("CASEBUILD")) omega_root = Path(case.get_value("SRCROOT")) / "components/omega" + din_loc_root = Path(case.get_value("DIN_LOC_ROOT")) + + case_name = case.get_value("CASE") + calendar = case.get_value("CALENDAR") + mesh_name = case.get_value("OCN_GRID") + continue_run = case.get_value("CONTINUE_RUN") + + # read input_files.yaml and find input files needed for this mesh + input_files = read_input_files_config(mesh_name=mesh_name) + streams_files = resolve_streams_files( + input_files=input_files, mesh_name=mesh_name, din_loc_root=din_loc_root + ) + + # convert CIME case configuration options to Omega config snippets + runtime_overrides = build_runtime_overrides( + calendar=calendar, + continue_run=continue_run, + case_name=case_name, + streams_files=streams_files, + ) - path_to_default_config = omega_root / "configs/Default.yml" + # read Omega's Defaults.yml, which serve as base configuration + defaults = read_default_config(omega_root / "configs/Default.yml") - config = _build_omega_config(path_to_default_config, case) + # read coupled simulation and mesh-specific Omega config snippets + config_overrides = read_config_overrides(mesh_name=mesh_name) - # write yaml file and copy to CaseDocs within a SharedArea context so that - # correct file permissions are set + # read user_nl_omega for any user-specified overrides + user_overrides = read_user_overrides(case_root / "user_nl_omega") + + # Using the defaults and overrides, build the final Omega config dictionary + config = build_omega_config( + defaults=defaults, + coupled_overrides=config_overrides["coupled"], + mesh_overrides=config_overrides.get("meshes", {}).get(mesh_name, {}), + runtime_overrides=runtime_overrides, + user_overrides=user_overrides, + ) + + # write files within SharedArea context so correct file permissions are set with SharedArea(): output_file = run_dir / "omega.yml" # write the updated yaml file to the run directory - with open(output_file, "w") as file: - yaml.dump(config, file, Dumper=OrderedSafeDumper) + write_yaml_mapping(config, output_file) # archive yaml file in CaseDocs case_docs = case_root / "CaseDocs" @@ -73,217 +97,7 @@ def buildnml(case, caseroot, compname): safe_copy(output_file, case_docs / "omega.yml") # needs to be called within buildnml - build_input_data_list(case, compname) - - -def build_input_data_list(case, compname): - """ - Add files to omega.input_data_list so they can be automatically downloaded - """ - expect(compname == "omega", compname) - - casebuild = case.get_value("CASEBUILD") - input_data_list = [] - - # Add the initial state file to the input data list - initial_state_file = _get_file_name(case, "InitialState") - input_data_list.append(f"initial_state = {initial_state_file}") - - # TODO: For now all information is contained in a single input file - # # Add the horizontal mesh file to the input data list - # horz_mesh_file = _get_file_name(case, "HorzMeshIn") - # input_data_list.append(horz_mesh_file) - - # # Add the initial vertical coordinate file to the input data list - # initial_vert_coord_file = _get_file_name(case, "InitialVertCoord") - # input_data_list.append(initial_vert_coord_file) - - with open(Path(casebuild) / "omega.input_data_list", "w") as f: - for input_file in input_data_list: - f.write(f"{input_file}\n") - - -def _build_omega_config(path_to_default, case): - """ - Build the omega config file based on case configuration - - Parameters - ---------- - path_to_default : str - Path to the default omega config file. - case : Case - The case object containing the configuration for the current run. - - Returns - ------- - dict - The updated omega config dictionary. - """ - with open(path_to_default, "r") as file: - default_config = yaml.load(file, Loader=yaml.SafeLoader) - # unpack the Omega section of the default config - default_config = default_config["Omega"] - - # start with a copy of the default config, unpack the Omega section - config = default_config.copy() - - # update the TimeIntegration section of the config - config["TimeIntegration"] = _build_time_integration_config( - default_config["TimeIntegration"], case - ) - - # update the IOStreams section of the config - config["IOStreams"] = _build_io_streams_config( - default_config["IOStreams"], case - ) - - # repack everything withn a top-level Omega section - return {"Omega": config} - - -# Helper function to build the various sections of the omega streams -def _get_file_name(case, stream): - """ """ - # TODO: this will be conditional set based on grid and compset. It - # also needs to check RUN_TYPE: for "hybrid"/"branch" runs this must - # point at the RUN_REFCASE restart file (see mpas-ocean's - # cime_config/buildnml, which swaps in run_refcase's restart file - # when run_type in ("hybrid", "branch")), not the default IC file. - - din_loc = Path(case.get_value("DIN_LOC_ROOT")) - ocn_grid = case.get_value("OCN_GRID") - - if ocn_grid not in omega_supported_grids: - msg = f"Unsupported OCN_GRID: {ocn_grid}" - raise ValueError(msg) - - # TODO: this will return a different value based on the stream name - # For now mesh file contains all information needed in single file - omega_dir = din_loc / "ocn/omega" - - if ocn_grid == "oQU240": - mesh_dir = omega_dir / "oQU240" - return str(mesh_dir / "ocean.QU.240km.151209.teos10.260807.nc") - elif ocn_grid == "EC30to60E2r2": - mesh_dir = omega_dir / "EC30to60E2r2" - return str(mesh_dir / "ocean.EC30to60E2r2.200908.teos10.260807.nc") - else: - msg = f"Unsupported OCN_GRID: {ocn_grid}" - raise ValueError(msg) - - -def _update_read_io_streams(config, case, stream): - """Update the read streams - - RUN_TYPE stays fixed for the life of a case, so it cannot tell us - whether *this* job submission should read the initial state or a - restart file. CONTINUE_RUN is FALSE only for the very first job of - a run (including hybrid/branch runs) and TRUE for every subsequent - continuation, so it is the correct flag to gate on here. - """ - continue_run = case.get_value("CONTINUE_RUN") - - # RestartRead has inverse logic to InitialState so negate value - if stream == "RestartRead": - continue_run = not continue_run - - config[stream]["Filename"] = _get_file_name(case, stream) - - # Leave `OnstartUp` default for HorzMeshIn and InitialVertCoord streams - if stream not in ["HorzMeshIn", "InitialVertCoord"]: - config[stream]["FreqUnits"] = "Never" if continue_run else "OnStartUp" - - # Make sure pointer, with the correct filename, is used by restart reads - if stream == "RestartRead": - config[stream]["UsePointerFile"] = True - config[stream]["PointerFilename"] = "rpointer.ocn" - - # Default.yml sets UseStartEnd/StartTime/EndTime for RestartRead as a - # workaround so the standalone driver (whose FreqUnits is always - # "OnStartup") can still skip reading a restart on the very first - # startup job. The coupled driver already gates this explicitly via - # CONTINUE_RUN above, so the UseStartEnd gate must be disabled here: - # otherwise it depends on the model clock's StartAlarm/EndAlarm - # already reflecting the coupler-provided current time, which is not - # guaranteed to be true this early in ocnInit1. - config[stream]["UseStartEnd"] = False - - return config - - -def _build_restart_write_io_stream(case): - """Coupler dictates when to write restart files, so we just need to set up - the stream here. - """ - case_name = case.get_value("CASE") - - return { - "UsePointerFile": True, - "PointerFilename": "rpointer.ocn", - "Filename": f"{case_name}.omega.r.$Y-$M-$D_$h.$m.$s", - "Mode": "write", - "IfExists": "replace", - "Precision": "double", - "Freq": 1, - "FreqUnits": "OnDemand", - "UseStartEnd": False, - "Contents": ["Restart"], - } - - -def _build_io_streams_config(default_io_streams, case): - """Build the io_streams section of the omega config file""" - io_streams = default_io_streams.copy() - - # Update InitialState, HorzMeshIn, InitialVertCoord, and RestartRead - io_streams = _update_read_io_streams(io_streams, case, "InitialState") - io_streams = _update_read_io_streams(io_streams, case, "HorzMeshIn") - io_streams = _update_read_io_streams(io_streams, case, "InitialVertCoord") - io_streams = _update_read_io_streams(io_streams, case, "RestartRead") - - # Add restart write stream - io_streams["RestartWrite"] = _build_restart_write_io_stream(case) - - # Never read forcing stream in coupled runs - io_streams["Forcing"]["FreqUnits"] = "Never" - - return io_streams - - -def _build_time_integration_config(default_time_integration, case): - """Start and stop time are set by the coupler. - The time step is mesh dependent. - """ - ocn_grid = case.get_value("OCN_GRID") - calendar = case.get_value("CALENDAR") - - time_integration = default_time_integration.copy() - - # convert from CIME calendar to Omega calendar names - if calendar == "NO_LEAP": - time_integration["CalendarType"] = "No Leap" - elif calendar == "GREGORIAN": - time_integration["CalendarType"] = "Gregorian" - else: - msg = f"Unsupported calendar type: {calendar}" - raise ValueError(msg) - - # Set start and stop time to invalid values in order to make clear to users - # these values come from the coupler at runtime - time_integration["StartTime"] = nan_time_instant_str - time_integration["StopTime"] = nan_time_instant_str - # RunDuration is not used in coupled runs - time_integration["RunDuration"] = "none" - - # All meshes use RK4 timestepping (for now) - time_integration["TimeStepper"] = "RungeKutta4" - - if ocn_grid == "oQU240": - time_integration["TimeStep"] = "0000_00:05:00" - elif ocn_grid == "EC30to60E2r2": - time_integration["TimeStep"] = "0000_00:01:00" - - return time_integration + write_input_data_list(streams_files, case_build) def _main_func(): diff --git a/components/omega/cime_config/config_compsets.xml b/components/omega/cime_config/config_compsets.xml index 303a3a59e2d2..57ec451d46da 100644 --- a/components/omega/cime_config/config_compsets.xml +++ b/components/omega/cime_config/config_compsets.xml @@ -16,9 +16,19 @@ - COMEGA-JRA1p4 - 2000_DATM%JRA-1p4-2018_SLND_DICE%SSMI_OMEGA%DATMFORCED_DROF%JRA-1p4-2018_SGLC_SWAV + COMEGA-JRA1p5 + 2000_DATM%JRA-1p5_SLND_DICE%SSMI_OMEGA%DATMFORCED_DROF%JRA-1p5_SGLC_SWAV Experimental, under development + + GOMEGA-IAF + 2000_DATM%IAF_SLND_MPASSI_OMEGA%DATMFORCED_DROF%IAF_SGLC_SWAV + + + + GOMEGA-JRA1p5 + 2000_DATM%JRA-1p5_SLND_MPASSI_OMEGA%DATMFORCED_DROF%JRA-1p5_SGLC_SWAV + + diff --git a/components/omega/cime_config/omega_buildnml/README.md b/components/omega/cime_config/omega_buildnml/README.md new file mode 100644 index 000000000000..3990c294df6e --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/README.md @@ -0,0 +1,66 @@ +# `omega_buildnml` + +Python package that builds Omega's runtime configuration (`omega.yml`) during +`case.setup`/`case.submit`. + +Special care is taken to validate the input configuration files against +`configs/Default.yml`, so mistakes (typos, unknown options, options CIME +controls exclusively) are caught early rather than surfacing as a confusing +runtime error. `buildnml` only validates the mesh requested by the case +being built; CI (`validate_config.py`) validates every mesh and override +entry in the packaged configuration files. See render documentation for +more details. + +> [!CAUTION] +> **No external dependencies are allowed in this package** +> +> Dependencies **must** be limited to the python (3.9+) standard library and `PyYAML`. +> CIME does not provide its own environment, it just assumes the case's python +> environment it is being called from already has PyYAML available. In practice +> this is often handled by loading a python module on an HPC system that +> includes `PyYAML`. + +## Modules + +``` +omega_buildnml +├── __init__.py # public API is defined by contents of __all__ +├── _types.py # shared type aliases +├── config.py # merges configuration layers and resolves mesh input files +├── read_write.py # reads/writes YAML config files, including packaged data +├── validate.py # validates configuration layers +├── data/ # packaged YAML data (see data/README.md) +└── tests/ # unit tests +``` + +## Development environment + +For now, we use the same `conda` [environment](../../dev-conda.txt) used for +linting Omega source code. The development environment primarily provides +`pytest` needed for running the unit tests. The actual code defined in this +package is limited to the python standard library and PyYAML. + +```bash +conda create -n omega_dev --file ../../dev-conda.txt +conda activate omega_dev +pre-commit install +``` + +## Testing + +To validate the configuration files and run unit tests, do: +```bash +cd components/omega/cime_config +./validate_config.py +pytest omega_buildnml/tests +``` + +## TODO + +- [ ] Validate `IOStreams` conditionally on what's set within a stream + - For example: if `UseStartEnd` is true, both `StartTime` and `EndTime` + must be provided. +- [ ] Use `difflib` to suggest close matches for unknown/misspelled options. +- [ ] Check that mesh input files exist in the remote input data database + - Gate failures on whether database is reachable at all. + (Prevents false positives when LCRC is down for monthly maintenance) diff --git a/components/omega/cime_config/omega_buildnml/__init__.py b/components/omega/cime_config/omega_buildnml/__init__.py new file mode 100644 index 000000000000..394c187dc5e8 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/__init__.py @@ -0,0 +1,25 @@ +from .config import ( + build_omega_config, + build_runtime_overrides, + resolve_streams_files, +) +from .read_write import ( + read_config_overrides, + read_default_config, + read_input_files_config, + read_user_overrides, + write_input_data_list, + write_yaml_mapping, +) + +__all__ = [ + "build_omega_config", + "build_runtime_overrides", + "read_config_overrides", + "read_default_config", + "read_input_files_config", + "read_user_overrides", + "resolve_streams_files", + "write_input_data_list", + "write_yaml_mapping", +] diff --git a/components/omega/cime_config/omega_buildnml/_types.py b/components/omega/cime_config/omega_buildnml/_types.py new file mode 100644 index 000000000000..4491a1b2713e --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/_types.py @@ -0,0 +1,5 @@ +from os import PathLike as OsPathLike +from typing import Any, Union + +YamlMapping = dict[str, Any] +PathLike = Union[str, OsPathLike[str]] diff --git a/components/omega/cime_config/omega_buildnml/config.py b/components/omega/cime_config/omega_buildnml/config.py new file mode 100644 index 000000000000..be92808dd8eb --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/config.py @@ -0,0 +1,184 @@ +from collections.abc import Mapping +from copy import deepcopy +from pathlib import Path +from typing import Any + +from ._types import PathLike, YamlMapping + + +def build_omega_config( + defaults: YamlMapping, + coupled_overrides: YamlMapping, + mesh_overrides: YamlMapping, + runtime_overrides: YamlMapping, + user_overrides: YamlMapping, +) -> YamlMapping: + """ + Build the Omega configuration dictionary. + + Parameters: + ----------- + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + coupled_overrides : dict[str, Any] + Coupled model overrides, loaded from cime_config/config_overrides.yaml + mesh_overrides: dict[str, Any] + Mesh specific overrides, loaded from cime_config/mesh_overrides.yaml + runtime_overrides : dict[str, Any] + Runtime specific overrides, based on CIME case configuration. + user_overrides : dict[str, Any] + User specified overrides, loaded from user_nl_omega. + + Returns: + -------- + dict[str, Any]: + Final Omega configuration dictionary. + """ + if "Omega" in defaults: + defaults = defaults["Omega"] + + config = deepcopy(defaults) + config = _deep_merge(config, coupled_overrides) + config = _deep_merge(config, mesh_overrides) + config = _deep_merge(config, runtime_overrides) + config = _deep_merge(config, user_overrides) + + return {"Omega": config} + + +def resolve_streams_files( + input_files: YamlMapping, mesh_name: str, din_loc_root: PathLike +) -> dict[str, str]: + """ + Resolve input filenames for all Omega IOStreams + + Uses mesh name and input_files.yaml to determine the correct input files + for each IOStream. + + Parameters: + ----------- + input_files : dict[str, str] + Parsed content of cime_config/input_files.yaml. + mesh_name : str + CIME ocean grid name + din_loc_root : Path + Path to root of E3SM input data directory. + + Returns: + -------- + dict[str, str] + Mapping of Omega IOStream names to resolved input filenames. + """ + meshes: YamlMapping = input_files["meshes"] + + mesh_dir = Path(din_loc_root) / "ocn" / "omega" / mesh_name + streams_files = {} + + for input_group in meshes[mesh_name]["inputs"]: + resolved_file_path = mesh_dir / input_group["file"] + + for stream in input_group["streams"]: + streams_files[stream] = str(resolved_file_path) + + return streams_files + + +def build_runtime_overrides( + calendar: str, + continue_run: bool, + case_name: str, + streams_files: dict[str, str], +) -> YamlMapping: + """ + Build the runtime override dictionary from the CIME case configuration. + + Parameters: + ----------- + calendar : str + CIME calendar type (i.e., "NO_LEAP", "GREGORIAN"). + continue_run : bool + Whether to continue a previous run. + case_name : str + CIME case name. + streams_files : dict[str, str] + Resolved input filenames keyed by Omega IOStream name. + + Returns: + -------- + dict[str, Any] + Runtime overrides dictionary. + """ + io_overrides = { + stream_name: {"Filename": filename} + for stream_name, filename in streams_files.items() + } + + if continue_run: + io_overrides["InitialState"]["FreqUnits"] = "Never" + io_overrides["RestartRead"] = {"FreqUnits": "OnStartup"} + else: + io_overrides["InitialState"]["FreqUnits"] = "OnStartup" + io_overrides["RestartRead"] = {"FreqUnits": "Never"} + + io_overrides["RestartWrite"] = { + "Filename": f"{case_name}.omega.r.$Y-$M-$D_$h.$m.$s" + } + + return { + "TimeIntegration": {"CalendarType": _to_omega_calendar(calendar)}, + "IOStreams": io_overrides, + } + + +def _deep_merge( + base: Mapping[str, Any], + override: Mapping[str, Any], +) -> YamlMapping: + """ + Recursively merge two mappings without modifying either input. + + Nested mappings are merged. All other values, including lists, are replaced + by the override value. + + Parameters: + ----------- + base (dict[str, Any]) + The base mapping + override (dict[str, Any]): + The mapping with the override values + + Returns: + dict[str, Any]: Merged mapping + """ + merged: YamlMapping = deepcopy(dict(base)) + + for key, override_value in override.items(): + base_value = merged.get(key) + + if isinstance(base_value, Mapping) and isinstance( + override_value, Mapping + ): + merged[key] = _deep_merge(base_value, override_value) + else: + merged[key] = deepcopy(override_value) + + return merged + + +def _to_omega_calendar(calendar: str) -> str: + """ + Convert a CIME calendar string to corresponding Omega calendar string + + Args: + calendar (str): CIME calendar string. + + Returns: + str: Omega calendar string. + """ + if calendar == "NO_LEAP": + return "No Leap" + elif calendar == "GREGORIAN": + return "Gregorian" + else: + msg = f"Unsupported calendar type: {calendar}" + raise ValueError(msg) diff --git a/components/omega/cime_config/omega_buildnml/data/README.md b/components/omega/cime_config/omega_buildnml/data/README.md new file mode 100644 index 000000000000..ea37f688d8ec --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/data/README.md @@ -0,0 +1,14 @@ +# `data` + +Packaged YAML data read via `importlib.resources` in `read_write.py`. See +render documentation for how to add a new mesh. + +## Files + +``` +data +├── input_files.yaml # maps each mesh to its input files +└── config_overrides.yaml # automatically applied configuration overrides +``` + +Both files are validated against `configs/Default.yml` when read. diff --git a/components/omega/cime_config/omega_buildnml/data/config_overrides.yaml b/components/omega/cime_config/omega_buildnml/data/config_overrides.yaml new file mode 100644 index 000000000000..2c32a1faa4d9 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/data/config_overrides.yaml @@ -0,0 +1,41 @@ +coupled: + TimeIntegration: + TimeStepper: RungeKutta4 + StartTime: 9999-12-31_00:00:00 + StopTime: 9999-12-31_00:00:00 + RunDuration: none + + Tendencies: + SfcStressForcingTendencyEnable: true + SfcThicknessForcingTendencyEnable: true + SfcTracerForcingTendencyEnable: true + + IOStreams: + Forcing: + FreqUnits: Never + + RestartRead: + UsePointerFile: true + PointerFilename: rpointer.ocn + UseStartEnd: false + + RestartWrite: + UsePointerFile: true + PointerFilename: rpointer.ocn + Mode: write + IfExists: replace + Precision: double + Freq: 1 + FreqUnits: OnDemand + UseStartEnd: false + Contents: [Restart] + +meshes: + # Per-mesh overrides may set any option except `IOStreams`. + oQU240: + TimeIntegration: + TimeStep: "0000_00:05:00" + + EC30to60E2r2: + TimeIntegration: + TimeStep: "0000_00:01:00" diff --git a/components/omega/cime_config/omega_buildnml/data/input_files.yaml b/components/omega/cime_config/omega_buildnml/data/input_files.yaml new file mode 100644 index 000000000000..0764c6cebcb6 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/data/input_files.yaml @@ -0,0 +1,41 @@ +# ``meshes``, keyed by mesh name, is currently the only key read from this file. +# What feels like an unnecessary ``meshes``` key was intentional included to +# leave room for the additional keys (e.g. ``experiments``) in the future. +# For additional top-level keys to be supported, additional validation logic +# will need to be added to `omega_buildnml/validate.py` + +meshes: + oQU240: + inputs: + - file: ocean.QU.240km.151209.teos10.260807.nc + streams: [HorzMeshIn, InitialVertCoord, InitialState] + + EC30to60E2r2: + inputs: + - file: ocean.EC30to60E2r2.200908.teos10.260807.nc + streams: [HorzMeshIn, InitialVertCoord, InitialState] + +# A mesh may list multiple ``inputs`` entries when its streams are split +# across separate files: +# +# meshes: +# example_mesh: +# inputs: +# - file: example_mesh.horiz_mesh_and_vert_coord.nc +# streams: [HorzMeshIn, InitialVertCoord] +# - file: example_mesh.initial_condition.nc +# streams: [InitialState] + +# Sketch of what an experiments top-level key might look like. +# +# meshes: +# ... +# experiments: +# example_spinup: +# when: +# mesh: example_mesh +# compset: EXAMPLE_COMPSET +# +# inputs: +# - file: example_mesh.spinup_initial_condition.nc +# streams: [InitialState] diff --git a/components/omega/cime_config/omega_buildnml/read_write.py b/components/omega/cime_config/omega_buildnml/read_write.py new file mode 100644 index 000000000000..48a2fee0b2a2 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/read_write.py @@ -0,0 +1,295 @@ +from importlib import resources +from pathlib import Path +from typing import IO, Any, Hashable, Optional + +import yaml + +from ._types import PathLike, YamlMapping +from .validate import ( + validate_config_overrides, + validate_input_files_config, + validate_user_overrides, +) + +#: Path to Omega's default configuration file, relative to this package +DEFAULT_CONFIG_PATH = Path(__file__).parents[2] / "configs" / "Default.yml" + + +def read_default_config(path: PathLike) -> YamlMapping: + """ + Read Omega's default configuration from a YAML file on disk. + + Parameters: + ----------- + path : PathLike + Path to default config file (i.e. components/omega/config/Defaults.yml) + + Returns: + -------- + dict[str, Any] + The read configuration as a dictionary. + """ + path = Path(path) + + if not path.is_file(): + err_msg = f"{path} does not exist or is not a file" + raise FileNotFoundError(err_msg) + + with path.open("r", encoding="utf-8") as f: + config = _read_yaml_file(f) + + if "Omega" not in config: + err_msg = f"{path} does not contain a top-level 'Omega' section" + raise ValueError(err_msg) + + # TODO: check that the config is valid (e.g., required keys are present) + return config["Omega"] + + +def read_input_files_config(mesh_name: Optional[str] = None) -> YamlMapping: + """ + Read the input_files.yaml configuration file from the package resources. + + Parameters: + ----------- + mesh_name : str, optional + The name of the mesh to validate. If not provided, all mesh entries + will be validated. + + Returns: + -------- + dict[str, Any] + The read configuration as a dictionary. + """ + input_files = _read_packaged_yaml_file("input_files.yaml") + defaults = read_default_config(DEFAULT_CONFIG_PATH) + + return validate_input_files_config( + input_files, defaults, mesh_name=mesh_name + ) + + +def read_config_overrides(mesh_name: Optional[str] = None) -> YamlMapping: + """ + Read config_overrides.yaml configuration file from the package resources. + + Parameters: + ----------- + mesh_name : str, optional + The name of the mesh to validate. If not provided, all mesh entries + will be validated. + + Returns: + -------- + dict[str, Any] + The read configuration as a dictionary. + """ + config_overrides = _read_packaged_yaml_file("config_overrides.yaml") + input_files = read_input_files_config(mesh_name=mesh_name) + defaults = read_default_config(DEFAULT_CONFIG_PATH) + + return validate_config_overrides( + config_overrides, input_files, defaults, mesh_name=mesh_name + ) + + +def read_user_overrides(path: PathLike) -> YamlMapping: + """ + Read a case's user_nl_omega file. + + The overrides are expected to be a YAML snippet mirroring the structure + of Omega's default configuration. An empty file is not an error, it just + means the user has not overridden anything. + + Parameters: + ----------- + path : PathLike + Path to the case's ``user_nl_omega`` file. + + Returns: + -------- + dict[str, Any] + The read overrides as a dictionary. + """ + path = Path(path) + + if not path.is_file(): + err_msg = f"{path} does not exist or is not a file" + raise FileNotFoundError(err_msg) + + with path.open("r", encoding="utf-8") as f: + user_overrides = _read_yaml_file(f) + + # a user_nl_omega with no overrides in it parses as ``None`` + if user_overrides is None: + return {} + + if not isinstance(user_overrides, dict): + err_msg = f"{path} is not a mapping." + raise ValueError(err_msg) + + user_overrides = _unwrap_omega_section(user_overrides) + defaults = read_default_config(DEFAULT_CONFIG_PATH) + + return validate_user_overrides(user_overrides, defaults) + + +def write_yaml_mapping(mapping: YamlMapping, file_path: PathLike) -> None: + """ + Write a mapping to a YAML file. + + Parameters: + ----------- + mapping : dict[str, Any] + The mapping to write to the YAML file. + file_path : PathLike + The path to the output YAML file. + + Returns: + -------- + None + """ + with Path(file_path).open("w", encoding="utf-8") as f: + yaml.safe_dump(mapping, f, sort_keys=False, default_flow_style=False) + + +def write_input_data_list( + streams_files: dict[str, str], casebuild: PathLike +) -> None: + """ + Build omega.input_data_list + + Enables automatic retrieval of missing input files + + Parameters + ---------- + streams_files : dict[str, str] + Dict mapping stream names to their corresponding input file paths. + casebuild : Path + Path to the case build directory (i.e. CASEBUILD) + + Returns + ------- + None + """ + + unique_files = list(dict.fromkeys(streams_files.values())) + + input_data_list = [ + f"omega_input_{index} = {filename}" + for index, filename in enumerate(unique_files, start=1) + ] + + path = Path(casebuild) / "omega.input_data_list" + with path.open("w", encoding="utf-8") as f: + for input_file in input_data_list: + f.write(f"{input_file}\n") + + +class _UniqueKeyLoader(yaml.SafeLoader): + """ + YAML loader that rejects mappings containing duplicate keys. + + ``yaml.SafeLoader`` silently keeps the last value when a key is repeated, + so a copy/pasted or misspelled entry would quietly override an earlier + one instead of being reported. + """ + + def construct_mapping( + self, node: yaml.MappingNode, deep: bool = False + ) -> dict[Hashable, Any]: + keys = set() + + for key_node, _ in node.value: + key = self.construct_object(key_node, deep=deep) + + if key in keys: + err_msg = ( + f"Duplicate key {key!r} found:\n{key_node.start_mark}" + ) + raise ValueError(err_msg) + + keys.add(key) + + return super().construct_mapping(node, deep=deep) + + +def _unwrap_omega_section(config: YamlMapping) -> YamlMapping: + """ + Return the ``Omega`` section of a config, if one is present. + + Users are expected to wrap their overrides in a top-level ``Omega`` key, + matching Omega's default configuration, but the key is optional. When it + is present it must be the only top-level key, otherwise under indented + sections would be silently dropped. + + Parameters + ---------- + config : dict[str, Any] + Parsed content of a YAML file. + + Returns: + -------- + dict[str, Any] + The ``Omega`` section, or the config itself when there isn't one. + """ + if "Omega" not in config: + return config + + siblings = sorted(set(config) - {"Omega"}) + if siblings: + err_msg = ( + f"`Omega` must be the only top-level key. Found " + f"{', '.join(siblings)} alongside it, please check the " + f"indentation of your overrides." + ) + raise ValueError(err_msg) + + omega_section: YamlMapping = config["Omega"] + + # an ``Omega`` key with nothing under it parses as ``None`` + if omega_section is None: + return {} + + return omega_section + + +def _read_yaml_file(f: IO[str]) -> YamlMapping: + """ + Read a YAML mapping from a file. + + Parameters + ---------- + f : TextIO + File-like object to read the YAML mapping from. + Returns: + ------- + dict[str, Any] + Read YAML mapping. + """ + + return yaml.load(f, Loader=_UniqueKeyLoader) + + +def _read_packaged_yaml_file(file_name: str) -> YamlMapping: + """ + Read a YAML mapping packaged within omega_buildnml/data + + Parameters + ---------- + file_path : str + Name of the packaged YAML file to read. + + Returns: + -------- + dict[str, Any] + Read YAML mapping. + """ + resource = resources.files("omega_buildnml.data").joinpath(file_name) + + if not resource.is_file(): + err_msg = f"Packaged configuration file '{file_name}' not found." + raise FileNotFoundError(err_msg) + + with resource.open("r", encoding="utf-8") as f: + return _read_yaml_file(f) diff --git a/components/omega/cime_config/omega_buildnml/tests/__init__.py b/components/omega/cime_config/omega_buildnml/tests/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/components/omega/cime_config/omega_buildnml/tests/test_config.py b/components/omega/cime_config/omega_buildnml/tests/test_config.py new file mode 100644 index 000000000000..56ac06e8095a --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/tests/test_config.py @@ -0,0 +1,158 @@ +import pytest + +from omega_buildnml.config import build_omega_config +from omega_buildnml.read_write import DEFAULT_CONFIG_PATH, read_default_config +from omega_buildnml.validate import BLOCKED_OPTIONS, validate_user_overrides + + +@pytest.fixture +def defaults(): + """A minimal stand in for Omega's default configuration.""" + return { + "Omega": { + "TimeIntegration": { + "TimeStepper": "default", + "TimeStep": "0000_00:30:00", + }, + "IOStreams": { + "History": {"Freq": 1, "FreqUnits": "months"}, + }, + } + } + + +@pytest.fixture +def omega_defaults(): + """Omega's actual default configuration.""" + return read_default_config(DEFAULT_CONFIG_PATH) + + +def _time_stepper(source): + """Build an override setting ``TimeStepper`` to ``source``.""" + return {"TimeIntegration": {"TimeStepper": source}} + + +def _nested_option(option, value): + """Build a nested mapping from a dotted config path.""" + override = value + + for key in reversed(option.split(".")): + override = {key: override} + + return override + + +def _build(defaults, **layers): + """Build a config, defaulting any unspecified override layer to empty.""" + return build_omega_config( + defaults=defaults, + coupled_overrides=layers.get("coupled_overrides", {}), + mesh_overrides=layers.get("mesh_overrides", {}), + runtime_overrides=layers.get("runtime_overrides", {}), + user_overrides=layers.get("user_overrides", {}), + ) + + +def test_defaults_are_used_when_nothing_is_overridden(defaults): + config = _build(defaults) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "default" + + +def test_coupled_overrides_the_defaults(defaults): + config = _build(defaults, coupled_overrides=_time_stepper("coupled")) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "coupled" + + +def test_mesh_overrides_the_coupled_overrides(defaults): + config = _build( + defaults, + coupled_overrides=_time_stepper("coupled"), + mesh_overrides=_time_stepper("mesh"), + ) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "mesh" + + +def test_runtime_overrides_the_mesh_overrides(defaults): + config = _build( + defaults, + mesh_overrides=_time_stepper("mesh"), + runtime_overrides=_time_stepper("runtime"), + ) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "runtime" + + +def test_user_overrides_the_runtime_overrides(defaults): + config = _build( + defaults, + runtime_overrides=_time_stepper("runtime"), + user_overrides=_time_stepper("user"), + ) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "user" + + +def test_user_overrides_are_applied_last(defaults): + config = _build( + defaults, + coupled_overrides=_time_stepper("coupled"), + mesh_overrides=_time_stepper("mesh"), + runtime_overrides=_time_stepper("runtime"), + user_overrides=_time_stepper("user"), + ) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "user" + + +def test_options_that_are_not_overridden_are_preserved(defaults): + config = _build(defaults, user_overrides=_time_stepper("user")) + + time_integration = config["Omega"]["TimeIntegration"] + + assert time_integration["TimeStep"] == "0000_00:30:00" + + +def test_overrides_do_not_modify_their_inputs(defaults): + user_overrides = _time_stepper("user") + + _build(defaults, user_overrides=user_overrides) + + assert defaults["Omega"]["TimeIntegration"]["TimeStepper"] == "default" + assert user_overrides == _time_stepper("user") + + +def test_defaults_are_accepted_without_an_omega_section(defaults): + config = _build(defaults["Omega"], user_overrides=_time_stepper("user")) + + assert config["Omega"]["TimeIntegration"]["TimeStepper"] == "user" + + +@pytest.mark.parametrize("blocked_option", sorted(BLOCKED_OPTIONS)) +def test_blocked_options_are_rejected_before_they_are_merged( + omega_defaults, blocked_option +): + """ + Blocked options are enforced by validation, not by the merge order. + + User overrides are merged last, so a blocked option would win if it ever + reached ``build_omega_config``. + """ + user_overrides = _nested_option(blocked_option, "user") + + with pytest.raises(ValueError, match="cannot be overridden"): + validate_user_overrides(user_overrides, omega_defaults) + + +@pytest.mark.parametrize("blocked_option", sorted(BLOCKED_OPTIONS)) +def test_blocked_options_are_defined_in_the_defaults( + omega_defaults, blocked_option +): + """``BLOCKED_OPTIONS`` can drift from the options Omega defines.""" + option = omega_defaults + + for key in blocked_option.split("."): + assert key in option, f"{blocked_option} is not in the defaults" + option = option[key] diff --git a/components/omega/cime_config/omega_buildnml/tests/test_read_write.py b/components/omega/cime_config/omega_buildnml/tests/test_read_write.py new file mode 100644 index 000000000000..5809af5020cd --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/tests/test_read_write.py @@ -0,0 +1,158 @@ +import pytest +import yaml + +from omega_buildnml.read_write import ( + _read_yaml_file, + _unwrap_omega_section, + read_user_overrides, +) + + +@pytest.fixture +def user_nl(tmp_path): + """Write a mapping to a user_nl_omega file, returning the path.""" + + def _write(overrides): + path = tmp_path / "user_nl_omega" + + with path.open("w", encoding="utf-8") as f: + yaml.safe_dump(overrides, f) + + return path + + return _write + + +@pytest.fixture +def malformed_yaml(tmp_path): + """ + Write a YAML file that a mapping cannot express. + + Duplicate keys, comments, and indentation mistakes are all lost when a + mapping is dumped, so they have to be written out as text. + """ + + def _write(contents): + path = tmp_path / "sample.yaml" + path.write_text(contents, encoding="utf-8") + + return path + + return _write + + +def _read_yaml(path): + """Read a YAML file the way the package reads its config files.""" + with path.open("r", encoding="utf-8") as f: + return _read_yaml_file(f) + + +def test_mappings_without_duplicates_are_read(user_nl): + overrides = { + "TimeIntegration": { + "TimeStep": "0000_00:30:00", + "TimeStepper": "Forward-Backward", + } + } + + assert _read_yaml(user_nl(overrides)) == overrides + + +def test_duplicate_top_level_keys_are_rejected(malformed_yaml): + path = malformed_yaml( + "TimeIntegration:\n TimeStep: a\nTimeIntegration:\n TimeStep: b\n" + ) + + with pytest.raises(ValueError, match="Duplicate key 'TimeIntegration'"): + _read_yaml(path) + + +def test_duplicate_nested_keys_are_rejected(malformed_yaml): + path = malformed_yaml("TimeIntegration:\n TimeStep: a\n TimeStep: b\n") + + with pytest.raises(ValueError, match="Duplicate key 'TimeStep'"): + _read_yaml(path) + + +def test_duplicate_keys_within_a_list_are_rejected(malformed_yaml): + path = malformed_yaml("inputs:\n - file: a.nc\n file: b.nc\n") + + with pytest.raises(ValueError, match="Duplicate key 'file'"): + _read_yaml(path) + + +def test_duplicate_key_errors_report_where_they_were_found(malformed_yaml): + path = malformed_yaml("Tracers:\n Base: a\n Base: b\n") + + with pytest.raises(ValueError, match="sample.yaml"): + _read_yaml(path) + + +def test_overrides_are_returned_when_there_is_no_omega_section(): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:30:00"}} + + assert _unwrap_omega_section(overrides) == overrides + + +def test_the_omega_section_is_unwrapped(): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:30:00"}} + + assert _unwrap_omega_section({"Omega": overrides}) == overrides + + +def test_an_empty_omega_section_is_not_an_error(): + assert _unwrap_omega_section({"Omega": None}) == {} + + +def test_keys_alongside_the_omega_section_are_rejected(): + overrides = { + "Omega": {"IOStreams": {}}, + "TimeIntegration": {"TimeStep": "0000_00:30:00"}, + } + + with pytest.raises(ValueError, match="only top-level key"): + _unwrap_omega_section(overrides) + + +def test_an_empty_user_nl_omega_is_not_an_error(malformed_yaml): + assert read_user_overrides(malformed_yaml("")) == {} + + +def test_a_user_nl_omega_with_only_comments_is_not_an_error(malformed_yaml): + assert read_user_overrides(malformed_yaml("# no overrides here\n")) == {} + + +def test_user_overrides_are_read_without_an_omega_section(user_nl): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:10:00"}} + + assert read_user_overrides(user_nl(overrides)) == overrides + + +def test_user_overrides_are_read_with_an_omega_section(user_nl): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:10:00"}} + + assert read_user_overrides(user_nl({"Omega": overrides})) == overrides + + +def test_unknown_user_overrides_are_rejected(user_nl): + overrides = {"TimeIntegration": {"TimeStepp": "0000_00:10:00"}} + + with pytest.raises(ValueError, match="Unknown option"): + read_user_overrides(user_nl(overrides)) + + +def test_blocked_user_overrides_are_rejected(user_nl): + overrides = {"IOStreams": {"RestartWrite": {"Precision": "single"}}} + + with pytest.raises(ValueError, match="cannot be overridden"): + read_user_overrides(user_nl(overrides)) + + +def test_user_overrides_that_are_not_a_mapping_are_rejected(user_nl): + with pytest.raises(ValueError, match="is not a mapping"): + read_user_overrides(user_nl(["TimeIntegration", "IOStreams"])) + + +def test_a_missing_user_nl_omega_is_reported(tmp_path): + with pytest.raises(FileNotFoundError, match="does not exist"): + read_user_overrides(tmp_path / "user_nl_omega") diff --git a/components/omega/cime_config/omega_buildnml/tests/test_validate_input_files.py b/components/omega/cime_config/omega_buildnml/tests/test_validate_input_files.py new file mode 100644 index 000000000000..d17d6c70a4e2 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/tests/test_validate_input_files.py @@ -0,0 +1,246 @@ +from copy import deepcopy + +import pytest + +from omega_buildnml.validate import validate_input_files_config + + +@pytest.fixture +def input_files(): + """A minimal, valid ``input_files.yaml`` configuration.""" + return { + "meshes": { + "Icos10": { + "inputs": [ + { + "file": "ocean.Icos10.nc", + "streams": [ + "HorzMeshIn", + "InitialVertCoord", + "InitialState", + ], + } + ] + } + } + } + + +@pytest.fixture +def defaults(): + """A minimal stand-in for the parsed ``Default.yml`` configuration.""" + return { + "IOStreams": { + "HorzMeshIn": {}, + "InitialVertCoord": {}, + "InitialState": {}, + "Forcing": {}, + } + } + + +@pytest.fixture +def mesh(input_files): + """The single mesh entry of the valid configuration.""" + return input_files["meshes"]["Icos10"] + + +@pytest.fixture +def input_group(mesh): + """The single input group of the valid configuration.""" + return mesh["inputs"][0] + + +def test_a_valid_configuration_is_returned(input_files, defaults): + assert validate_input_files_config(input_files, defaults) == input_files + + +def test_a_valid_mesh_is_returned(input_files, defaults): + validated = validate_input_files_config( + input_files, defaults, mesh_name="Icos10" + ) + + assert validated == input_files + + +@pytest.mark.parametrize("config", [{}, [], "meshes", None]) +def test_configurations_that_are_not_mappings_are_rejected(config, defaults): + with pytest.raises(ValueError, match="empty or is not a mapping"): + validate_input_files_config(config, defaults) + + +def test_unknown_top_level_keys_are_rejected(input_files, defaults): + input_files["mesh"] = {} + + with pytest.raises(ValueError, match=r"Unknown top-level key\(s\): mesh"): + validate_input_files_config(input_files, defaults) + + +@pytest.mark.parametrize("meshes", [{}, [], "Icos10", None]) +def test_meshes_that_are_not_mappings_are_rejected( + input_files, meshes, defaults +): + input_files["meshes"] = meshes + + with pytest.raises(ValueError, match="`meshes` is missing, empty"): + validate_input_files_config(input_files, defaults) + + +def test_an_unsupported_mesh_is_rejected(input_files, defaults): + with pytest.raises(ValueError, match="Unsupported OCN_GRID"): + validate_input_files_config(input_files, defaults, mesh_name="Icos30") + + +@pytest.mark.parametrize("entry", [{}, [], "inputs", None]) +def test_mesh_entries_that_are_not_mappings_are_rejected( + input_files, entry, defaults +): + input_files["meshes"]["Icos10"] = entry + + with pytest.raises(ValueError, match="is empty or is not a mapping"): + validate_input_files_config(input_files, defaults) + + +def test_unknown_keys_in_a_mesh_entry_are_rejected( + input_files, mesh, defaults +): + mesh["input"] = [] + + with pytest.raises(ValueError, match=r"Unknown key\(s\) input"): + validate_input_files_config(input_files, defaults) + + +@pytest.mark.parametrize("inputs", [{}, [], "file", None]) +def test_inputs_that_are_not_lists_are_rejected( + input_files, mesh, inputs, defaults +): + mesh["inputs"] = inputs + + with pytest.raises(ValueError, match="Missing inputs"): + validate_input_files_config(input_files, defaults) + + +def test_input_groups_that_are_not_mappings_are_rejected( + input_files, mesh, defaults +): + mesh["inputs"] = ["ocean.Icos10.nc"] + + with pytest.raises(ValueError, match="Input group 0 is not a mapping"): + validate_input_files_config(input_files, defaults) + + +def test_unknown_keys_in_an_input_group_are_rejected( + input_files, input_group, defaults +): + input_group["stream"] = [] + + with pytest.raises(ValueError, match=r"Unknown key\(s\) stream"): + validate_input_files_config(input_files, defaults) + + +@pytest.mark.parametrize("file_name", ["", [], None]) +def test_input_groups_without_a_file_are_rejected( + input_files, input_group, file_name, defaults +): + input_group["file"] = file_name + + with pytest.raises(ValueError, match="Missing file in input group 0"): + validate_input_files_config(input_files, defaults) + + +@pytest.mark.parametrize("streams", [[], {}, "InitialState", None]) +def test_input_groups_without_streams_are_rejected( + input_files, input_group, streams, defaults +): + input_group["streams"] = streams + + with pytest.raises(ValueError, match="Missing streams in input group 0"): + validate_input_files_config(input_files, defaults) + + +@pytest.mark.parametrize("stream", ["", [], None]) +def test_stream_names_that_are_not_strings_are_rejected( + input_files, input_group, stream, defaults +): + input_group["streams"].append(stream) + + with pytest.raises(ValueError, match="must be non-empty strings"): + validate_input_files_config(input_files, defaults) + + +def test_unknown_stream_names_are_rejected(input_files, input_group, defaults): + input_group["streams"].append("MyCustomStream") + + with pytest.raises(ValueError, match="Unknown IOStream 'MyCustomStream'"): + validate_input_files_config(input_files, defaults) + + +def test_streams_assigned_more_than_once_are_rejected( + input_files, mesh, defaults +): + mesh["inputs"].append({"file": "state.nc", "streams": ["InitialState"]}) + + with pytest.raises(ValueError, match="assigned more than once"): + validate_input_files_config(input_files, defaults) + + +def test_streams_may_be_split_across_input_groups(input_files, mesh, defaults): + mesh["inputs"] = [ + {"file": "mesh.nc", "streams": ["HorzMeshIn", "InitialVertCoord"]}, + {"file": "state.nc", "streams": ["InitialState"]}, + ] + + assert validate_input_files_config(input_files, defaults) == input_files + + +def test_missing_required_streams_are_rejected( + input_files, input_group, defaults +): + input_group["streams"].remove("InitialVertCoord") + + with pytest.raises( + ValueError, match=r"Missing required IOStream\(s\) InitialVertCoord" + ): + validate_input_files_config(input_files, defaults) + + +def test_optional_streams_may_be_assigned(input_files, mesh, defaults): + mesh["inputs"].append({"file": "forcing.nc", "streams": ["Forcing"]}) + + assert validate_input_files_config(input_files, defaults) == input_files + + +def test_every_mesh_is_validated_when_no_mesh_is_given(input_files, defaults): + broken = deepcopy(input_files["meshes"]["Icos10"]) + broken["inputs"][0]["streams"].remove("HorzMeshIn") + input_files["meshes"]["Icos30"] = broken + + with pytest.raises(ValueError, match="mesh: Icos30"): + validate_input_files_config(input_files, defaults) + + +def test_only_the_given_mesh_is_validated(input_files, defaults): + input_files["meshes"]["Icos30"] = {"inputs": []} + + validated = validate_input_files_config( + input_files, defaults, mesh_name="Icos10" + ) + + assert validated == input_files + + +def test_problems_are_reported_together( + input_files, mesh, input_group, defaults +): + mesh["input"] = [] + input_group["streams"].remove("InitialVertCoord") + input_group["streams"].append(1) + + with pytest.raises(ValueError) as error: + validate_input_files_config(input_files, defaults) + + reported = str(error.value) + + assert "Unknown key(s) input" in reported + assert "must be non-empty strings" in reported + assert "Missing required IOStream(s) InitialVertCoord" in reported diff --git a/components/omega/cime_config/omega_buildnml/tests/test_validate_overrides.py b/components/omega/cime_config/omega_buildnml/tests/test_validate_overrides.py new file mode 100644 index 000000000000..a605fa082449 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/tests/test_validate_overrides.py @@ -0,0 +1,280 @@ +import pytest + +from omega_buildnml.validate import ( + validate_blocked_options, + validate_config_overrides, + validate_overrides, +) + + +@pytest.fixture +def defaults(): + """A minimal stand in for Omega's default configuration.""" + return { + "TimeIntegration": { + "TimeStepper": "Forward-Backward", + "TimeStep": "0000_00:30:00", + }, + "Tendencies": {"SurfaceTracerRestoringEnable": False}, + "IOStreams": { + "InitialState": {"Filename": "ocean.nc"}, + "History": {"Freq": 1, "FreqUnits": "months"}, + }, + } + + +@pytest.fixture +def input_files(): + """An ``input_files.yaml`` defining the meshes Omega supports.""" + return {"meshes": {"Icos10": {}, "Icos30": {}}} + + +@pytest.fixture +def config_overrides(): + """A minimal, valid ``config_overrides.yaml`` configuration.""" + return { + "coupled": { + "TimeIntegration": {"TimeStepper": "RungeKutta4"}, + }, + "meshes": { + "Icos10": {"TimeIntegration": {"TimeStep": "0000_00:05:00"}}, + }, + } + + +def test_a_valid_configuration_is_returned( + config_overrides, input_files, defaults +): + validated = validate_config_overrides( + config_overrides, input_files, defaults + ) + + assert validated == config_overrides + + +@pytest.mark.parametrize("config", [{}, [], "coupled", None]) +def test_configurations_that_are_not_mappings_are_rejected( + config, input_files, defaults +): + with pytest.raises(ValueError, match="empty or is not a mapping"): + validate_config_overrides(config, input_files, defaults) + + +def test_unknown_top_level_keys_are_rejected( + config_overrides, input_files, defaults +): + config_overrides["mesh"] = {} + + with pytest.raises(ValueError, match=r"Unknown top-level key\(s\): mesh"): + validate_config_overrides(config_overrides, input_files, defaults) + + +@pytest.mark.parametrize("coupled", [{}, [], "TimeIntegration", None]) +def test_coupled_overrides_are_required( + config_overrides, input_files, defaults, coupled +): + config_overrides["coupled"] = coupled + + with pytest.raises(ValueError, match="`coupled` is missing, empty"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_mesh_overrides_are_optional(config_overrides, input_files, defaults): + del config_overrides["meshes"] + + validated = validate_config_overrides( + config_overrides, input_files, defaults + ) + + assert validated == config_overrides + + +@pytest.mark.parametrize("meshes", [[], "Icos10"]) +def test_meshes_that_are_not_mappings_are_rejected( + config_overrides, input_files, defaults, meshes +): + config_overrides["meshes"] = meshes + + with pytest.raises(ValueError, match="`meshes` is not a mapping"): + validate_config_overrides(config_overrides, input_files, defaults) + + +@pytest.mark.parametrize("overrides", [{}, [], "TimeIntegration", None]) +def test_mesh_entries_that_are_not_mappings_are_rejected( + config_overrides, input_files, defaults, overrides +): + config_overrides["meshes"]["Icos10"] = overrides + + with pytest.raises(ValueError, match="is empty or is not a mapping"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_overrides_for_unsupported_meshes_are_rejected( + config_overrides, input_files, defaults +): + config_overrides["meshes"]["Icos120"] = { + "TimeIntegration": {"TimeStep": "0000_01:00:00"} + } + + with pytest.raises(ValueError, match="Unsupported mesh: Icos120"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_unknown_coupled_overrides_are_rejected( + config_overrides, input_files, defaults +): + config_overrides["coupled"]["TimeIntegration"]["TimeSteper"] = "RK4" + + with pytest.raises(ValueError, match="in coupled overrides"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_unknown_mesh_overrides_are_rejected( + config_overrides, input_files, defaults +): + config_overrides["meshes"]["Icos10"]["TimeIntegration"]["Step"] = "0" + + with pytest.raises(ValueError, match="in overrides for mesh: Icos10"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_iostreams_under_a_mesh_entry_are_rejected( + config_overrides, input_files, defaults +): + config_overrides["meshes"]["Icos10"]["IOStreams"] = { + "MyStream": {"Freq": 1} + } + + with pytest.raises( + ValueError, match="`IOStreams` is not permitted under mesh: Icos10" + ): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_iostreams_are_allowed_under_coupled( + config_overrides, input_files, defaults +): + config_overrides["coupled"]["IOStreams"] = {"MyStream": {"Freq": 1}} + + validated = validate_config_overrides( + config_overrides, input_files, defaults + ) + + assert validated == config_overrides + + +def test_every_mesh_is_validated_when_no_mesh_is_given( + config_overrides, input_files, defaults +): + config_overrides["meshes"]["Icos30"] = {"TimeIntegration": {"Step": "0"}} + + with pytest.raises(ValueError, match="in overrides for mesh: Icos30"): + validate_config_overrides(config_overrides, input_files, defaults) + + +def test_only_the_given_mesh_is_validated( + config_overrides, input_files, defaults +): + config_overrides["meshes"]["Icos30"] = {"TimeIntegration": {"Step": "0"}} + + validated = validate_config_overrides( + config_overrides, input_files, defaults, mesh_name="Icos10" + ) + + assert validated == config_overrides + + +def test_problems_are_reported_together( + config_overrides, input_files, defaults +): + config_overrides["mesh"] = {} + config_overrides["coupled"]["Tendencies"] = {"Restoring": True} + config_overrides["meshes"]["Icos120"] = {"TimeIntegration": {}} + + with pytest.raises(ValueError) as error: + validate_config_overrides(config_overrides, input_files, defaults) + + reported = str(error.value) + + assert "Unknown top-level key(s): mesh" in reported + assert "Tendencies.Restoring" in reported + assert "Unsupported mesh: Icos120" in reported + + +def test_overrides_of_known_options_are_valid(defaults): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:05:00"}} + + assert validate_overrides(overrides, defaults, "test") == [] + + +def test_overrides_of_unknown_options_are_reported(defaults): + overrides = {"TimeIntegration": {"TimeSteps": "0000_00:05:00"}} + + assert validate_overrides(overrides, defaults, "test") == [ + "Unknown option(s) TimeIntegration.TimeSteps in test." + ] + + +def test_unknown_sections_are_reported(defaults): + overrides = {"Tendancies": {"SurfaceTracerRestoringEnable": True}} + + assert validate_overrides(overrides, defaults, "test") == [ + "Unknown option(s) Tendancies in test." + ] + + +def test_options_that_are_not_sections_are_reported(defaults): + overrides = {"TimeIntegration": "0000_00:05:00"} + + assert validate_overrides(overrides, defaults, "test") == [ + "Unknown option(s) TimeIntegration in test." + ] + + +def test_sections_that_are_not_options_are_reported(defaults): + overrides = {"TimeIntegration": {"TimeStep": {"Value": "0000_00:05:00"}}} + + assert validate_overrides(overrides, defaults, "test") == [ + "Unknown option(s) TimeIntegration.TimeStep in test." + ] + + +def test_new_streams_are_allowed(defaults): + overrides = {"IOStreams": {"MyStream": {"Freq": 1}}} + + assert validate_overrides(overrides, defaults, "test") == [] + + +def test_unknown_options_are_reported_together(defaults): + overrides = { + "TimeIntegration": {"TimeSteps": "0000_00:05:00"}, + "Tendancies": {"SurfaceTracerRestoringEnable": True}, + } + + assert validate_overrides(overrides, defaults, "test") == [ + "Unknown option(s) TimeIntegration.TimeSteps, Tendancies in test." + ] + + +def test_options_that_are_not_blocked_are_allowed(): + overrides = {"TimeIntegration": {"TimeStep": "0000_00:05:00"}} + + assert validate_blocked_options(overrides, "test") == [] + + +def test_blocked_options_are_reported(): + overrides = {"TimeIntegration": {"StartTime": "0001-01-01_00:00:00"}} + + assert validate_blocked_options(overrides, "test") == [ + "Option(s) TimeIntegration.StartTime in test are set by CIME and " + "cannot be overridden." + ] + + +def test_blocked_sections_are_reported_rather_than_their_options(): + overrides = {"IOStreams": {"RestartWrite": {"Precision": "single"}}} + + assert validate_blocked_options(overrides, "test") == [ + "Option(s) IOStreams.RestartWrite in test are set by CIME and " + "cannot be overridden." + ] diff --git a/components/omega/cime_config/omega_buildnml/tests/test_validate_user_overrides.py b/components/omega/cime_config/omega_buildnml/tests/test_validate_user_overrides.py new file mode 100644 index 000000000000..7ca59a65e545 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/tests/test_validate_user_overrides.py @@ -0,0 +1,87 @@ +import pytest + +from omega_buildnml.validate import validate_user_overrides + + +@pytest.fixture +def defaults(): + """A minimal stand in for Omega's default configuration.""" + return { + "TimeIntegration": { + "TimeStepper": "Forward-Backward", + "TimeStep": "0000_00:30:00", + "StartTime": "0001-01-01_00:00:00", + }, + "Tendencies": {"SurfaceTracerRestoringEnable": False}, + "IOStreams": { + "InitialState": {"Filename": "ocean.nc"}, + "History": {"Freq": 1, "FreqUnits": "months"}, + }, + } + + +def test_an_empty_configuration_is_returned(defaults): + assert validate_user_overrides({}, defaults) == {} + + +def test_a_valid_configuration_is_returned(defaults): + user_overrides = {"Tendencies": {"SurfaceTracerRestoringEnable": True}} + + validated = validate_user_overrides(user_overrides, defaults) + + assert validated == user_overrides + + +@pytest.mark.parametrize("config", [[], "TimeIntegration", None]) +def test_configurations_that_are_not_mappings_are_rejected(config, defaults): + with pytest.raises(ValueError, match="is not a mapping"): + validate_user_overrides(config, defaults) + + +def test_unknown_options_are_rejected(defaults): + user_overrides = {"TimeIntegration": {"TimeSteps": "0000_00:05:00"}} + + with pytest.raises(ValueError, match="Unknown option"): + validate_user_overrides(user_overrides, defaults) + + +def test_blocked_options_are_rejected(defaults): + user_overrides = {"TimeIntegration": {"StartTime": "0002-01-01_00:00:00"}} + + with pytest.raises(ValueError, match="cannot be overridden"): + validate_user_overrides(user_overrides, defaults) + + +def test_custom_iostreams_are_allowed(defaults): + """ + A user can add a wholly new IOStream in ``user_nl_omega``, so long as + the stream itself is a valid configuration (e.g. for debugging). + + Regression test for a bug where the now-removed ``KNOWN_STREAMS`` list + was mistakenly thought to gate this; ``IOStreams`` is an + ``OPEN_SECTIONS`` entry in :func:`validate_overrides`, so custom stream + names were never actually blocked here. + """ + user_overrides = { + "IOStreams": { + "MyCustomHiFreq": { + "Filename": "ocn.hifreq.$Y-$M", + "Mode": "write", + "Freq": 1, + "FreqUnits": "second", + "Contents": ["State"], + } + } + } + + validated = validate_user_overrides(user_overrides, defaults) + + assert validated == user_overrides + + +def test_overriding_an_existing_stream_is_allowed(defaults): + user_overrides = {"IOStreams": {"History": {"Freq": 5}}} + + validated = validate_user_overrides(user_overrides, defaults) + + assert validated == user_overrides diff --git a/components/omega/cime_config/omega_buildnml/validate.py b/components/omega/cime_config/omega_buildnml/validate.py new file mode 100644 index 000000000000..514c024c4170 --- /dev/null +++ b/components/omega/cime_config/omega_buildnml/validate.py @@ -0,0 +1,600 @@ +from typing import Optional + +from ._types import YamlMapping + +DATA_PATH = "components/omega/cime_config/omega_buildnml/data" + +INPUT_FILES_PATH = f"{DATA_PATH}/input_files.yaml" + +OVERRIDES_PATH = f"{DATA_PATH}/config_overrides.yaml" + +DEFAULTS_PATH = "components/omega/configs/Default.yml" + +#: IOStreams that every mesh must provide an input file for +REQUIRED_STREAMS = frozenset( + {"HorzMeshIn", "InitialVertCoord", "InitialState"} +) + +#: Keys allowed in a single mesh entry +MESH_KEYS = frozenset({"inputs"}) + +#: Keys allowed in a single entry of a mesh's ``inputs`` list +INPUT_GROUP_KEYS = frozenset({"file", "streams"}) + +#: Keys allowed at the top level of ``config_overrides.yaml`` +OVERRIDES_KEYS = frozenset({"coupled", "meshes"}) + +#: Sections not strictly validated against the defaults. IOStreams entries may +#: define new streams, and have requirements that depend on the values passed, +#: so they need their own validation. +OPEN_SECTIONS = frozenset({"IOStreams"}) + +#: IOStreams wholly controlled by CIME and the coupler. The required streams +#: are staged through ``omega.input_data_list``, the restart streams are named +#: and scheduled by the coupler, and forcing is provided by the coupler. +BLOCKED_STREAMS = frozenset( + REQUIRED_STREAMS | {"Forcing", "RestartRead", "RestartWrite"} +) + +#: Config options set by CIME, which a user is not permitted to override. +#: Matched as prefixes, so naming a section blocks everything below it. +BLOCKED_OPTIONS = frozenset( + {f"IOStreams.{stream}" for stream in BLOCKED_STREAMS} + | { + # start, stop, and duration are provided by the coupler at runtime + "TimeIntegration.StartTime", + "TimeIntegration.StopTime", + "TimeIntegration.RunDuration", + # calendar must agree with the CIME ``CALENDAR`` setting + "TimeIntegration.CalendarType", + } +) + + +def validate_input_files_config( + input_files: YamlMapping, + defaults: YamlMapping, + mesh_name: Optional[str] = None, +) -> YamlMapping: + """ + Validate the contents of the ``input_files.yaml`` configuration. + + All problems found are collected and reported together, rather than + raising on the first one encountered. + + Parameters: + ----------- + input_files : dict[str, Any] + Parsed content of ``cime_config/omega_buildnml/data/input_files.yaml`` + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. Used + to check that streams assigned an input file are streams Omega + actually defines, since ``input_files.yaml`` only supplies a + ``Filename`` for an existing ``IOStreams`` entry. + mesh_name : str, optional + The name of the mesh to validate. If not provided, all mesh entries + will be validated. + + Returns: + -------- + dict[str, Any] + The validated configuration. + + Raises: + ------- + ValueError + If any required keys are missing or if any values are invalid. + """ + if not isinstance(input_files, dict) or not input_files: + err_msg = f"`{INPUT_FILES_PATH}` is empty or is not a mapping." + raise ValueError(err_msg) + + unknown_keys = sorted(set(input_files) - {"meshes"}) + if unknown_keys: + _raise( + [f"Unknown top-level key(s): {', '.join(unknown_keys)}."], + INPUT_FILES_PATH, + ) + + meshes: YamlMapping = input_files.get("meshes", {}) + if not isinstance(meshes, dict) or not meshes: + _raise( + ["`meshes` is missing, empty, or is not a mapping."], + INPUT_FILES_PATH, + ) + + known_streams = frozenset(defaults.get("IOStreams", {})) + + if mesh_name is None: + errors = [] + for name in meshes: + errors.extend( + _validate_input_files_entry(input_files, known_streams, name) + ) + _raise(errors, INPUT_FILES_PATH) + return input_files + + if mesh_name not in meshes: + err_msg = ( + f"Unsupported OCN_GRID for Omega: {mesh_name}. \n" + f"Could not find entry in `{INPUT_FILES_PATH}`" + ) + raise ValueError(err_msg) + + _raise( + _validate_input_files_entry(input_files, known_streams, mesh_name), + INPUT_FILES_PATH, + ) + + return input_files + + +def validate_config_overrides( + config_overrides: YamlMapping, + input_files: YamlMapping, + defaults: YamlMapping, + mesh_name: Optional[str] = None, +) -> YamlMapping: + """ + Validate the contents of the ``config_overrides.yaml`` configuration. + + All problems found are collected and reported together, rather than + raising on the first one encountered. + + Mesh specific overrides are optional, so a mesh without an entry is not + an error. + + Parameters: + ----------- + config_overrides : dict[str, Any] + Parsed content of + ``cime_config/omega_buildnml/data/config_overrides.yaml`` + input_files : dict[str, Any] + Parsed content of ``cime_config/omega_buildnml/data/input_files.yaml``, + which defines the meshes Omega supports. + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + mesh_name : str, optional + The name of the mesh to validate. If not provided, all mesh entries + will be validated. + + Returns: + -------- + dict[str, Any] + The validated configuration. + + Raises: + ------- + ValueError + If any required keys are missing or if any values are invalid. + """ + if not isinstance(config_overrides, dict) or not config_overrides: + err_msg = f"`{OVERRIDES_PATH}` is empty or is not a mapping." + raise ValueError(err_msg) + + errors: list[str] = [] + + unknown_keys = sorted(set(config_overrides) - OVERRIDES_KEYS) + if unknown_keys: + errors.append(f"Unknown top-level key(s): {', '.join(unknown_keys)}.") + + coupled = config_overrides.get("coupled") + if not isinstance(coupled, dict) or not coupled: + errors.append("`coupled` is missing, empty, or is not a mapping.") + else: + errors.extend( + validate_overrides(coupled, defaults, "coupled overrides") + ) + + meshes: YamlMapping = config_overrides.get("meshes", {}) + if not isinstance(meshes, dict): + errors.append("`meshes` is not a mapping.") + _raise(errors, OVERRIDES_PATH) + + if mesh_name is None: + for name in meshes: + errors.extend( + _validate_config_overrides_entry( + config_overrides, input_files, defaults, name + ) + ) + elif mesh_name in meshes: + errors.extend( + _validate_config_overrides_entry( + config_overrides, input_files, defaults, mesh_name + ) + ) + + _raise(errors, OVERRIDES_PATH) + + return config_overrides + + +def validate_user_overrides( + user_overrides: YamlMapping, defaults: YamlMapping +) -> YamlMapping: + """ + Validate the contents of a case's ``user_nl_omega`` file. + + All problems found are collected and reported together, rather than + raising on the first one encountered. + + An empty ``user_nl_omega`` is not an error, it just means the user has + not overridden anything. + + Parameters: + ----------- + user_overrides : dict[str, Any] + Parsed content of the case's ``user_nl_omega``. + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + + Returns: + -------- + dict[str, Any] + The validated user overrides. + + Raises: + ------- + ValueError + If any unknown or blocked options are set. + """ + if not isinstance(user_overrides, dict): + err_msg = "`user_nl_omega` is not a mapping." + raise ValueError(err_msg) + + errors = validate_overrides(user_overrides, defaults, "user overrides") + errors.extend(validate_blocked_options(user_overrides, "user overrides")) + + _raise(errors, "user_nl_omega") + + return user_overrides + + +def _raise(errors: list[str], config_path: str) -> None: + """ + Raise a single ``ValueError`` describing all accumulated errors. + + Does nothing when ``errors`` is empty. + + Parameters: + ----------- + errors : list[str] + Error messages collected during validation. + config_path : str + Path of the configuration file the errors were found in. + + Raises: + ------- + ValueError + If ``errors`` is non-empty. + """ + if not errors: + return + + details = "\n".join(f" - {error}" for error in errors) + err_msg = ( + f"Invalid Omega configuration:\n{details}\n" + f"Please check your setting in `{config_path}`" + ) + raise ValueError(err_msg) + + +def _validate_input_files_entry( + input_files: YamlMapping, known_streams: frozenset, mesh_name: str +) -> list[str]: + """ + Validate that the specified mesh has a valid configuration in input_files. + + Parameters: + ----------- + input_files : dict[str, Any] + Parsed content of ``cime_config/omega_buildnml/data/input_files.yaml`` + known_streams : frozenset[str] + Names of the ``IOStreams`` Omega defines in ``configs/Default.yml``. + mesh_name : str + The name of the mesh to validate. + + Returns: + -------- + list[str] + Error messages describing any problems found. Empty when the entry is + valid. + """ + meshes: YamlMapping = input_files["meshes"] + mesh: YamlMapping = meshes[mesh_name] + + if not isinstance(mesh, dict) or not mesh: + return [f"Mesh: {mesh_name} is empty or is not a mapping."] + + errors: list[str] = [] + streams_files = {} + + unknown_keys = sorted(set(mesh) - MESH_KEYS) + if unknown_keys: + errors.append( + f"Unknown key(s) {', '.join(unknown_keys)} for mesh: {mesh_name}." + ) + + inputs = mesh.get("inputs") + if not isinstance(inputs, list) or not inputs: + errors.append( + f"Missing inputs, or inputs is not a non-empty list, for " + f"mesh: {mesh_name}." + ) + return errors + + missing_msg = "Missing {key} in input group {index} for mesh: {mesh_name}." + + for index, input_group in enumerate(inputs): + if not isinstance(input_group, dict): + errors.append( + f"Input group {index} is not a mapping for mesh: {mesh_name}." + ) + continue + + unknown_keys = sorted(set(input_group) - INPUT_GROUP_KEYS) + if unknown_keys: + errors.append( + f"Unknown key(s) {', '.join(unknown_keys)} in input group " + f"{index} for mesh: {mesh_name}." + ) + + file_name = input_group.get("file") + streams = input_group.get("streams") + + if not isinstance(file_name, str) or not file_name: + errors.append( + missing_msg.format( + key="file", index=index, mesh_name=mesh_name + ) + ) + + if not isinstance(streams, list) or not streams: + errors.append( + missing_msg.format( + key="streams", index=index, mesh_name=mesh_name + ) + ) + continue + + for stream in streams: + if not isinstance(stream, str) or not stream: + errors.append( + f"Stream names must be non-empty strings, got " + f"'{stream}' in input group {index} for " + f"mesh: {mesh_name}." + ) + continue + + if stream in streams_files: + errors.append( + f"Stream '{stream}' is assigned more than once for " + f"mesh: {mesh_name}." + ) + continue + + if stream not in known_streams: + errors.append( + f"Unknown IOStream '{stream}' in input group {index} " + f"for mesh: {mesh_name}. Streams referenced in " + f"input_files.yaml must already be defined in " + f"Default.yml." + ) + continue + + # just store filename; the full path will be resolved later + # must store something to test for duplicates + streams_files[stream] = str(file_name) + + missing_streams = REQUIRED_STREAMS - set(streams_files) + if missing_streams: + errors.append( + f"Missing required IOStream(s) " + f"{', '.join(sorted(missing_streams))} for mesh: {mesh_name}." + ) + + return errors + + +def _validate_config_overrides_entry( + config_overrides: YamlMapping, + input_files: YamlMapping, + defaults: YamlMapping, + mesh_name: str, +) -> list[str]: + """ + Validate the overrides of a single mesh in config_overrides. + + Parameters: + ----------- + config_overrides : dict[str, Any] + Parsed content of + ``cime_config/omega_buildnml/data/config_overrides.yaml`` + input_files : dict[str, Any] + Parsed content of ``cime_config/omega_buildnml/data/input_files.yaml``, + which defines the meshes Omega supports. + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + mesh_name : str + The name of the mesh to validate. + + Returns: + -------- + list[str] + Error messages describing any problems found. Empty when the entry is + valid. + """ + meshes: YamlMapping = config_overrides["meshes"] + supported_meshes: YamlMapping = input_files.get("meshes", {}) + overrides: YamlMapping = meshes[mesh_name] + + errors: list[str] = [] + + if not isinstance(overrides, dict) or not overrides: + errors.append(f"Mesh: {mesh_name} is empty or is not a mapping.") + else: + errors.extend( + validate_overrides( + overrides, defaults, f"overrides for mesh: {mesh_name}" + ) + ) + if "IOStreams" in overrides: + errors.append( + f"`IOStreams` is not permitted under mesh: {mesh_name}. " + f"IOStreams shared by every mesh belong under `coupled`, " + f"and case-specific IOStreams belong in `user_nl_omega`." + ) + + if mesh_name not in supported_meshes: + errors.append( + f"Unsupported mesh: {mesh_name}. Could not find entry in " + f"`{INPUT_FILES_PATH}`." + ) + + return errors + + +def validate_overrides( + overrides: YamlMapping, defaults: YamlMapping, source: str +) -> list[str]: + """ + Validate that overrides only set options defined in Omega's defaults. + + Overrides that do not appear in the defaults would silently add new + options, rather than overriding an existing one. Sections listed in + ``OPEN_SECTIONS`` are not checked. + + Parameters: + ----------- + overrides : dict[str, Any] + Override options to validate. + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + source : str + Description of where the overrides came from, used in error messages. + + Returns: + -------- + list[str] + Error messages describing any problems found. Empty when the overrides + are valid. + """ + unknown_options = _unknown_override_options(overrides, defaults) + + if not unknown_options: + return [] + + return [f"Unknown option(s) {', '.join(unknown_options)} in {source}."] + + +def validate_blocked_options(overrides: YamlMapping, source: str) -> list[str]: + """ + Validate that overrides do not set options controlled by CIME. + + Options listed in ``BLOCKED_OPTIONS`` are set from the case configuration, + or by the coupler at runtime, so overriding them would either be silently + discarded or leave the run inconsistent with the rest of the case. + + Parameters: + ----------- + overrides : dict[str, Any] + Override options to validate. + source : str + Description of where the overrides came from, used in error messages. + + Returns: + -------- + list[str] + Error messages describing any problems found. Empty when the overrides + are valid. + """ + blocked_options = _blocked_override_options(overrides) + + if not blocked_options: + return [] + + return [ + f"Option(s) {', '.join(blocked_options)} in {source} are set by CIME " + f"and cannot be overridden." + ] + + +def _unknown_override_options( + overrides: YamlMapping, defaults: YamlMapping, prefix: str = "" +) -> list[str]: + """ + Find override options that are not defined in Omega's defaults. + + Parameters: + ----------- + overrides : dict[str, Any] + Override options to validate. + defaults : dict[str, Any] + Default configuration values, loaded from configs/Default.yml. + prefix : str, optional + Dotted path of the parent section, used when recursing. + + Returns: + -------- + list[str] + Dotted paths of any options not defined in the defaults. + """ + unknown_options: list[str] = [] + + for key, value in overrides.items(): + option = f"{prefix}{key}" + + if option in OPEN_SECTIONS: + continue + + if key not in defaults: + unknown_options.append(option) + continue + + default_value = defaults[key] + + if isinstance(value, dict) != isinstance(default_value, dict): + unknown_options.append(option) + elif isinstance(value, dict): + unknown_options.extend( + _unknown_override_options(value, default_value, f"{option}.") + ) + + return unknown_options + + +def _blocked_override_options( + overrides: YamlMapping, prefix: str = "" +) -> list[str]: + """ + Find override options that are controlled by CIME. + + Every level of the override tree is checked on the way down, so listing a + section in ``BLOCKED_OPTIONS`` blocks all the options below it. + + Parameters: + ----------- + overrides : dict[str, Any] + Override options to validate. + prefix : str, optional + Dotted path of the parent section, used when recursing. + + Returns: + -------- + list[str] + Dotted paths of any options that are controlled by CIME. + """ + blocked_options: list[str] = [] + + for key, value in overrides.items(): + option = f"{prefix}{key}" + + if option in BLOCKED_OPTIONS: + blocked_options.append(option) + elif isinstance(value, dict): + blocked_options.extend( + _blocked_override_options(value, f"{option}.") + ) + + return blocked_options diff --git a/components/omega/cime_config/user_nl_omega b/components/omega/cime_config/user_nl_omega index c5186fe22b7f..974e00ac7e3e 100644 --- a/components/omega/cime_config/user_nl_omega +++ b/components/omega/cime_config/user_nl_omega @@ -1,11 +1,24 @@ -!---------------------------------------------------------------------------------- -! Users should add all user specific namelist changes after these comments -! in the form of -! namelist_var = new_namelist_value -! *** EXCEPT FOR *** -! TODO: Change to OMEGA config names -! 1. DO NOT CHANGE config_start_time, config_run_duration, config_stop_time, -! config_do_restart, config_Restart_timestamp_filename, config_calendar_type -! 2. To preview the namelists, invoke $CASEROOT preview-namelists and look at -! $CASEROOT/CaseDocs/omega_in -!---------------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# Users should add all user specific configuration changes after these +# comments, as a YAML snippet mirroring the structure of Omega's default +# configuration (components/omega/configs/Default.yml). Only the options +# being changed need to be listed. For example: +# +# Omega: +# TimeIntegration: +# TimeStep: "0000_00:10:00" +# IOStreams: +# History: +# Freq: 1 +# FreqUnits: months +# +# *** THE FOLLOWING CANNOT BE CHANGED *** +# They are set from the case configuration, or by the coupler at runtime: +# 1. TimeIntegration: StartTime, StopTime, RunDuration, CalendarType +# 2. IOStreams: HorzMeshIn, InitialVertCoord, InitialState, Forcing, +# RestartRead, and RestartWrite. The History and Highfreq streams may +# be changed, as may any new stream you define. +# +# To preview the configuration, invoke $CASEROOT/preview_namelists and look +# at $CASEROOT/CaseDocs/omega.yml +# ---------------------------------------------------------------------------- diff --git a/components/omega/cime_config/validate_config.py b/components/omega/cime_config/validate_config.py new file mode 100755 index 000000000000..ae45c3aeac58 --- /dev/null +++ b/components/omega/cime_config/validate_config.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 + +""" +Validate the configuration files used to build Omega's runtime config + +Every entry in the packaged configuration files is validated, rather than only +the entries needed by a single mesh. This ensures the configuration files are +valid and complete, when called from ./case.setup. + +The following files are validated: + + components/omega/cime_config/omega_buildnml/data/input_files.yaml + Checks that every mesh assigns an input file to each of the required + IOStreams (HorzMeshIn, InitialVertCoord, InitialState), and that no + IOStream is assigned more than one file. Every stream name must already + be defined in configs/Default.yml, since input_files.yaml only supplies + a Filename for an existing IOStreams entry. + + components/omega/cime_config/omega_buildnml/data/config_overrides.yaml + Checks that the coupled overrides are present, and that every mesh with + overrides is a mesh defined in input_files.yaml. The coupled and mesh + specific options are checked against configs/Default.yml, so that an + override cannot silently add a new option, rather than setting an existing + one. The coupled section's IOStreams overrides are not checked against + Default.yml, since they are allowed to define new streams; per-mesh + IOStreams overrides are not permitted at all. + +Exits non-zero, and reports all the problems found, if any file is invalid. +""" + +import argparse + +from omega_buildnml import read_config_overrides, read_input_files_config + + +def main() -> None: + """ + Validate all of Omega's packaged configuration files. + """ + parser = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + parser.parse_args() + + # validate all the mesh entries in input_files.yaml + read_input_files_config() + + # validate the coupled and all the mesh entries in config_overrides.yaml + read_config_overrides() + + print("PASS: Omega configuration files are valid") + + +if __name__ == "__main__": + main() diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..e33403c7cc91 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -58,6 +58,8 @@ Omega: Mode: Implicit Type: Constant BottomDragCoeff: 1.0e-3 + SfcThicknessForcingTendencyEnable: false + SfcTracerForcingTendencyEnable: false TracerHorzAdvTendencyEnable: true TracerDiffTendencyEnable: true EddyDiff2: 10.0 diff --git a/components/omega/dev-conda.txt b/components/omega/dev-conda.txt index 525462603405..f73607572fa3 100644 --- a/components/omega/dev-conda.txt +++ b/components/omega/dev-conda.txt @@ -23,6 +23,10 @@ ruff flynt mypy +# omega_buildnml testing +pytest +pyyaml + # documentation sphinx sphinx_rtd_theme diff --git a/components/omega/doc/devGuide/BuildNml.md b/components/omega/doc/devGuide/BuildNml.md new file mode 100644 index 000000000000..b765cb99f047 --- /dev/null +++ b/components/omega/doc/devGuide/BuildNml.md @@ -0,0 +1,69 @@ +(omega-dev-buildnml)= + +# CIME `buildnml` and Configuration Validation + +`components/omega/cime_config/buildnml` generates a case's `omega.yml` by +layering Omega's defaults, coupled/mesh overrides, runtime overrides derived +from the case, and `user_nl_omega` (see {ref}`omega-user-config`). This +page covers the Python package behind `buildnml`, how its configuration is +validated, and how to add support for a new mesh. + +## `omega_buildnml` package + +`components/omega/cime_config/omega_buildnml/` implements the merging and +validation logic used by `buildnml`: + +- `read_write.py` reads/writes YAML config files and packaged data. +- `config.py` merges the configuration layers and resolves mesh input files. +- `validate.py` validates overrides and `IOStreams` against `Default.yml`. +- `data/input_files.yaml` maps each supported mesh to its input file(s). +- `data/config_overrides.yaml` holds the coupled and mesh-specific overrides. + +Validation runs whenever these files are read, so a bad edit fails fast at +`case.setup` rather than surfacing as a confusing runtime error. Stream +names are checked dynamically against `Default.yml` rather than a +hardcoded list: the `coupled` section of `config_overrides.yaml` and a +case's `user_nl_omega` are free to define brand-new `IOStreams` entries, +but a stream referenced in `input_files.yaml` must already exist in +`Default.yml`, since `input_files.yaml` only supplies a `Filename` +override for an existing `IOStreams` entry. `IOStreams` is not permitted +under a mesh entry in `config_overrides.yaml`'s `meshes` section, since +per-mesh IOStreams overrides aren't a supported use case; put IOStreams +that apply to every mesh under `coupled` instead. The required streams +(`HorzMeshIn`, `InitialVertCoord`, `InitialState`) are always enforced. + +## Validation and CI + +`components/omega/cime_config/validate_config.py` runs the same validation +against every mesh and override entry (rather than just the ones needed for +one case), so it catches problems anywhere in the packaged configuration: + +```bash +cd components/omega/cime_config +./validate_config.py +``` + +The `omega-buildnml` GitHub Actions workflow runs this script, along with +`omega_buildnml`'s pytest unit tests across supported Python versions, on +every pull request touching `cime_config/` or `Default.yml`. Run both +locally before opening a PR that changes either: + +```bash +cd components/omega/cime_config +./validate_config.py +python -m pytest omega_buildnml/tests -v +``` + +## Adding a supported mesh + +1. Confirm the mesh's grid alias is already defined for E3SM in + `cime_config/config_grids.xml` at the repository root; `buildnml` looks + up input files by the case's `OCN_GRID` value. +2. Add an entry for the mesh to `data/input_files.yaml`, listing the input + file(s) that provide its `HorzMeshIn`, `InitialVertCoord`, and + `InitialState` streams (split across multiple `inputs` entries if the + initial condition is a separate file from the mesh). +3. If the mesh needs overrides beyond Omega's defaults and the coupled + overrides, add a `meshes.` entry to `data/config_overrides.yaml`. +4. Run `./validate_config.py` to confirm the new entries are complete and + consistent with `Default.yml`. diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 97d38d8ae42d..7e69dc302bc6 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -6,7 +6,8 @@ This page describes design and implementation details for forcing-related pathways in Omega, currently this includes: - Surface stress forcing (e.g. wind stress) -- Surface tracer restoring +- Surface thickness and tracer flux forcing (actively coupled or data-forced) +- Surface tracer restoring (soon to be ported as a field originating from the coupler) ## Surface stress forcing design @@ -37,6 +38,70 @@ pathways in Omega, currently this includes: - `Omega.Tendencies.SfcStressForcingTendencyEnable` - gates execution of surface stress forcing tendency kernel +## Surface thickness and tracer flux forcing design + +### Surface thickness and tracer flux forcing data flow + +**Thickness equation pathway:** + +1. External fields provide freshwater and salt flux components: + - `SnowFlux`, `RainFlux`, `EvaporationFlux` + - `SeaIceFreshWaterFlux`, `IceRunoffFlux`, `RiverRunoffFlux` + - `SeaIceSaltFlux` +2. `Forcing` stores the flux fields in `TracerForcingVars` +3. The tendency term `SfcThicknessForcingOnCell` sums the freshwater and salt mass fluxes and applies them to +the surface layer pseudo-thickness. + +**Tracer equation pathway:** + +1. External fields provide heat and salt flux components: + - `LatentHeatFlux`, `SensibleHeatFlux` + - `LongWaveHeatFluxUp`, `LongWaveHeatFluxDown` + - `SeaIceHeatFlux`, `ShortWaveHeatFlux` + - mass fluxes which add energy changes (`SnowFlux`, `RainFlux`, `IceRunoffFlux`, `RiverRunoffFlux`) + - `SeaIceSaltFlux` +2. `Forcing` stores the flux fields in `TracerForcingVars` +3. The tendency term `SfcTracerForcingOnCell` converts the summed external heat fluxes to a conservative-temperature tendency, + and applies the external sea-ice salt flux to the top layer salt content thus impacting salinity. + +### Surface thickness and tracer flux forcing key classes/components + +- `TracerForcingVars` + - Stores 13 coupled flux cell-centered fields: 6 freshwater fluxes, 6 heat + fluxes, and 1 salt flux component + - Fields initialized to zero and registered in `Forcing` field group +- `SfcThicknessForcingOnCell` tendency term + - Computes the layer mass contribution (converted to pseudo-thickness): $\sum (\text{SnowFlux} + \text{RainFlux} + \text{EvaporationFlux} + \text{SeaIceFreshWaterFlux} + \text{IceRunoffFlux} + \text{RiverRunoffFlux} + \text{SeaIceSaltFlux}) / \rho_{sw}$ + - Applied only at surface layer (top active layer) using `MinLayerCell` +- `SfcTracerForcingOnCell` tendency term + - For temperature: adds the direct heat fluxes + $Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{lw,up}} + Q_{\text{lw,down}} + Q_{\text{ice}} + Q_{\text{sw}}$ +, the phase change and enthalpy of added mass $(\text{RainFlux} + \text{RiverRunoffFlux}) c^0_{p,sw} C_T^{\text{top}} + (\text{SnowFlux} + \text{IceRunoffFlux})(c^0_{p,sw} C_T^{\text{frz}} - L_{\text{ice}})$, + (where $C_T^{\text{frz}}$ is from EOS at top-layer salinity and pressure), + and scales by $H_{\text{FluxFac}}$. + - For salinity: applies salt flux with unit conversion: $\text{SeaIceSaltFlux} \times S_{\text{FluxFac}}$ + - Applied only at surface layer using `MinLayerCell` + - Uses tracer index validation to apply to specific tracers only +- `Forcing` + - Manages `TracerForcingVars` instance +- `Tendencies` + - Calls `SfcThicknessForcingOnCell` in `computePseudoThicknessTendenciesOnly` + - Calls `SfcTracerForcingOnCell` in `computeTracerTendenciesOnly` after surface tracer restoring + +### Surface thickness and tracer flux forcing config coupling + +- `Omega.Tendencies.SfcThicknessForcingTendencyEnable` + - gates execution of coupled flux thickness kernel + - controls freshwater and salt flux forcing on sea surface height +- `Omega.Tendencies.SfcTracerForcingTendencyEnable` + - gates execution of coupled flux tracer kernel + - controls direct heat flux forcing on temperature and salt flux forcing on salinity + +## Notes + +- Currently all forcing is applied to the surface layer only. In the future, vertical spreading of river runoff contributions will be needed. +- `SeaIceFreshWaterFlux` is the pure freshwater mass from sea ice. The full mass flux from sea ice is `SeaIceFreshWaterFlux + SeaIceSaltFlux` + ## Surface tracer restoring design ### Surface tracer restoring data flow diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 5fa72197132f..2fd762e28b3e 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -41,9 +41,11 @@ implemented: - `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` +- `SfcThicknessForcingOnCell` +- `SfcTracerForcingOnCell` - `SurfaceTracerRestoringOnCell` ## See Also -Additional information on forcing (currently wind forcing and surface tracer -restoring) is detailed in [](omega-dev-forcing). +Additional information on forcing (surface stress, surface mass and tracer flux forcing, and +surface tracer restoring) is detailed in [](omega-dev-forcing). diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index dfc81a284670..76d45b9803be 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -76,6 +76,7 @@ devGuide/Driver devGuide/EOS devGuide/Broadcast devGuide/CMakeBuild +devGuide/BuildNml devGuide/Logging devGuide/Decomp devGuide/Dimension diff --git a/components/omega/doc/userGuide/Config.md b/components/omega/doc/userGuide/Config.md index f0fdbf60f49c..c00c251d25e1 100644 --- a/components/omega/doc/userGuide/Config.md +++ b/components/omega/doc/userGuide/Config.md @@ -88,3 +88,62 @@ input YAML-formatted file. Details of the implementation within omega can be found in the [Developer's Guide](#omega-dev-config) and the actual interfaces for extracting configuration variables into the modules that own them are described there as well. + +### Coupled Run Configuration (`user_nl_omega`) + +When Omega runs as the `ocn` component of a coupled E3SM case, `omega.yml` +is generated by `case.setup`/`case.submit` rather than written by hand, by +layering configuration from lowest to highest precedence: + +| Precedence | Layer | Source | +| :--------: | ----- | ------ | +| 1 (lowest) | Omega's defaults | `Default.yml` | +| 2 | coupled overrides | applied to every coupled run | +| 3 | mesh overrides | specific to the case's ocean mesh (`OCN_GRID`) | +| 4 | runtime overrides | derived from the case | +| 5 (highest) | user overrides | `user_nl_omega` | + +Runtime overrides (4) are derived from the case configuration: start/stop time, +calendar, restart behavior, and the resolved input file names for the case's +mesh. + +Each layer is merged on top of the ones below it, so `user_nl_omega` is the +way to customize an individual case without editing Omega's packaged +configuration. + +Edit `user_nl_omega` in the case directory before `case.setup`, wrapping +your overrides in a top-level `Omega:` key that matches `Default.yml`, e.g.: + +```yaml +Omega: + Tendencies: + SurfaceTracerRestoringEnable: true + + IOStreams: + History: + Freq: 5 + FreqUnits: Days +``` + +Unknown sections/options are rejected, as are typos. A small set of options +are set from the case configuration or by the coupler at runtime, and +cannot be overridden in `user_nl_omega`: + +```yaml +Omega: + TimeIntegration: + StartTime: ... # blocked: set from the case's run start date + StopTime: ... # blocked: controlled by the coupler + RunDuration: ... # blocked: controlled by the coupler + CalendarType: ... # blocked: must match the case's CALENDAR + + IOStreams: + HorzMeshIn: {} # blocked: mesh file, resolved from OCN_GRID + InitialVertCoord: {} # blocked: mesh file, resolved from OCN_GRID + InitialState: {} # blocked: initial condition, resolved from OCN_GRID + RestartRead: {} # blocked: managed by the coupler + RestartWrite: {} # blocked: managed by the coupler + Forcing: {} # blocked: provided by the coupler +``` + +Attempting to set one of these is a `case.setup` error. diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index d8fac4383730..01eeff4b5e3d 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -5,6 +5,7 @@ This page documents the user-facing configuration and behavior for current forcing in Omega: - Surface stress forcing (e.g. wind stress) +- Coupled tracer flux forcing (mass, energy and salt) - Surface tracer restoring ## Surface stress forcing @@ -39,6 +40,74 @@ Surface stress forcing uses surface stress input fields: These are stored in forcing variables and used to form edge-normal stress (`NormalStressEdge`) that enters momentum tendencies. +## Surface thickness and tracer flux forcing + +Surface thickness and tracer flux forcing applies ocean-atmosphere and ocean-sea ice fluxes from the other model +components (atmosphere, sea ice) to the thickness and tracer equations. This enables +the ocean to respond to heat, freshwater, and salt exchanges at the surface. These fluxes can be from data or (active) coupled components. + +### Surface thickness and tracer flux forcing configuration + +Surface thickness and tracer flux forcing is controlled by two configuration flags: + +```yaml +Omega: + Tendencies: + SfcThicknessForcingTendencyEnable: false + SfcTracerForcingTendencyEnable: false +``` + +- `Tendencies.SfcThicknessForcingTendencyEnable`: enables coupled freshwater and salt flux forcing on thickness +- `Tendencies.SfcTracerForcingTendencyEnable`: enables coupled heat and salt flux forcing on tracers + + +### Required input fields + +Coupled tracer flux forcing uses 13 auxiliary fields organized by type: + +**Freshwater mass fluxes (kg m⁻² s⁻¹):** +- `SnowFlux`: precipitation from snow +- `RainFlux`: precipitation from rain +- `EvaporationFlux`: evaporative water loss +- `SeaIceFreshWaterFlux`: freshwater mass flux from sea-ice melt or formation +- `IceRunoffFlux`: runoff from land ice +- `RiverRunoffFlux`: runoff from rivers + +**Heat fluxes (W m⁻²):** +- `LatentHeatFlux`: latent heat transfer +- `SensibleHeatFlux`: sensible heat transfer +- `LongWaveHeatFluxUp`: upward longwave radiation +- `LongWaveHeatFluxDown`: downward longwave radiation +- `SeaIceHeatFlux`: heat/energy from sea-ice interaction (incl. enthalpy of meltwater) +- `ShortWaveHeatFlux`: shortwave (solar) radiation + +**Salt mass flux (kg m⁻² s⁻¹):** +- `SeaIceSaltFlux`: salt flux from sea-ice formation/melt processes + +These fields are populated by external coupling components (typically atmosphere +and ice models). Omega assumes the incoming values match the documented units. +For now, there are assumed to come from a `forcing.nc` file, but later will be provided +by the equivalent `ocn_comp_mct.F`. + +### Notes + +- Coupled fluxes are applied only at the surface layer (top active layer) for each cell. In the future, vertical spreading of contributions from river runoff will be needed. +- Pseudo-thickness tendency is computed from the (six) freshwater mass fluxes and the salt mass flux + `SeaIceSaltFlux`, converted to a pseudo-thickness change. +- Temperature tendency is computed from direct heat flux plus + mass-flux enthalpy terms, converted to conservative-temperature tendency via + $H_{\text{FluxFac}} = 1.0 / (\rho_{sw} c^0_{p,sw})$ where $c^0_{p,sw}$ is the reference + specific heat of seawater defined by TEOS-10. The enthalpy associated with mass fluxes is currently hard-coded to SST for liquid fluxes and the freezing temperature for solid fluxes (which are melted using a constant latent heat of fusion). Note that the enthalpy of liquid meltwater from sea ice is already included in `SeaIceHeatFlux`. +- Salinity tendency from `SeaIceSaltFlux` is scaled by + $S_{\text{FluxFac}} = 1.0e3 / \rho_{sw}$ to account for unit conversion from + kg/(m²·s) to salinity units (g/kg). +- Fluxes are assumed to be in the documented units (i.e. net mass fluxes); + any unit conversion should be performed by the coupling component before providing flux + values to Omega. +- The reference density used here ($\rho_{sw}$) is not a Boussinesq density, it is the + conversion factor from mass to pseudo-thickness. +- No iceberg fluxes are included for now. + ## Surface tracer restoring Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 3db2b4c10098..c2cd4b3b4327 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -20,6 +20,8 @@ tendency terms are currently implemented: | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers | SfcStressForcingOnEdge | forcing by surface stress (e.g. wind), defined on edges | BottomDragOnEdge | bottom drag, defined on edges +| SfcThicknessForcingOnCell | surface pseudo-thickness forcing from coupled freshwater and salt fluxes, defined on cells +| SfcTracerForcingOnCell | surface tracer forcing from coupled heat and salt fluxes, with direct heat always and mass-flux enthalpy terms gated by thickness forcing, defined on cells | SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells Among the internal data stored by each functor is a `bool` which can enable or @@ -57,6 +59,8 @@ the currently available tendency terms: | | BottomDragTendency:Mode | bottom drag mode; `Implicit` or `Explicit` | | BottomDragTendency:Type | bottom drag type; `Constant` | | BottomDragTendency:BottomDragCoeff | bottom drag coefficient +| SfcThicknessForcingOnCell | SfcThicknessForcingTendencyEnable | enable/disable term +| SfcTracerForcingOnCell | SfcTracerForcingTendencyEnable | enable/disable term | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term ## Second Order Horizontal Advection Algorithm @@ -142,5 +146,6 @@ Tracer higer order convergence example of a cosine bell advected on a sphere sho ## See Also -Additional information on forcing (currently wind forcing and surface tracer -restoring) is detailed in [](omega-user-forcing). +Additional information on forcing, including surface stress forcing, +surface thickness and tracer flux forcing, and surface tracer restoring, is detailed in +[](omega-user-forcing). diff --git a/components/omega/src/CMakeLists.txt b/components/omega/src/CMakeLists.txt index dc7e764da3ab..a2f1ae537f58 100644 --- a/components/omega/src/CMakeLists.txt +++ b/components/omega/src/CMakeLists.txt @@ -52,9 +52,9 @@ target_link_libraries( OmegaLibFlags INTERFACE Kokkos::kokkos - spdlog + spdlog::spdlog pioc - yaml-cpp + yaml-cpp::yaml-cpp parmetis metis gptl diff --git a/components/omega/src/drivers/coupled/omega_cpl_indices.F90 b/components/omega/src/drivers/coupled/omega_cpl_indices.F90 index 0b22be5918c6..d05da77f49e6 100644 --- a/components/omega/src/drivers/coupled/omega_cpl_indices.F90 +++ b/components/omega/src/drivers/coupled/omega_cpl_indices.F90 @@ -5,8 +5,8 @@ module omega_cpl_indices implicit none private - integer, parameter, public :: num_omega_imports = 2 - integer, parameter, public :: num_omega_exports = 5 + integer, parameter, public :: num_omega_imports = 15 + integer, parameter, public :: num_omega_exports = 11 integer, public :: num_coupler_imports, num_coupler_exports ! Names of import/export fields as defined by seq_flds_mod @@ -43,6 +43,19 @@ subroutine omega_set_cpl_indices() ! Import (x2o) Coupler field names import_field_names(1) = "Foxx_taux" import_field_names(2) = "Foxx_tauy" + import_field_names(3) = "Foxx_swnet" + import_field_names(4) = "Foxx_sen" + import_field_names(5) = "Foxx_lat" + import_field_names(6) = "Foxx_lwup" + import_field_names(7) = "Faxa_lwdn" + import_field_names(8) = "Fioi_salt" + import_field_names(9) = "Fioi_melth" + import_field_names(10) = "Fioi_meltw" + import_field_names(11) = "Faxa_snow" + import_field_names(12) = "Faxa_rain" + import_field_names(13) = "Foxx_evap" + import_field_names(14) = "Foxx_rofl" + import_field_names(15) = "Foxx_rofi" ! get mct_avect_index value for each import field name call get_indices_from_names( & @@ -58,6 +71,12 @@ subroutine omega_set_cpl_indices() export_field_names(3) = "So_u" export_field_names(4) = "So_v" export_field_names(5) = "So_ssh" + export_field_names(6) = "So_dhdx" + export_field_names(7) = "So_dhdy" + export_field_names(8) = "So_fswpen" + export_field_names(9) = "Faoo_h2otemp" + export_field_names(10) = "Fioo_q" + export_field_names(11) = "Fioo_frazil" ! get mct_avect_index value for each export field name call get_indices_from_names( & diff --git a/components/omega/src/drivers/coupled/omega_cxx2f_interface.cpp b/components/omega/src/drivers/coupled/omega_cxx2f_interface.cpp index 2347521336f6..b527bbc69cf1 100644 --- a/components/omega/src/drivers/coupled/omega_cxx2f_interface.cpp +++ b/components/omega/src/drivers/coupled/omega_cxx2f_interface.cpp @@ -13,6 +13,8 @@ #include "TimeStepper.h" #include +bool iOwnKokkos = false; + // helper C++ functions namespace { @@ -67,7 +69,10 @@ void omega_ocn_init1( MPI_Comm Comm = MPI_Comm_f2c(FComm); // initialize Kokkos - Kokkos::initialize(); + if (!Kokkos::is_initialized() && !Kokkos::is_finalized()) { + Kokkos::initialize(); + iOwnKokkos = true; + } // initialize Pacer timing in coupled mode Pacer::initialize(Comm, Pacer::PACER_INTEGRATED); @@ -153,7 +158,9 @@ int omega_ocn_finalize() { // no Pacer::print or Pacer::finalize in coupled mode; cpl will handle it // finalize Kokkos - Kokkos::finalize(); + if (!Kokkos::is_finalized() && iOwnKokkos) { + Kokkos::finalize(); + } return ErrFinalize; } diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 16ea4da9f5a2..ba36a2bb61ac 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -327,22 +327,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } -Real Eos::calcPtFromCt(const Real &Sa, const Real &Ct) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct); - } - - return Ct; -} - -Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const { - if (EosChoice == EosType::Teos10Eos) { - return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); - } - - return Pt; -} - /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 5e4fd89cae21..2494d6a2d5d8 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -358,8 +358,8 @@ class Teos10Eos { /// (polynomial error in [-5e-4, 6e-4] K, from GSW package). /// P is relative pressure (gauge pressure in Pa, i.e., absolute pressure /// minus the standard atmosphere). - KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { + static KOKKOS_FUNCTION Real calcCtFreezingTeos10( + const Real Sa, const Real P, const Real SaturationFract) { constexpr Real Sso = 35.16504; constexpr Real C0 = 0.017947064327968736; constexpr Real C1 = -6.076099099929818; @@ -756,10 +756,47 @@ class Eos { const Array2DReal &SpecVol); /// Convert Conservative Temperature to potential temperature - Real calcPtFromCt(const Real &Sa, const Real &Ct) const; + /// For TEOS-10, uses the TEOS-10 polynomial + /// For other EOS choices, conservative temperature is equal to potential + /// temperature + KOKKOS_FUNCTION Real calcPtFromCt(const Real &Sa, const Real &Ct) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct); + } + return Ct; + } + + /// Convert potential temperature to Conservative Temperature. + /// For TEOS-10, uses the TEOS-10 polynomial + /// For other EOS choices, potential temperature equals conservative + /// temperature + KOKKOS_FUNCTION Real calcCtFromPt(const Real &Sa, const Real &Pt) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); + } + return Pt; + } - /// Convert potential temperature to Conservative Temperature - Real calcCtFromPt(const Real &Sa, const Real &Pt) const; + /// Calculate freezing temperature of seawater. + /// For TEOS-10, uses the Roquet et al. 75-term polynomial. + /// For LinearEos, uses a simple linear salinity-dependent approximation. + /// For ConstantEos, returns a constant approximate ocean freezing point. + static KOKKOS_FUNCTION Real calcCtFreezing(EosType Choice, const Real Sa, + const Real P, + const Real SaturationFract) { + if (Choice == EosType::Teos10Eos) { + return Teos10Eos::calcCtFreezingTeos10(Sa, P, SaturationFract); + } + if (Choice == EosType::LinearEos) { + // Linear salinity-dependent freezing point; coefficient -0.054 + // degC/PSU with absolute-to-practical salinity conversion (g/kg -> + // PSU). + constexpr Real Coeff = -0.054_Real; + return Coeff * Sa / Psu2Gpkg; + } + // ConstantEos: constant approximate ocean freezing point (degC) + return -1.9_Real; + } /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/Forcing.cpp b/components/omega/src/ocn/Forcing.cpp index 42643dba5385..46af26094113 100644 --- a/components/omega/src/ocn/Forcing.cpp +++ b/components/omega/src/ocn/Forcing.cpp @@ -29,18 +29,30 @@ static std::string stripDefault(const std::string &Name) { // mesh/halo. Forcing::Forcing(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo) : Name(stripDefault(Name)), SfcStressForcing(stripDefault(Name), Mesh), - Mesh(Mesh), MeshHalo(MeshHalo) {} + TracerForcing(stripDefault(Name), Mesh), Mesh(Mesh), MeshHalo(MeshHalo) {} // Destructor. Unregisters fields from IO streams. Forcing::~Forcing() { unregisterFields(); } // Register surface stress fields with IO streams for a given mesh. void Forcing::registerFields(const std::string &MeshName) const { - SfcStressForcing.registerFields(MeshName); + if (SfcStressFieldsEnabled) { + SfcStressForcing.registerFields(MeshName); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.registerFields(MeshName); + } } // Unregister surface stress fields from IO streams. -void Forcing::unregisterFields() const { SfcStressForcing.unregisterFields(); } +void Forcing::unregisterFields() const { + if (SfcStressFieldsEnabled) { + SfcStressForcing.unregisterFields(); + } + if (TracerForcingFieldsEnabled) { + TracerForcing.unregisterFields(); + } +} // Create and register a non-default forcing instance. Forcing *Forcing::create(const std::string &Name, const HorzMesh *Mesh, @@ -76,10 +88,9 @@ void Forcing::init() { ABORT_ERROR("Forcing: failed to initialize default forcing state"); } - DefaultForcing->registerFields(DefMesh->MeshName); - Config *OmegaConfig = Config::getOmegaConfig(); DefaultForcing->readConfigOptions(OmegaConfig); + DefaultForcing->registerFields(DefMesh->MeshName); // for now, forcing fields are read at start-up only. // to be extended to include switch from standalone to coupled. // to be moved to a Forcing->prepareForStep(SimTime) method later. @@ -136,10 +147,62 @@ void Forcing::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("Forcing: Unknown InterpType requested"); } + + Config TendConfig("Tendencies"); + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, "Forcing: Tendencies group not found in Config"); + + Err += + TendConfig.get("SfcStressForcingTendencyEnable", SfcStressFieldsEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcStressForcingTendencyEnable not found " + "in Tendencies config"); + + bool SfcThicknessForcingEnabled = false; + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + SfcThicknessForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcThicknessForcingTendencyEnable not " + "found in Tendencies config"); + + bool SfcTracerForcingEnabled = false; + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + SfcTracerForcingEnabled); + CHECK_ERROR_ABORT(Err, "Forcing: SfcTracerForcingTendencyEnable not found " + "in Tendencies config"); + + TracerForcingFieldsEnabled = + SfcThicknessForcingEnabled || SfcTracerForcingEnabled; } // Compute all forcing variables (dispatches to specific computations). -void Forcing::computeAll() const { computeSfcStressForcingOnEdge(); } +void Forcing::computeAll() const { + exchangeHalo(); + computeSfcStressForcingOnEdge(); +} + +// Reset forcing arrays so omitted optional fields remain zero after read. +void Forcing::resetArrays() { + if (SfcStressFieldsEnabled) { + deepCopy(SfcStressForcing.NormalStressEdge, 0.0_Real); + deepCopy(SfcStressForcing.ZonalStressCell, 0.0_Real); + deepCopy(SfcStressForcing.MeridStressCell, 0.0_Real); + } + + if (TracerForcingFieldsEnabled) { + deepCopy(TracerForcing.SnowFluxCell, 0.0_Real); + deepCopy(TracerForcing.RainFluxCell, 0.0_Real); + deepCopy(TracerForcing.EvaporationFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(TracerForcing.IceRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.RiverRunoffFluxCell, 0.0_Real); + deepCopy(TracerForcing.LatentHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SensibleHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(TracerForcing.LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(TracerForcing.SeaIceSaltFluxCell, 0.0_Real); + } +} // Compute edge-normal stress from cell-center zonal and meridional components. void Forcing::computeSfcStressForcingOnEdge() const { @@ -153,14 +216,17 @@ void Forcing::computeSfcStressForcingOnEdge() const { Pacer::stop("Forcing:edge1", 2); } -// Exchange halo for surface stress cell fields. +// Exchange halo for surface stress cell fields. Only needed for variables that +// need information beyond cell-centered values. I4 Forcing::exchangeHalo() const { I4 Err = 0; - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, - OnCell); - Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, - OnCell); + if (SfcStressFieldsEnabled) { + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.ZonalStressCell, + OnCell); + Err += MeshHalo->exchangeFullArrayHalo(SfcStressForcing.MeridStressCell, + OnCell); + } return Err; } @@ -172,13 +238,20 @@ void Forcing::readStreamIntoArrays() { std::string StreamName = "Forcing"; + resetArrays(); + + // Nothing to read if neither stress nor tracer forcing tendencies are + // enabled. + if (!SfcStressFieldsEnabled && !TracerForcingFieldsEnabled) { + return; + } + // Attempt to read stream; if unavailable, log and fall back to zero forcing. Err = IOStream::read(StreamName); if (Err.isFail()) { LOG_INFO("Forcing: Error while reading {} stream, using zero forcing", StreamName); - deepCopy(SfcStressForcing.ZonalStressCell, 0._Real); - deepCopy(SfcStressForcing.MeridStressCell, 0._Real); + resetArrays(); } I4 HaloErr = exchangeHalo(); diff --git a/components/omega/src/ocn/Forcing.h b/components/omega/src/ocn/Forcing.h index 5fdae7e550b8..b061bdcf25d7 100644 --- a/components/omega/src/ocn/Forcing.h +++ b/components/omega/src/ocn/Forcing.h @@ -17,6 +17,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "forcingVars/SfcStressForcingVars.h" +#include "forcingVars/TracerForcingVars.h" #include #include @@ -32,6 +33,7 @@ class Forcing { std::string Name; ///< Name identifier for this forcing instance SfcStressForcingVars SfcStressForcing; ///< Surface stress forcing variables + TracerForcingVars TracerForcing; ///< Tracer forcing vars (thickness and T,S) ~Forcing(); @@ -69,6 +71,9 @@ class Forcing { /// Read forcing fields from input stream at startup void readStreamIntoArrays(); + /// Reset all forcing arrays to zero before reading optional fields + void resetArrays(); + /// Compute all forcing variables void computeAll() const; @@ -86,6 +91,8 @@ class Forcing { const HorzMesh *Mesh; Halo *MeshHalo; + bool SfcStressFieldsEnabled = false; + bool TracerForcingFieldsEnabled = false; static Forcing *DefaultForcing; static std::map> AllForcing; diff --git a/components/omega/src/ocn/GlobalConstants.h b/components/omega/src/ocn/GlobalConstants.h index d31bae11f489..0fe2adb1f54e 100644 --- a/components/omega/src/ocn/GlobalConstants.h +++ b/components/omega/src/ocn/GlobalConstants.h @@ -115,11 +115,12 @@ constexpr Real Pa2Db = 1.0e-4; // Pascal to Decibar constexpr Real Cm2M = 1.0e-2; // Centimeters to meters constexpr Real M2Cm = 1.0e2; // Meters to centimeters constexpr Real HFluxFac = - 1.0 / (RhoSw * CpSw); // Heat flux (W/m^2) to temp flux (C*m/s) + 1.0 / (RhoSw * Cp0Sw); // Heat flux (W/m^2) to Conserv Temp flux (C*m/s) constexpr Real FwFluxFac = 1.e-6; // Fw flux (kg/m^2/s) to salt((msu/psu)*m/s) constexpr Real SaltFac = - -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s) -constexpr Real SFluxFac = 1.0; // Salt flux (kg/m^2/s) to salt flux (msu*m/s) + -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s) +constexpr Real SFluxFac = + 1.e3 / RhoSw; // Salt flux (kg/m^2/s) to salinity flux (m*(g/kg)/s) } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/SfcCoupling.cpp b/components/omega/src/ocn/SfcCoupling.cpp index ecb28d465553..2d4d35c58dc4 100644 --- a/components/omega/src/ocn/SfcCoupling.cpp +++ b/components/omega/src/ocn/SfcCoupling.cpp @@ -1,6 +1,7 @@ #include "SfcCoupling.h" #include "Eos.h" #include "GlobalConstants.h" +#include "HorzOperators.h" #include "Logging.h" #include "OceanState.h" #include "OmegaKokkos.h" @@ -66,8 +67,9 @@ SfcCoupling::SfcCoupling(const std::string &Name_, const HorzMesh *Mesh, ImportIdxMap(ImportIdxMap), ExportIdxMap(ExportIdxMap), CplToOcn(Name_, Mesh), OcnToCpl(Name_, Mesh), Layout(Layout) { - // Retrieve mesh cell count + // Retrieve mesh cell/edge count NCellsOwned = Mesh->NCellsOwned; + NEdgesAll = Mesh->NEdgesAll; NAccumSteps = 0; @@ -187,14 +189,40 @@ void SfcCoupling::importFromCoupler() { "method must be called before importing data from the coupler."); } - // Get import field indices for surface stress components - int TauxIdx = ImportIdxMap.at("Foxx_taux"); - int TauyIdx = ImportIdxMap.at("Foxx_tauy"); + // Get import field indices + int TauxIdx = ImportIdxMap.at("Foxx_taux"); + int TauyIdx = ImportIdxMap.at("Foxx_tauy"); + int SwnetIdx = ImportIdxMap.at("Foxx_swnet"); + int SenIdx = ImportIdxMap.at("Foxx_sen"); + int LatIdx = ImportIdxMap.at("Foxx_lat"); + int LwupIdx = ImportIdxMap.at("Foxx_lwup"); + int LwdnIdx = ImportIdxMap.at("Faxa_lwdn"); + int SaltIdx = ImportIdxMap.at("Fioi_salt"); + int MelthIdx = ImportIdxMap.at("Fioi_melth"); + int MeltwIdx = ImportIdxMap.at("Fioi_meltw"); + int SnowIdx = ImportIdxMap.at("Faxa_snow"); + int RainIdx = ImportIdxMap.at("Faxa_rain"); + int EvapIdx = ImportIdxMap.at("Foxx_evap"); + int RoflIdx = ImportIdxMap.at("Foxx_rofl"); + int RofiIdx = ImportIdxMap.at("Foxx_rofi"); // Copy Kokkos view handles - auto CplToOcnView_ = CplToOcnView; - auto SfcStressZonal_ = CplToOcn.SfcStressZonal; - auto SfcStressMerid_ = CplToOcn.SfcStressMerid; + auto CplToOcnView_ = CplToOcnView; + auto SfcStressZonal_ = CplToOcn.SfcStressZonal; + auto SfcStressMerid_ = CplToOcn.SfcStressMerid; + auto SnowFlux_ = CplToOcn.SnowFlux; + auto RainFlux_ = CplToOcn.RainFlux; + auto EvaporationFlux_ = CplToOcn.EvaporationFlux; + auto SeaIceFreshWaterFlux_ = CplToOcn.SeaIceFreshWaterFlux; + auto IceRunoffFlux_ = CplToOcn.IceRunoffFlux; + auto RiverRunoffFlux_ = CplToOcn.RiverRunoffFlux; + auto LatentHeatFlux_ = CplToOcn.LatentHeatFlux; + auto SensibleHeatFlux_ = CplToOcn.SensibleHeatFlux; + auto LongWaveHeatFluxUp_ = CplToOcn.LongWaveHeatFluxUp; + auto LongWaveHeatFluxDown_ = CplToOcn.LongWaveHeatFluxDown; + auto SeaIceHeatFlux_ = CplToOcn.SeaIceHeatFlux; + auto ShortWaveHeatFlux_ = CplToOcn.ShortWaveHeatFlux; + auto SeaIceSaltFlux_ = CplToOcn.SeaIceSaltFlux; /// TODO: Shouldn't be making direct calls to Kokkos here. /// How often is threading used? Becuase this will be a serial loop @@ -202,8 +230,21 @@ void SfcCoupling::importFromCoupler() { auto Policy = Kokkos::RangePolicy>( 0, NCellsOwned); Kokkos::parallel_for("importFromCoupler", Policy, [=](int Idx) { - SfcStressZonal_(Idx) = CplToOcnView_(TauxIdx, Idx); - SfcStressMerid_(Idx) = CplToOcnView_(TauyIdx, Idx); + SfcStressZonal_(Idx) = CplToOcnView_(TauxIdx, Idx); + SfcStressMerid_(Idx) = CplToOcnView_(TauyIdx, Idx); + SnowFlux_(Idx) = CplToOcnView_(SnowIdx, Idx); + RainFlux_(Idx) = CplToOcnView_(RainIdx, Idx); + EvaporationFlux_(Idx) = CplToOcnView_(EvapIdx, Idx); + SeaIceFreshWaterFlux_(Idx) = CplToOcnView_(MeltwIdx, Idx); + IceRunoffFlux_(Idx) = CplToOcnView_(RofiIdx, Idx); + RiverRunoffFlux_(Idx) = CplToOcnView_(RoflIdx, Idx); + LatentHeatFlux_(Idx) = CplToOcnView_(LatIdx, Idx); + SensibleHeatFlux_(Idx) = CplToOcnView_(SenIdx, Idx); + LongWaveHeatFluxUp_(Idx) = CplToOcnView_(LwupIdx, Idx); + LongWaveHeatFluxDown_(Idx) = CplToOcnView_(LwdnIdx, Idx); + SeaIceHeatFlux_(Idx) = CplToOcnView_(MelthIdx, Idx); + ShortWaveHeatFlux_(Idx) = CplToOcnView_(SwnetIdx, Idx); + SeaIceSaltFlux_(Idx) = CplToOcnView_(SaltIdx, Idx); }); } @@ -223,6 +264,8 @@ void SfcCoupling::exportToCoupler() { int VelUIdx = ExportIdxMap.at("So_u"); int VelVIdx = ExportIdxMap.at("So_v"); int SshIdx = ExportIdxMap.at("So_ssh"); + int DhdxIdx = ExportIdxMap.at("So_dhdx"); + int DhdyIdx = ExportIdxMap.at("So_dhdy"); // Copy Kokkos view handles auto OcnToCplView_ = OcnToCplView; @@ -230,6 +273,8 @@ void SfcCoupling::exportToCoupler() { auto AvgSfcSalinity_ = OcnToCpl.AvgSfcSalinityH; auto AvgSfcVelocityZonal_ = OcnToCpl.AvgSfcVelocityZonalH; auto AvgSfcVelocityMerid_ = OcnToCpl.AvgSfcVelocityMeridH; + auto AvgSfcSshGradZonal_ = OcnToCpl.AvgSfcSshGradZonalH; + auto AvgSfcSshGradMerid_ = OcnToCpl.AvgSfcSshGradMeridH; auto InstSshCellH_ = OcnToCpl.InstSshCellH; // Initalize all o2x fields to 0.0 for next coupling interval @@ -243,6 +288,8 @@ void SfcCoupling::exportToCoupler() { OcnToCplView_(SalinIdx, Idx) = AvgSfcSalinity_(Idx); OcnToCplView_(VelUIdx, Idx) = AvgSfcVelocityZonal_(Idx); OcnToCplView_(VelVIdx, Idx) = AvgSfcVelocityMerid_(Idx); + OcnToCplView_(DhdxIdx, Idx) = AvgSfcSshGradZonal_(Idx); + OcnToCplView_(DhdyIdx, Idx) = AvgSfcSshGradMerid_(Idx); OcnToCplView_(SshIdx, Idx) = InstSshCellH_(Idx); }); @@ -258,63 +305,109 @@ void SfcCoupling::applyImportFields(Forcing *Forcing) { CplToOcn.SfcStressZonal); deepCopy(ownedSubView(Forcing->SfcStressForcing.MeridStressCell), CplToOcn.SfcStressMerid); + + deepCopy(ownedSubView(Forcing->TracerForcing.SnowFluxCell), + CplToOcn.SnowFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.RainFluxCell), + CplToOcn.RainFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.EvaporationFluxCell), + CplToOcn.EvaporationFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.SeaIceFreshWaterFluxCell), + CplToOcn.SeaIceFreshWaterFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.IceRunoffFluxCell), + CplToOcn.IceRunoffFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.RiverRunoffFluxCell), + CplToOcn.RiverRunoffFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.LatentHeatFluxCell), + CplToOcn.LatentHeatFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.SensibleHeatFluxCell), + CplToOcn.SensibleHeatFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.LongWaveHeatFluxUpCell), + CplToOcn.LongWaveHeatFluxUp); + deepCopy(ownedSubView(Forcing->TracerForcing.LongWaveHeatFluxDownCell), + CplToOcn.LongWaveHeatFluxDown); + deepCopy(ownedSubView(Forcing->TracerForcing.SeaIceHeatFluxCell), + CplToOcn.SeaIceHeatFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.ShortWaveHeatFluxCell), + CplToOcn.ShortWaveHeatFlux); + deepCopy(ownedSubView(Forcing->TracerForcing.SeaIceSaltFluxCell), + CplToOcn.SeaIceSaltFlux); }; void SfcCoupling::updateExportFields(const OceanState *State, const Array3DReal &TracerArray) { - OcnToCpl.updateFields(State, TracerArray, NAccumSteps, NCellsOwned); + OcnToCpl.updateFields(State, TracerArray, NAccumSteps, NCellsOwned, + NEdgesAll); NAccumSteps++; } CplToOcnFields::CplToOcnFields(const std::string &Suffix, const HorzMesh *Mesh) : SfcStressZonal("SfcStressZonal" + Suffix, Mesh->NCellsOwned), - SfcStressMerid("SfcStressMeridional" + Suffix, Mesh->NCellsOwned) {} + SfcStressMerid("SfcStressMeridional" + Suffix, Mesh->NCellsOwned), + SnowFlux("SnowFlux" + Suffix, Mesh->NCellsOwned), + RainFlux("RainFlux" + Suffix, Mesh->NCellsOwned), + EvaporationFlux("EvaporationFlux" + Suffix, Mesh->NCellsOwned), + SeaIceFreshWaterFlux("SeaIceFreshWaterFlux" + Suffix, Mesh->NCellsOwned), + IceRunoffFlux("IceRunoffFlux" + Suffix, Mesh->NCellsOwned), + RiverRunoffFlux("RiverRunoffFlux" + Suffix, Mesh->NCellsOwned), + LatentHeatFlux("LatentHeatFlux" + Suffix, Mesh->NCellsOwned), + SensibleHeatFlux("SensibleHeatFlux" + Suffix, Mesh->NCellsOwned), + LongWaveHeatFluxUp("LongWaveHeatFluxUp" + Suffix, Mesh->NCellsOwned), + LongWaveHeatFluxDown("LongWaveHeatFluxDown" + Suffix, Mesh->NCellsOwned), + SeaIceHeatFlux("SeaIceHeatFlux" + Suffix, Mesh->NCellsOwned), + ShortWaveHeatFlux("ShortWaveHeatFlux" + Suffix, Mesh->NCellsOwned), + SeaIceSaltFlux("SeaIceSaltFlux" + Suffix, Mesh->NCellsOwned) {} OcnToCplFields::OcnToCplFields(const std::string &Suffix, const HorzMesh *Mesh) : AvgSfcTemperature("AvgSfcTemperature" + Suffix, Mesh->NCellsOwned), AvgSfcSalinity("AvgSfcSalinity" + Suffix, Mesh->NCellsOwned), - AvgSfcVelocityZonal("AvgSfcVelocityZonal" + Suffix, Mesh->NCellsOwned), - AvgSfcVelocityMerid("AvgSfcVelocityMeridional" + Suffix, + AvgSfcNormalVelocity("AvgSfcNormalVelocity" + Suffix, Mesh->NEdgesSize), + AvgSfcVelocityZonalH("AvgSfcVelocityZonal" + Suffix, Mesh->NCellsOwned), + AvgSfcVelocityMeridH("AvgSfcVelocityMeridional" + Suffix, + Mesh->NCellsOwned), + AvgSfcSshGrad("AvgSfcSshGrad" + Suffix, Mesh->NEdgesSize), + AvgSfcSshGradZonalH("AvgSfcSshGradZonal" + Suffix, Mesh->NCellsOwned), + AvgSfcSshGradMeridH("AvgSfcSshGradMeridional" + Suffix, Mesh->NCellsOwned), InstSshCellH("InstSshCellH" + Suffix, Mesh->NCellsOwned), - InSituTempScratch("InSituTempScratch" + Suffix, Mesh->NCellsOwned) { + InSituTempScratch("InSituTempScratch" + Suffix, Mesh->NCellsOwned), + PracSalinityScratch("PracSalinityScratch" + Suffix, Mesh->NCellsOwned), + ReconZonalScratch("ReconZonalScratch" + Suffix, Mesh->NCellsOwned), + ReconMeridScratch("ReconMeridScratch" + Suffix, Mesh->NCellsOwned) { // Kokkok views created with a label are zero-initialized by default. // We reset the fields here anyway to be explicit about the fact that the // OcnToCpl fields need to begin a coupling interval with all zeros. resetFields(); - AvgSfcTemperatureH = createHostMirrorCopy(AvgSfcTemperature); - AvgSfcSalinityH = createHostMirrorCopy(AvgSfcSalinity); - AvgSfcVelocityZonalH = createHostMirrorCopy(AvgSfcVelocityZonal); - AvgSfcVelocityMeridH = createHostMirrorCopy(AvgSfcVelocityMerid); + AvgSfcTemperatureH = createHostMirrorCopy(AvgSfcTemperature); + AvgSfcSalinityH = createHostMirrorCopy(AvgSfcSalinity); } void OcnToCplFields::updateFields(const OceanState *State, const Array3DReal &TracerArray, - const I4 NAccumSteps, const I4 NCellsOwned) { + const I4 NAccumSteps, const I4 NCellsOwned, + const I4 NEdgesAll) { I4 TemperatureIdx, SalinityIdx; Tracers::getIndex(TemperatureIdx, "Temperature"); Tracers::getIndex(SalinityIdx, "Salinity"); + // get normal velocity at current time level + auto NormalVel = State->getNormalVelocity(0); auto Temperature = Kokkos::subview(TracerArray, TemperatureIdx, Kokkos::ALL, Kokkos::ALL); auto Salinity = Kokkos::subview(TracerArray, SalinityIdx, Kokkos::ALL, Kokkos::ALL); + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); VertCoord *DefVertCoord = VertCoord::getDefault(); OMEGA_SCOPE(LocMinLayerCell, DefVertCoord->MinLayerCell); OMEGA_SCOPE(LocAvgSfcSalinity, AvgSfcSalinity); OMEGA_SCOPE(LocAvgSfcTemp, AvgSfcTemperature); - OMEGA_SCOPE(LocAvgSfcVelZonal, AvgSfcVelocityZonal); - OMEGA_SCOPE(LocAvgSfcVelMerid, AvgSfcVelocityMerid); - - // TODO: Implement vector reconsturction for velocity field. - constexpr Real ConstSfcVelocity = 1e-4; parallelFor( {NCellsOwned}, KOKKOS_LAMBDA(int ICell) { @@ -326,12 +419,35 @@ void OcnToCplFields::updateFields(const OceanState *State, LocAvgSfcSalinity(ICell) = updateAverage( LocAvgSfcSalinity(ICell), Salinity(ICell, KSfc), NAccumSteps); + }); - LocAvgSfcVelZonal(ICell) = updateAverage( - LocAvgSfcVelZonal(ICell), ConstSfcVelocity, NAccumSteps); + OMEGA_SCOPE(LocCellsOnEdge, DefHorzMesh->CellsOnEdge); + OMEGA_SCOPE(LocDcEdge, DefHorzMesh->DcEdge); + OMEGA_SCOPE(LocSshCell, DefVertCoord->SshCell); + OMEGA_SCOPE(LocMinLayerEdgeBot, DefVertCoord->MinLayerEdgeBot); + OMEGA_SCOPE(LocMaxLayerEdgeTop, DefVertCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocAvgSfcSshGrad, AvgSfcSshGrad); + OMEGA_SCOPE(LocAvgSfcNormalVel, AvgSfcNormalVelocity); - LocAvgSfcVelMerid(ICell) = updateAverage( - LocAvgSfcVelMerid(ICell), ConstSfcVelocity, NAccumSteps); + parallelFor( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + + // exclude outer halo edges whose layer range is invalid + if (KMin <= KMax) { + LocAvgSfcNormalVel(IEdge) = + updateAverage(LocAvgSfcNormalVel(IEdge), + NormalVel(IEdge, KMin), NAccumSteps); + + const int ICell0 = LocCellsOnEdge(IEdge, 0); + const int ICell1 = LocCellsOnEdge(IEdge, 1); + const Real SshGrad = + (LocSshCell(ICell1) - LocSshCell(ICell0)) / LocDcEdge(IEdge); + + LocAvgSfcSshGrad(IEdge) = + updateAverage(LocAvgSfcSshGrad(IEdge), SshGrad, NAccumSteps); + } }); } @@ -355,21 +471,55 @@ void OcnToCplFields::copyToHost() { OMEGA_SCOPE(LocAvgSfcTemp, AvgSfcTemperature); OMEGA_SCOPE(LocAvgSfcSalinity, AvgSfcSalinity); OMEGA_SCOPE(LocInSituTemp, InSituTempScratch); - + OMEGA_SCOPE(LocPracSalinity, PracSalinityScratch); + + // TEOS-10 conversion is applied once per coupling interval to the averaged + // conservative temperature and absolute salinity. Therefore this computes + // PtFromCt(mean(Sa), mean(Ct)), not mean(PtFromCt(Sa, Ct)); these are not + // generally equivalent because the conversion is nonlinear. A true time + // average of the converted quantity would require converting each timestep + // before accumulating it. The same consideration applies to any future + // nonlinear absolute-to-practical salinity conversion. parallelFor( {(int)AvgSfcTemperature.extent(0)}, KOKKOS_LAMBDA(int Cell) { - const Real Ct = LocAvgSfcTemp(Cell); - const Real Sa = LocAvgSfcSalinity(Cell); - const Real Pt = LocEosChoice == EosType::Teos10Eos - ? LocTeos10.calcPtFromCt(Sa, Ct) - : Ct; - LocInSituTemp(Cell) = Pt + TkFrz; // C to K temperature conversion + const Real Ct = LocAvgSfcTemp(Cell); + const Real Sa = LocAvgSfcSalinity(Cell); + const Real Pt = LocEosChoice == EosType::Teos10Eos + ? LocTeos10.calcPtFromCt(Sa, Ct) + : Ct; + LocInSituTemp(Cell) = Pt + TkFrz; // C to K temperature conversion + LocPracSalinity(Cell) = Sa / Psu2Gpkg; // abs to prac sal. conversion }); deepCopy(AvgSfcTemperatureH, InSituTempScratch); - deepCopy(AvgSfcSalinityH, AvgSfcSalinity); - deepCopy(AvgSfcVelocityZonalH, AvgSfcVelocityZonal); - deepCopy(AvgSfcVelocityMeridH, AvgSfcVelocityMerid); + deepCopy(AvgSfcSalinityH, PracSalinityScratch); + + // Retrieve the default horizontal mesh + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); + + OMEGA_SCOPE(LocNormalVelocity, AvgSfcNormalVelocity); + OMEGA_SCOPE(LocSshGradEdge, AvgSfcSshGrad); + OMEGA_SCOPE(LocReconZonal, ReconZonalScratch); + OMEGA_SCOPE(LocReconMerid, ReconMeridScratch); + + VectorReconOnCell ReconCell(DefHorzMesh); + + parallelFor( + {DefHorzMesh->NCellsOwned}, KOKKOS_LAMBDA(int ICell) { + ReconCell(LocReconZonal, LocReconMerid, ICell, LocNormalVelocity); + }); + + deepCopy(AvgSfcVelocityZonalH, ReconZonalScratch); + deepCopy(AvgSfcVelocityMeridH, ReconMeridScratch); + + // Reuse the same scratch arrays; VectorReconOnCell overwrites every cell + parallelFor( + {DefHorzMesh->NCellsOwned}, KOKKOS_LAMBDA(int ICell) { + ReconCell(LocReconZonal, LocReconMerid, ICell, LocSshGradEdge); + }); + + deepCopy(AvgSfcSshGradZonalH, ReconZonalScratch); + deepCopy(AvgSfcSshGradMeridH, ReconMeridScratch); // SSH is an instantaneous field, so we don't bother with a device mirror of // our own. Instead, copy from the VertCoord, which owns SSH, host array. @@ -385,7 +535,7 @@ void OcnToCplFields::copyToHost() { void OcnToCplFields::resetFields() { deepCopy(AvgSfcTemperature, 0.0_Real); deepCopy(AvgSfcSalinity, 0.0_Real); - deepCopy(AvgSfcVelocityZonal, 0.0_Real); - deepCopy(AvgSfcVelocityMerid, 0.0_Real); + deepCopy(AvgSfcNormalVelocity, 0.0_Real); + deepCopy(AvgSfcSshGrad, 0.0_Real); } } // namespace OMEGA diff --git a/components/omega/src/ocn/SfcCoupling.h b/components/omega/src/ocn/SfcCoupling.h index 4158e4861c81..0913b023cffb 100644 --- a/components/omega/src/ocn/SfcCoupling.h +++ b/components/omega/src/ocn/SfcCoupling.h @@ -50,6 +50,22 @@ class CplToOcnFields { HostArray1DReal SfcStressZonal; ///< Foxx_taux [N m^-2] HostArray1DReal SfcStressMerid; ///< Foxx_tauy [N m^-2] + HostArray1DReal SnowFlux; ///< Faxa_snow [kg m^-2 s^-1] + HostArray1DReal RainFlux; ///< Faxa_rain [kg m^-2 s^-1] + HostArray1DReal EvaporationFlux; ///< Foxx_evap [kg m^-2 s^-1] + HostArray1DReal SeaIceFreshWaterFlux; ///< Fioi_meltw [kg m^-2 s^-1] + HostArray1DReal IceRunoffFlux; ///< Foxx_rofi [kg m^-2 s^-1] + HostArray1DReal RiverRunoffFlux; ///< Foxx_rofl [kg m^-2 s^-1] + + HostArray1DReal LatentHeatFlux; ///< Foxx_lat [W m^-2] + HostArray1DReal SensibleHeatFlux; ///< Foxx_sen [W m^-2] + HostArray1DReal LongWaveHeatFluxUp; ///< Foxx_lwup [W m^-2] + HostArray1DReal LongWaveHeatFluxDown; ///< Faxa_lwdn [W m^-2] + HostArray1DReal SeaIceHeatFlux; ///< Fioi_melth [W m^-2] + HostArray1DReal ShortWaveHeatFlux; ///< Foxx_swnet [W m^-2] + + HostArray1DReal SeaIceSaltFlux; ///< Fioi_salt [kg m^-2 s^-1] + CplToOcnFields(const std::string &Suffix, const HorzMesh *Mesh); }; @@ -59,8 +75,7 @@ class OcnToCplFields { ///< So_t [K], in-situ approx (potential temp at P=0) HostArray1DReal AvgSfcTemperatureH; - /// TODO: Export practical salinity (unitless) to coupler - ///< So_s [g kg^-1], absolute salinity + ///< So_s [psu], paractical salinity HostArray1DReal AvgSfcSalinityH; ///< So_u [m s^-1] @@ -69,13 +84,18 @@ class OcnToCplFields { ///< So_v [m s^-1] HostArray1DReal AvgSfcVelocityMeridH; + ///< So_dhdx [m m-1], zonal sea surface slope + HostArray1DReal AvgSfcSshGradZonalH; + ///< So_dhdy [m m-1], meridional sea surface slope + HostArray1DReal AvgSfcSshGradMeridH; + ///< So_ssh [m] /// instantaneous field, so no device mirror is needed HostArray1DReal InstSshCellH; // Accumulate one ocean timestep's contribution to the running averages void updateFields(const OceanState *State, const Array3DReal &TracerArray, - I4 NAccumSteps, I4 NCellsOwned); + I4 NAccumSteps, I4 NCellsOwned, I4 NEdgesAll); // Copy device arrays into their host mirrors and do unit conversion. void copyToHost(); @@ -92,11 +112,17 @@ class OcnToCplFields { // the rest of the code. Array1DReal AvgSfcTemperature; // [C], conservative temperature Array1DReal AvgSfcSalinity; // [g kg^-1], absolute salinity - Array1DReal AvgSfcVelocityZonal; - Array1DReal AvgSfcVelocityMerid; - // Scratch buffer for the in-situ Kelvin conversion in copyToHost() + Array1DReal AvgSfcNormalVelocity; // [m s^-1], velocity normal to edge + Array1DReal AvgSfcSshGrad; // [m m^-1], ssh gradient normal to edge + + // Scratch buffers for the in-situ and Kelvin temperature conversion and + // paractical salinity conversio done in copyToHost() Array1DReal InSituTempScratch; // [K], in-situ approx (potential temp at P=0) + Array1DReal PracSalinityScratch; // [Psu], Parctical salinity + // Scratch arrays for edge normal vector field reconstructed to cell centers + Array1DReal ReconZonalScratch; + Array1DReal ReconMeridScratch; }; /// A class for interfacing with the coupler @@ -142,6 +168,7 @@ class SfcCoupling { std::string Name; I4 NCellsOwned; ///< Number of cells owned by this task + I4 NEdgesAll; ///< Total number (owned+halo) of local edges // The values below will be larger than InportIdx.size() and // ExportIdxMap.size() because omega does not ingest all cpl fields (e.g. diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index d8273a844de0..419c105292a6 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -283,6 +283,18 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config } } + Err += TendConfig.get("SfcThicknessForcingTendencyEnable", + this->SfcThicknessForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcThicknessForcingTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("SfcTracerForcingTendencyEnable", + this->SfcTracerForcing.Enabled); + CHECK_ERROR_ABORT( + Err, + "Tendencies: SfcTracerForcingTendencyEnable not found in TendConfig"); + if (this->TracerDiffusion.Enabled) { Err += TendConfig.get("EddyDiff2", this->TracerDiffusion.EddyDiff2); CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff2 not found in TendConfig"); @@ -461,6 +473,9 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord), VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), ExplicitBottomDrag(Mesh, VCoord), + SfcThicknessForcing(Mesh, VCoord), + SfcTracerForcing(Mesh, VCoord, Tracers::IndxTemp, Tracers::IndxSalt, + EqState), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), @@ -510,6 +525,7 @@ void Tendencies::computePseudoThicknessTendenciesOnly( OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); OMEGA_SCOPE(LocThicknessFluxDiv, PseudoThicknessFluxDiv); + OMEGA_SCOPE(LocSfcThicknessForcing, SfcThicknessForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -553,6 +569,32 @@ void Tendencies::computePseudoThicknessTendenciesOnly( VAdv->computePseudoThicknessVAdvTend(PseudoThicknessTend); Pacer::stop("Tend:computePseudoThicknessVAdvTend", 2); + if (LocSfcThicknessForcing.Enabled) { + Pacer::start("Tend:sfcThicknessForcing", 2); + const auto *ForcingState = Forcing::getDefault(); + + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &EvaporationFlux = + ForcingState->TracerForcing.EvaporationFluxCell; + const auto &SeaIceFreshWaterFlux = + ForcingState->TracerForcing.SeaIceFreshWaterFluxCell; + const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &RiverRunoffFlux = + ForcingState->TracerForcing.RiverRunoffFluxCell; + const auto &SeaIceSaltFlux = + ForcingState->TracerForcing.SeaIceSaltFluxCell; + + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + LocSfcThicknessForcing(LocPseudoThicknessTend, ICell, SnowFlux, + RainFlux, EvaporationFlux, + SeaIceFreshWaterFlux, IceRunoffFlux, + RiverRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcThicknessForcing", 2); + } + if (CustomThicknessTend) { Pacer::start("Tend:customThicknessTend", 2); CustomThicknessTend(LocPseudoThicknessTend, State, AuxState, @@ -772,6 +814,7 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); + OMEGA_SCOPE(LocSfcTracerForcing, SfcTracerForcing); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -893,6 +936,42 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:surfaceTracerRestoring", 2); } + // compute tracer forcing tendency + if (LocSfcTracerForcing.Enabled) { + Pacer::start("Tend:sfcTracerForcing", 2); + const auto *ForcingState = Forcing::getDefault(); + const auto &LatentHeatFlux = + ForcingState->TracerForcing.LatentHeatFluxCell; + const auto &SensibleHeatFlux = + ForcingState->TracerForcing.SensibleHeatFluxCell; + const auto &LongWaveHeatFluxUp = + ForcingState->TracerForcing.LongWaveHeatFluxUpCell; + const auto &LongWaveHeatFluxDown = + ForcingState->TracerForcing.LongWaveHeatFluxDownCell; + const auto &SeaIceHeatFlux = + ForcingState->TracerForcing.SeaIceHeatFluxCell; + const auto &ShortWaveHeatFlux = + ForcingState->TracerForcing.ShortWaveHeatFluxCell; + const auto &SnowFlux = ForcingState->TracerForcing.SnowFluxCell; + const auto &RainFlux = ForcingState->TracerForcing.RainFluxCell; + const auto &IceRunoffFlux = ForcingState->TracerForcing.IceRunoffFluxCell; + const auto &RiverRunoffFlux = + ForcingState->TracerForcing.RiverRunoffFluxCell; + const auto &SeaIceSaltFlux = + ForcingState->TracerForcing.SeaIceSaltFluxCell; + const auto &PressureMid = VCoord->PressureMid; + + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + LocSfcTracerForcing( + LocTracerTend, ICell, TracerArray, PressureMid, LatentHeatFlux, + SensibleHeatFlux, LongWaveHeatFluxUp, LongWaveHeatFluxDown, + SeaIceHeatFlux, ShortWaveHeatFlux, SnowFlux, RainFlux, + IceRunoffFlux, RiverRunoffFlux, SeaIceSaltFlux); + }); + Pacer::stop("Tend:sfcTracerForcing", 2); + } + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 877ba22faf2b..c60a2783ddb3 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -74,6 +74,8 @@ class Tendencies { VelocityHyperDiffOnEdge VelocityHyperDiff; SfcStressForcingOnEdge SfcStressForcing; BottomDragOnEdge ExplicitBottomDrag; + SfcThicknessForcingOnCell SfcThicknessForcing; + SfcTracerForcingOnCell SfcTracerForcing; TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 5b142888142c..2353f38049bd 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -11,6 +11,7 @@ #include "TendencyTerms.h" #include "AuxiliaryState.h" #include "DataTypes.h" +#include "Eos.h" #include "HorzMesh.h" #include "HorzOperators.h" #include "OceanState.h" @@ -71,6 +72,19 @@ BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh, NVertLayers(VCoord->NVertLayers), EdgeMask(VCoord->EdgeMask), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +SfcThicknessForcingOnCell::SfcThicknessForcingOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + +SfcTracerForcingOnCell::SfcTracerForcingOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord, + I4 TempTracerIndex, + I4 SaltTracerIndex, + const Eos *EosInst) + : TempIndex(TempTracerIndex), SaltIndex(SaltTracerIndex), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + EosChoice(EosInst->EosChoice) {} + TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) : HorzontalMesh(Mesh), VerticalCoord(VCoord), diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 548832b290bf..a2a8dc09492c 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -11,6 +11,7 @@ //===----------------------------------------------------------------------===// #include "AuxiliaryState.h" +#include "Eos.h" #include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" @@ -372,6 +373,106 @@ class BottomDragOnEdge { Array1DI4 MaxLayerEdgeTop; }; +/// Coupled freshwater flux forcing for thickness equation. +class SfcThicknessForcingOnCell { + public: + bool Enabled = false; + + SfcThicknessForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 ICell, + const Array1DReal &SnowFlux, + const Array1DReal &RainFlux, + const Array1DReal &EvaporationFlux, + const Array1DReal &SeaIceFreshWaterFlux, + const Array1DReal &IceRunoffFlux, + const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { + + const I4 KTop = MinLayerCell(ICell); + if (KTop > MaxLayerCell(ICell)) { + return; + } + + const Real FreshWaterFlux = SnowFlux(ICell) + RainFlux(ICell) + + EvaporationFlux(ICell) + + SeaIceFreshWaterFlux(ICell) + + IceRunoffFlux(ICell) + RiverRunoffFlux(ICell); + + Tend(ICell, KTop) += (FreshWaterFlux + SeaIceSaltFlux(ICell)) / RhoSw; + } + + private: + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + +/// Coupled surface flux forcing for active tracers. +class SfcTracerForcingOnCell { + public: + bool Enabled = false; + + SfcTracerForcingOnCell(const HorzMesh *Mesh, const VertCoord *VCoord, + I4 TempTracerIndex, I4 SaltTracerIndex, + const Eos *EosInst); + + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 ICell, const Array3DReal &TracerCell, + const Array2DReal &PressureMid, const Array1DReal &LatentHeatFlux, + const Array1DReal &SensibleHeatFlux, + const Array1DReal &LongWaveHeatFluxUp, + const Array1DReal &LongWaveHeatFluxDown, + const Array1DReal &SeaIceHeatFlux, + const Array1DReal &ShortWaveHeatFlux, const Array1DReal &SnowFlux, + const Array1DReal &RainFlux, const Array1DReal &IceRunoffFlux, + const Array1DReal &RiverRunoffFlux, + const Array1DReal &SeaIceSaltFlux) const { + + const I4 KTop = MinLayerCell(ICell); + if (KTop > MaxLayerCell(ICell)) { + return; + } + + if (TempIndex >= 0) { + const Real PTopDb = PressureMid(ICell, KTop) * Pa2Db; + const Real SaTop = SaltIndex >= 0 + ? TracerCell(SaltIndex, ICell, KTop) + : 0.0_Real; // not sure we want zero here? + const Real CtFrz = + Eos::calcCtFreezing(EosChoice, SaTop, PTopDb, 0.0_Real); + const Real CtTop = TracerCell(TempIndex, ICell, KTop); + + // Heat tendencies are due to direct heat fluxes + enthalpy fluxes + // The enthalpy of liquid water is assumed to be: + // - local SST for liquid mass fluxes (rain, rivers) + // - local freezing point for solid --> liq mass fluxes (snow, frozen + // runoff) + // - solid mass fluxes are locally melted by the ocean (constant Lat + // heat of fusion) + const Real HeatFlux = + LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) + + LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) + + SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) + + (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop + + (SnowFlux(ICell) + IceRunoffFlux(ICell)) * + (Cp0Sw * CtFrz - LatIce); + + Tend(TempIndex, ICell, KTop) += HeatFlux * HFluxFac; + } + + if (SaltIndex >= 0) { + Tend(SaltIndex, ICell, KTop) += SeaIceSaltFlux(ICell) * SFluxFac; + } + } + + private: + I4 TempIndex; + I4 SaltIndex; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + EosType EosChoice; +}; + // Tracer horizontal advection term class TracerHorzAdvOnCell { public: diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp new file mode 100644 index 000000000000..a6478bb612c4 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.cpp @@ -0,0 +1,155 @@ +#include "TracerForcingVars.h" +#include "Eos.h" +#include "Field.h" +#include "Tracers.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +TracerForcingVars::TracerForcingVars(const std::string &Suffix, + const HorzMesh *Mesh) + : SnowFluxCell("SnowFlux" + Suffix, Mesh->NCellsSize), + RainFluxCell("RainFlux" + Suffix, Mesh->NCellsSize), + EvaporationFluxCell("EvaporationFlux" + Suffix, Mesh->NCellsSize), + SeaIceFreshWaterFluxCell("SeaIceFreshWaterFlux" + Suffix, + Mesh->NCellsSize), + IceRunoffFluxCell("IceRunoffFlux" + Suffix, Mesh->NCellsSize), + RiverRunoffFluxCell("RiverRunoffFlux" + Suffix, Mesh->NCellsSize), + LatentHeatFluxCell("LatentHeatFlux" + Suffix, Mesh->NCellsSize), + SensibleHeatFluxCell("SensibleHeatFlux" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxUpCell("LongWaveHeatFluxUp" + Suffix, Mesh->NCellsSize), + LongWaveHeatFluxDownCell("LongWaveHeatFluxDown" + Suffix, + Mesh->NCellsSize), + SeaIceHeatFluxCell("SeaIceHeatFlux" + Suffix, Mesh->NCellsSize), + ShortWaveHeatFluxCell("ShortWaveHeatFlux" + Suffix, Mesh->NCellsSize), + SeaIceSaltFluxCell("SeaIceSaltFlux" + Suffix, Mesh->NCellsSize) { + deepCopy(SnowFluxCell, 0.0_Real); + deepCopy(RainFluxCell, 0.0_Real); + deepCopy(EvaporationFluxCell, 0.0_Real); + deepCopy(SeaIceFreshWaterFluxCell, 0.0_Real); + deepCopy(IceRunoffFluxCell, 0.0_Real); + deepCopy(RiverRunoffFluxCell, 0.0_Real); + deepCopy(LatentHeatFluxCell, 0.0_Real); + deepCopy(SensibleHeatFluxCell, 0.0_Real); + deepCopy(LongWaveHeatFluxUpCell, 0.0_Real); + deepCopy(LongWaveHeatFluxDownCell, 0.0_Real); + deepCopy(SeaIceHeatFluxCell, 0.0_Real); + deepCopy(ShortWaveHeatFluxCell, 0.0_Real); + deepCopy(SeaIceSaltFluxCell, 0.0_Real); +} + +void TracerForcingVars::registerFields(const std::string &MeshName) const { + const int NDims = 1; + std::vector DimNames(NDims); + + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + DimNames[0] = "NCells" + DimSuffix; + + auto SnowFluxField = + Field::create(SnowFluxCell.label(), "snow freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto RainFluxField = + Field::create(RainFluxCell.label(), "rain freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto EvaporationFluxField = + Field::create(EvaporationFluxCell.label(), "evaporation freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto SeaIceFreshWaterFluxField = Field::create( + SeaIceFreshWaterFluxCell.label(), "sea-ice freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto IceRunoffFluxField = + Field::create(IceRunoffFluxCell.label(), "ice runoff freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto RiverRunoffFluxField = Field::create( + RiverRunoffFluxCell.label(), "river runoff freshwater flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + + auto LatentHeatFluxField = + Field::create(LatentHeatFluxCell.label(), "latent heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto SensibleHeatFluxField = + Field::create(SensibleHeatFluxCell.label(), "sensible heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto LongWaveHeatFluxUpField = Field::create( + LongWaveHeatFluxUpCell.label(), "upward longwave heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto LongWaveHeatFluxDownField = Field::create( + LongWaveHeatFluxDownCell.label(), "downward longwave heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto SeaIceHeatFluxField = + Field::create(SeaIceHeatFluxCell.label(), "sea-ice heat flux", "W m^-2", + "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + auto ShortWaveHeatFluxField = + Field::create(ShortWaveHeatFluxCell.label(), "shortwave heat flux", + "W m^-2", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + + auto SeaIceSaltFluxField = + Field::create(SeaIceSaltFluxCell.label(), "sea-ice salt flux", + "kg m^-2 s^-1", "", std::numeric_limits::lowest(), + std::numeric_limits::max(), NDims, DimNames); + + FieldGroup::addFieldToGroup(SnowFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(RainFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(EvaporationFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceFreshWaterFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(IceRunoffFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(RiverRunoffFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LatentHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SensibleHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LongWaveHeatFluxUpCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(LongWaveHeatFluxDownCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(ShortWaveHeatFluxCell.label(), "Forcing"); + FieldGroup::addFieldToGroup(SeaIceSaltFluxCell.label(), "Forcing"); + + SnowFluxField->attachData(SnowFluxCell); + RainFluxField->attachData(RainFluxCell); + EvaporationFluxField->attachData(EvaporationFluxCell); + SeaIceFreshWaterFluxField->attachData(SeaIceFreshWaterFluxCell); + IceRunoffFluxField->attachData(IceRunoffFluxCell); + RiverRunoffFluxField->attachData(RiverRunoffFluxCell); + LatentHeatFluxField->attachData(LatentHeatFluxCell); + SensibleHeatFluxField->attachData(SensibleHeatFluxCell); + LongWaveHeatFluxUpField->attachData(LongWaveHeatFluxUpCell); + LongWaveHeatFluxDownField->attachData(LongWaveHeatFluxDownCell); + SeaIceHeatFluxField->attachData(SeaIceHeatFluxCell); + ShortWaveHeatFluxField->attachData(ShortWaveHeatFluxCell); + SeaIceSaltFluxField->attachData(SeaIceSaltFluxCell); +} + +void TracerForcingVars::unregisterFields() const { + Field::destroy(SnowFluxCell.label()); + Field::destroy(RainFluxCell.label()); + Field::destroy(EvaporationFluxCell.label()); + Field::destroy(SeaIceFreshWaterFluxCell.label()); + Field::destroy(IceRunoffFluxCell.label()); + Field::destroy(RiverRunoffFluxCell.label()); + Field::destroy(LatentHeatFluxCell.label()); + Field::destroy(SensibleHeatFluxCell.label()); + Field::destroy(LongWaveHeatFluxUpCell.label()); + Field::destroy(LongWaveHeatFluxDownCell.label()); + Field::destroy(SeaIceHeatFluxCell.label()); + Field::destroy(ShortWaveHeatFluxCell.label()); + Field::destroy(SeaIceSaltFluxCell.label()); +} +} // namespace OMEGA diff --git a/components/omega/src/ocn/forcingVars/TracerForcingVars.h b/components/omega/src/ocn/forcingVars/TracerForcingVars.h new file mode 100644 index 000000000000..e38d9948f672 --- /dev/null +++ b/components/omega/src/ocn/forcingVars/TracerForcingVars.h @@ -0,0 +1,42 @@ +#ifndef OMEGA_TRACER_FORCING_H +#define OMEGA_TRACER_FORCING_H + +#include "DataTypes.h" +#include "HorzMesh.h" + +#include + +namespace OMEGA { + +// Forward declarations. Full definitions not needed in this header since only +// pointers are used. +class VertCoord; +class Eos; + +class TracerForcingVars { + public: + Array1DReal SnowFluxCell; + Array1DReal RainFluxCell; + Array1DReal EvaporationFluxCell; + Array1DReal SeaIceFreshWaterFluxCell; + Array1DReal IceRunoffFluxCell; + Array1DReal RiverRunoffFluxCell; + + Array1DReal LatentHeatFluxCell; + Array1DReal SensibleHeatFluxCell; + Array1DReal LongWaveHeatFluxUpCell; + Array1DReal LongWaveHeatFluxDownCell; + Array1DReal SeaIceHeatFluxCell; + Array1DReal ShortWaveHeatFluxCell; + + Array1DReal SeaIceSaltFluxCell; + + TracerForcingVars(const std::string &Suffix, const HorzMesh *Mesh); + + void registerFields(const std::string &MeshName) const; + void unregisterFields() const; +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 0e7eb2c5be76..f8ca084d4559 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -706,6 +706,42 @@ void testBruntVaisalaFreqSqTeos10() { return; } +/// Test all Eos::calcCtFreezing pathways (Teos10, Linear, Constant) +void testCalcCtFreezing() { + const Real RTol = 1e-10; + + constexpr Real SaturationFrac = 0.0; + constexpr Real PDb = 500.0; // pressure in dbar for GSW pathway + constexpr Real SaLocal = 32.0; + + const Real CtTeosExpected = + gsw_ct_freezing_poly(SaLocal, PDb, SaturationFrac); + const Real CtTeos = + Eos::calcCtFreezing(EosType::Teos10Eos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtTeos, CtTeosExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Teos10 FAIL, expected {}, got {}", + CtTeosExpected, CtTeos); + } + + const Real CtLinearExpected = -0.054_Real * SaLocal / Psu2Gpkg; + const Real CtLinear = + Eos::calcCtFreezing(EosType::LinearEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtLinear, CtLinearExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Linear FAIL, expected {}, got {}", + CtLinearExpected, CtLinear); + } + + const Real CtConstExpected = -1.9_Real; + const Real CtConst = + Eos::calcCtFreezing(EosType::ConstantEos, SaLocal, PDb, SaturationFrac); + if (!isApprox(CtConst, CtConstExpected, RTol)) { + ABORT_ERROR("testCalcCtFreezing: Constant FAIL, expected {}, got {}", + CtConstExpected, CtConst); + } + + return; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { Eos::destroyInstance(); @@ -767,22 +803,22 @@ void checkValueGswcN2() { } /// Test that the calcCtFreezing function returns the expected value -void checkValueCtFreezing() { +void checkValueGswcCtFreezing() { const Real RTol = 1e-10; - Teos10Eos TestEos(VertCoord::getDefault()); constexpr Real SaturationFrac = 0.0; constexpr Real P = 500.0 * Db2Pa; // Convert dbar to Pa constexpr Real Sa = 32.0; /// Get freezing temperature from GSW-C library double CtFreezGswc = gsw_ct_freezing_poly(Sa, P * Pa2Db, SaturationFrac); - double CtFreez = TestEos.calcCtFreezing(Sa, P * Pa2Db, SaturationFrac); + double CtFreez = + Teos10Eos::calcCtFreezingTeos10(Sa, P * Pa2Db, SaturationFrac); /// Check the value against the GSW-C value bool Check = isApprox(CtFreezGswc, CtFreez, RTol); if (!Check) { - ABORT_ERROR("checkValueCtFreezing: CtFreez FAIL, expected {}, got {}", + ABORT_ERROR("checkValueGswcCtFreezing: CtFreez FAIL, expected {}, got {}", CtFreezGswc, CtFreez); } return; @@ -837,7 +873,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); - checkValueCtFreezing(); + checkValueGswcCtFreezing(); checkValueGswcCtFromPt(); checkValueGswcPtFromCt(); @@ -848,6 +884,7 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { testEosTeos10(); testEosTeos10Displaced(); testBruntVaisalaFreqSqTeos10(); + testCalcCtFreezing(); finalizeEosTest(); diff --git a/components/omega/test/ocn/SfcCouplingTest.cpp b/components/omega/test/ocn/SfcCouplingTest.cpp index 2f20f41f67a6..94acc7e12a76 100644 --- a/components/omega/test/ocn/SfcCouplingTest.cpp +++ b/components/omega/test/ocn/SfcCouplingTest.cpp @@ -23,10 +23,15 @@ using namespace OMEGA; struct TestSetup { - std::map ImportIdxMap = {{"Foxx_taux", 3}, - {"Foxx_tauy", 8}}; + std::map ImportIdxMap = { + {"Foxx_taux", 3}, {"Foxx_tauy", 8}, {"Foxx_swnet", 0}, + {"Foxx_sen", 1}, {"Foxx_lat", 2}, {"Foxx_lwup", 4}, + {"Faxa_lwdn", 5}, {"Fioi_salt", 6}, {"Fioi_melth", 7}, + {"Fioi_meltw", 9}, {"Faxa_snow", 10}, {"Faxa_rain", 11}, + {"Foxx_evap", 12}, {"Foxx_rofl", 13}, {"Foxx_rofi", 14}}; std::map ExportIdxMap = { - {"So_t", 2}, {"So_s", 4}, {"So_u", 6}, {"So_v", 9}, {"So_ssh", 1}}; + {"So_t", 2}, {"So_s", 4}, {"So_u", 6}, {"So_v", 9}, + {"So_dhdx", 5}, {"So_dhdy", 3}, {"So_ssh", 1}}; }; CouplingInitParams mockCouplingInitParams( @@ -39,7 +44,7 @@ CouplingInitParams mockCouplingInitParams( TimeInterval CouplingTimeStep_ = CouplingTimeStep.value_or(DefTimeStepper->getTimeStep()); - CouplingInitParams CouplingParams{.NImportFields = 10, + CouplingInitParams CouplingParams{.NImportFields = 15, .NExportFields = 10, .ImportIdxMap = Setup.ImportIdxMap, .ExportIdxMap = Setup.ExportIdxMap, @@ -150,31 +155,64 @@ int testImportFromCoupler(const CouplingLayout Layout) { int NImports = DefCoupling->NImportFields; int NExports = DefCoupling->NExportFields; - int TauxIdx = CouplingParams.ImportIdxMap.at("Foxx_taux"); - int TauyIdx = CouplingParams.ImportIdxMap.at("Foxx_tauy"); - std::vector CplToOcnData(NCells * NImports, 0.0); std::vector OcnToCplData(NCells * NExports, 0.0); - HostArray1DReal ExpectedSfcStressZonal = - makeCellVarryingArray("ExpectedSfcStressZonal", NCells, Real(TauxIdx)); - HostArray1DReal ExpectedSfcStressMerid = - makeCellVarryingArray("ExpectedSfcStressMerid", NCells, Real(TauyIdx)); - - for (int Cell = 0; Cell < NCells; Cell++) { - CplToOcnData[flatIdx(Layout, Cell, TauxIdx, NCells, NImports)] = - ExpectedSfcStressZonal(Cell); - CplToOcnData[flatIdx(Layout, Cell, TauyIdx, NCells, NImports)] = - ExpectedSfcStressMerid(Cell); - } + auto fillImportField = [&](const std::string &Name) { + const int FieldIdx = CouplingParams.ImportIdxMap.at(Name); + for (int Cell = 0; Cell < NCells; Cell++) { + CplToOcnData[flatIdx(Layout, Cell, FieldIdx, NCells, NImports)] = + static_cast(FieldIdx + Cell); + } + }; + + fillImportField("Foxx_taux"); + fillImportField("Foxx_tauy"); + fillImportField("Foxx_swnet"); + fillImportField("Foxx_sen"); + fillImportField("Foxx_lat"); + fillImportField("Foxx_lwup"); + fillImportField("Faxa_lwdn"); + fillImportField("Fioi_salt"); + fillImportField("Fioi_melth"); + fillImportField("Fioi_meltw"); + fillImportField("Faxa_snow"); + fillImportField("Faxa_rain"); + fillImportField("Foxx_evap"); + fillImportField("Foxx_rofl"); + fillImportField("Foxx_rofi"); DefCoupling->attachData(CplToOcnData.data(), OcnToCplData.data()); DefCoupling->importFromCoupler(); - auto ImportPass = arraysEqual(DefCoupling->CplToOcn.SfcStressZonal, - ExpectedSfcStressZonal) && - arraysEqual(DefCoupling->CplToOcn.SfcStressMerid, - ExpectedSfcStressMerid); + auto checkImportField = [&](const HostArray1DReal &Field, + const std::string &Name) { + const int FieldIdx = CouplingParams.ImportIdxMap.at(Name); + HostArray1DReal Expected = + makeCellVarryingArray("Expected" + Name, NCells, Real(FieldIdx)); + return arraysEqual(Field, Expected); + }; + + auto ImportPass = + checkImportField(DefCoupling->CplToOcn.SfcStressZonal, "Foxx_taux") && + checkImportField(DefCoupling->CplToOcn.SfcStressMerid, "Foxx_tauy") && + checkImportField(DefCoupling->CplToOcn.ShortWaveHeatFlux, + "Foxx_swnet") && + checkImportField(DefCoupling->CplToOcn.SensibleHeatFlux, "Foxx_sen") && + checkImportField(DefCoupling->CplToOcn.LatentHeatFlux, "Foxx_lat") && + checkImportField(DefCoupling->CplToOcn.LongWaveHeatFluxUp, + "Foxx_lwup") && + checkImportField(DefCoupling->CplToOcn.LongWaveHeatFluxDown, + "Faxa_lwdn") && + checkImportField(DefCoupling->CplToOcn.SeaIceSaltFlux, "Fioi_salt") && + checkImportField(DefCoupling->CplToOcn.SeaIceHeatFlux, "Fioi_melth") && + checkImportField(DefCoupling->CplToOcn.SeaIceFreshWaterFlux, + "Fioi_meltw") && + checkImportField(DefCoupling->CplToOcn.SnowFlux, "Faxa_snow") && + checkImportField(DefCoupling->CplToOcn.RainFlux, "Faxa_rain") && + checkImportField(DefCoupling->CplToOcn.EvaporationFlux, "Foxx_evap") && + checkImportField(DefCoupling->CplToOcn.RiverRunoffFlux, "Foxx_rofl") && + checkImportField(DefCoupling->CplToOcn.IceRunoffFlux, "Foxx_rofi"); if (ImportPass) { LOG_INFO("SfcCouplingTest: importFromCoupler with {} layout PASS", @@ -204,31 +242,74 @@ int testApplyImportFields() { int NImports = DefCoupling->NImportFields; int NExports = DefCoupling->NExportFields; - int TauxIdx = CouplingParams.ImportIdxMap.at("Foxx_taux"); - int TauyIdx = CouplingParams.ImportIdxMap.at("Foxx_tauy"); - std::vector CplToOcnData(NCells * NImports, 0.0); std::vector OcnToCplData(NCells * NExports, 0.0); - int Offset = 27; - HostArray1DReal ExpectedSfcStressZonal = makeCellVarryingArray( - "ExpectedSfcStressZonal", NCells, Real(TauxIdx + Offset)); - HostArray1DReal ExpectedSfcStressMerid = makeCellVarryingArray( - "ExpectedSfcStressMerid", NCells, Real(TauyIdx + Offset)); - - // Copy the expected values into the CplToOcn fields directly - deepCopy(DefCoupling->CplToOcn.SfcStressZonal, ExpectedSfcStressZonal); - deepCopy(DefCoupling->CplToOcn.SfcStressMerid, ExpectedSfcStressMerid); + int Offset = 27; + + auto setImportField = [&](HostArray1DReal &Field, const std::string &Name) { + const int FieldIdx = CouplingParams.ImportIdxMap.at(Name); + HostArray1DReal Expected = makeCellVarryingArray( + "Expected" + Name, NCells, Real(FieldIdx + Offset)); + deepCopy(Field, Expected); + }; + + setImportField(DefCoupling->CplToOcn.SfcStressZonal, "Foxx_taux"); + setImportField(DefCoupling->CplToOcn.SfcStressMerid, "Foxx_tauy"); + setImportField(DefCoupling->CplToOcn.ShortWaveHeatFlux, "Foxx_swnet"); + setImportField(DefCoupling->CplToOcn.SensibleHeatFlux, "Foxx_sen"); + setImportField(DefCoupling->CplToOcn.LatentHeatFlux, "Foxx_lat"); + setImportField(DefCoupling->CplToOcn.LongWaveHeatFluxUp, "Foxx_lwup"); + setImportField(DefCoupling->CplToOcn.LongWaveHeatFluxDown, "Faxa_lwdn"); + setImportField(DefCoupling->CplToOcn.SeaIceSaltFlux, "Fioi_salt"); + setImportField(DefCoupling->CplToOcn.SeaIceHeatFlux, "Fioi_melth"); + setImportField(DefCoupling->CplToOcn.SeaIceFreshWaterFlux, "Fioi_meltw"); + setImportField(DefCoupling->CplToOcn.SnowFlux, "Faxa_snow"); + setImportField(DefCoupling->CplToOcn.RainFlux, "Faxa_rain"); + setImportField(DefCoupling->CplToOcn.EvaporationFlux, "Foxx_evap"); + setImportField(DefCoupling->CplToOcn.RiverRunoffFlux, "Foxx_rofl"); + setImportField(DefCoupling->CplToOcn.IceRunoffFlux, "Foxx_rofi"); DefCoupling->applyImportFields(DefForcing); - auto SfcStressZonalOwned = Kokkos::subview( - DefForcing->SfcStressForcing.ZonalStressCell, std::pair(0, NCells)); - auto SfcStressMeridOwned = Kokkos::subview( - DefForcing->SfcStressForcing.MeridStressCell, std::pair(0, NCells)); - - auto ApplyPass = arraysEqual(SfcStressZonalOwned, ExpectedSfcStressZonal) && - arraysEqual(SfcStressMeridOwned, ExpectedSfcStressMerid); + auto checkAppliedField = [&](const Array1DReal &Field, + const std::string &Name) { + const int FieldIdx = CouplingParams.ImportIdxMap.at(Name); + HostArray1DReal Expected = makeCellVarryingArray( + "Expected" + Name, NCells, Real(FieldIdx + Offset)); + auto Owned = Kokkos::subview(Field, std::pair(0, NCells)); + return arraysEqual(Owned, Expected); + }; + + auto ApplyPass = + checkAppliedField(DefForcing->SfcStressForcing.ZonalStressCell, + "Foxx_taux") && + checkAppliedField(DefForcing->SfcStressForcing.MeridStressCell, + "Foxx_tauy") && + checkAppliedField(DefForcing->TracerForcing.ShortWaveHeatFluxCell, + "Foxx_swnet") && + checkAppliedField(DefForcing->TracerForcing.SensibleHeatFluxCell, + "Foxx_sen") && + checkAppliedField(DefForcing->TracerForcing.LatentHeatFluxCell, + "Foxx_lat") && + checkAppliedField(DefForcing->TracerForcing.LongWaveHeatFluxUpCell, + "Foxx_lwup") && + checkAppliedField(DefForcing->TracerForcing.LongWaveHeatFluxDownCell, + "Faxa_lwdn") && + checkAppliedField(DefForcing->TracerForcing.SeaIceSaltFluxCell, + "Fioi_salt") && + checkAppliedField(DefForcing->TracerForcing.SeaIceHeatFluxCell, + "Fioi_melth") && + checkAppliedField(DefForcing->TracerForcing.SeaIceFreshWaterFluxCell, + "Fioi_meltw") && + checkAppliedField(DefForcing->TracerForcing.SnowFluxCell, "Faxa_snow") && + checkAppliedField(DefForcing->TracerForcing.RainFluxCell, "Faxa_rain") && + checkAppliedField(DefForcing->TracerForcing.EvaporationFluxCell, + "Foxx_evap") && + checkAppliedField(DefForcing->TracerForcing.RiverRunoffFluxCell, + "Foxx_rofl") && + checkAppliedField(DefForcing->TracerForcing.IceRunoffFluxCell, + "Foxx_rofi"); if (ApplyPass) { LOG_INFO("SfcCouplingTest: applyImportFields PASS"); @@ -246,11 +327,7 @@ int testUpdateExportFields(const I4 NSteps) { int Err = 0; - auto *DefStepper = TimeStepper::getDefault(); - Clock *ModelClock = DefStepper->getClock(); - - // Reset the shared clock - ModelClock->setCurrentTime(DefStepper->getStartTime()); + auto *DefStepper = TimeStepper::getDefault(); // Coupling interval spans NSteps ocean timesteps auto CouplingParams = mockCouplingInitParams( @@ -271,7 +348,10 @@ int testUpdateExportFields(const I4 NSteps) { Tracers::getIndex(TempIdx, "Temperature"); Tracers::getIndex(SalinIdx, "Salinity"); - while (!DefCoupling->getCouplingAlarm()->isRinging()) { + // Run exactly NSteps updates. Clock and alarm behavior is covered by + // TimeMgrTest; keeping this test focused on field accumulation avoids + // advancing a shared clock with alarms from cleared test objects. + for (I4 Step = 0; Step < NSteps; ++Step) { Real CurrStep = static_cast(DefCoupling->getNAccumSteps()); HostArray2DReal TempH = Tracers::getHostByIndex(0, TempIdx); @@ -286,11 +366,9 @@ int testUpdateExportFields(const I4 NSteps) { Tracers::copyToDevice(0); DefCoupling->updateExportFields(DefState, Tracers::getAll(0)); - - ModelClock->advance(); } - // Sanity check: alarm should ring after NSteps + // Sanity check: the expected number of updates was performed if (DefCoupling->getNAccumSteps() != NSteps) { Err++; LOG_ERROR("SfcCouplingTest: updateExportFields FAIL - " @@ -322,7 +400,7 @@ int testUpdateExportFields(const I4 NSteps) { TempErr++; } - if (!isApprox(AvgSalinH(Cell), ExpectedSalin(Cell), RTol)) { + if (!isApprox(AvgSalinH(Cell), ExpectedSalin(Cell) / Psu2Gpkg, RTol)) { SalinErr++; } } @@ -351,9 +429,6 @@ int testUpdateExportFields(const I4 NSteps) { RTol, SalinErr); } - // reset model clock to the start time for any subsequent tests - ModelClock->setCurrentTime(DefStepper->getStartTime()); - SfcCoupling::clear(); return Err; } @@ -413,11 +488,14 @@ int testExportToCoupler(const CouplingLayout Layout) { DefCoupling->exportToCoupler(); - // Check 1: exportToCoupler properly packs into OcnToCplView. Velocity - // is skipped here: its averaging is a hardcoded stub pending real vector - // reconstruction (see OcnToCplFields::updateAverages), not yet - // meaningful to check. - // copyToHost() converts temp to Kelvin (identity CT->PT w/ ConstantEos) + // Check 1: exportToCoupler properly packs into OcnToCplView. + // NormalVelocity is accumulated on edges and reconstructed at cell centers + // during copyToHost(). Recon correctness is tested in HorzOperatorsTest. + // copyToHost() converts conservative temperature to in-situ Kelvin and + // absolute salinity to practical salinity. + + // TODO: Add end-to-end SSH-gradient accumulation/export coverage if we + // decide to continue passing ssh grad (cf. ssh directly) to mpas-si int PackErr = 0; for (int Cell = 0; Cell < NCells; Cell++) { if (OcnToCplData[flatIdx(Layout, Cell, TempIdx, NCells, NExports)] != @@ -425,7 +503,7 @@ int testExportToCoupler(const CouplingLayout Layout) { PackErr++; } if (OcnToCplData[flatIdx(Layout, Cell, SalinIdx, NCells, NExports)] != - ExpectedSalin(Cell)) { + ExpectedSalin(Cell) / Psu2Gpkg) { PackErr++; } if (OcnToCplData[flatIdx(Layout, Cell, SshIdx, NCells, NExports)] != @@ -496,7 +574,7 @@ int testEraseAndGet() { TimeInterval TimeStep = DefStepper->getTimeStep(); // test creation of a non-default, named surface coupling object - SfcCoupling::create("AnotherSfcCoupling", DefMesh, 10, 12, + SfcCoupling::create("AnotherSfcCoupling", DefMesh, 15, 12, Setup.ImportIdxMap, Setup.ExportIdxMap, DefStepper, TimeStep, CouplingLayout::MCT); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5268d0a436ea..13f8ed2ba2a4 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -54,6 +54,9 @@ struct TestSetup { constexpr Geometry Geom = Geometry::Spherical; constexpr int NVertLayers = 60; +int testSfcTracerForcing(); +int testSfcThicknessForcing(); + int initState() { int Err = 0; @@ -303,11 +306,26 @@ int testTendencies() { "NormalVelocityTend"); } + const Real NormVelTendSum = + sum(DefTendencies->NormalVelocityTend, Mesh->NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); + if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { + Err++; + LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); + } + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + // Test surface tracer forcing with enthalpy terms (TEOS-10 CtFrz path) + const int TracerForcingErr = testSfcTracerForcing(); + Err += TracerForcingErr; + + // Test surface thickness forcing with freshwater terms + const int ThicknessForcingErr = testSfcThicknessForcing(); + Err += ThicknessForcingErr; + // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; - int NEdgesOwned = Mesh->NEdgesOwned; int NTracers = Tracers::getNumTracers(); const Real PseudoThickTendSum = @@ -318,14 +336,6 @@ int testTendencies() { LOG_ERROR("TendenciesTest: PseudoThickTend FAIL"); } - const Real NormVelTendSum = - sum(DefTendencies->NormalVelocityTend, NEdgesOwned, - VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); - if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { - Err++; - LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); - } - const Real TraceTendSum = sum(DefTendencies->TracerTend, NTracers, NCellsOwned, VCoord->MinLayerCell, VCoord->MaxLayerCell); @@ -335,6 +345,382 @@ int testTendencies() { } Tendencies::clear(); + return Err; +} + +int testSfcTracerForcing() { + int Err = 0; + + auto *VCoord = VertCoord::getDefault(); + auto *DefTendencies = Tendencies::getDefault(); + auto *State = OceanState::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + auto *DefForcing = Forcing::getDefault(); + auto *EosInst = Eos::getInstance(); + + Array3DReal TracerArray = Tracers::getAll(0); + + const I4 TempIndex = Tracers::IndxTemp; + const I4 SaltIndex = Tracers::IndxSalt; + + if (TempIndex < 0 || SaltIndex < 0) { + LOG_ERROR("TendenciesTest: Invalid tracer indices for SfcTracerForcing"); + return -1; + } + + deepCopy(DefTendencies->TracerTend, 0._Real); + + // Set up single test cell at top layer + const I4 ICellTest = 0; + const I4 KTop = VCoord->MinLayerCellH(ICellTest); + + if (KTop > VCoord->MaxLayerCellH(ICellTest)) { + LOG_ERROR("TendenciesTest: Test cell has no layers"); + return -1; + } + + // Known tracer values for testing + const Real CtTopValue = 15.0_Real; // °C (conservative temperature) + const Real SaTopValue = 35.0_Real; // g/kg (salinity) + + // Set tracer values at test cell + OMEGA_SCOPE(LocTracerArray, TracerArray); + Kokkos::parallel_for( + "SetTestTracersForcing", 1, KOKKOS_LAMBDA(int i) { + LocTracerArray(TempIndex, ICellTest, KTop) = CtTopValue; + LocTracerArray(SaltIndex, ICellTest, KTop) = SaTopValue; + }); + + // Retrieve forcing field views + auto &SensibleHeatFlux = DefForcing->TracerForcing.SensibleHeatFluxCell; + auto &LatentHeatFlux = DefForcing->TracerForcing.LatentHeatFluxCell; + auto &LongWaveHeatFluxUp = DefForcing->TracerForcing.LongWaveHeatFluxUpCell; + auto &LongWaveHeatFluxDown = + DefForcing->TracerForcing.LongWaveHeatFluxDownCell; + auto &SeaIceHeatFlux = DefForcing->TracerForcing.SeaIceHeatFluxCell; + auto &ShortWaveHeatFlux = DefForcing->TracerForcing.ShortWaveHeatFluxCell; + auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; + auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; + auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; + auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; + auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; + + // Initialize all fluxes to zero + deepCopy(SensibleHeatFlux, 0._Real); + deepCopy(LatentHeatFlux, 0._Real); + deepCopy(LongWaveHeatFluxUp, 0._Real); + deepCopy(LongWaveHeatFluxDown, 0._Real); + deepCopy(SeaIceHeatFlux, 0._Real); + deepCopy(ShortWaveHeatFlux, 0._Real); + deepCopy(RainFlux, 0._Real); + deepCopy(RiverRunoffFlux, 0._Real); + deepCopy(SnowFlux, 0._Real); + deepCopy(IceRunoffFlux, 0._Real); + deepCopy(SeaIceSaltFlux, 0._Real); + + // Set test forcing values + // Non-zero sensible heat: 100 W/m² + const Real TestSensibleHeat = 100.0_Real; + // Non-zero rain: 1e-8 kg/m²/s + const Real TestRain = 1.0e-8_Real; + // Non-zero snow: 5e-9 kg/m²/s + const Real TestSnow = 5.0e-9_Real; + // Sea ice salt flux: 1e-4 kg/m²/s + const Real TestSeaIceSaltFlux = 1.0e-4_Real; + + OMEGA_SCOPE(LocSensibleHeatFlux, SensibleHeatFlux); + OMEGA_SCOPE(LocRainFlux, RainFlux); + OMEGA_SCOPE(LocSnowFlux, SnowFlux); + OMEGA_SCOPE(LocSeaIceSaltFlux, SeaIceSaltFlux); + Kokkos::parallel_for( + "SetTestForcingTracer", 1, KOKKOS_LAMBDA(int i) { + LocSensibleHeatFlux(ICellTest) = TestSensibleHeat; + LocRainFlux(ICellTest) = TestRain; + LocSnowFlux(ICellTest) = TestSnow; + LocSeaIceSaltFlux(ICellTest) = TestSeaIceSaltFlux; + }); + + DefForcing->computeAll(); + + // Disable all tendencies except SfcTracerForcing + const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; + const bool OrigSfcThicknessEnabled = + DefTendencies->SfcThicknessForcing.Enabled; + const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; + const bool OrigPseudoThicknessDiv = + DefTendencies->PseudoThicknessFluxDiv.Enabled; + const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; + const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; + const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; + const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; + const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; + const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; + const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; + const bool OrigSurfaceTracerRestoring = + DefTendencies->SurfaceTracerRestoring.Enabled; + + DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->SfcThicknessForcing.Enabled = false; + DefTendencies->SfcTracerForcing.Enabled = false; + DefTendencies->PseudoThicknessFluxDiv.Enabled = false; + DefTendencies->PotentialVortHAdv.Enabled = false; + DefTendencies->KEGrad.Enabled = false; + DefTendencies->VelocityDiffusion.Enabled = false; + DefTendencies->VelocityHyperDiff.Enabled = false; + DefTendencies->TracerHorzAdv.Enabled = false; + DefTendencies->TracerDiffusion.Enabled = false; + DefTendencies->TracerHyperDiff.Enabled = false; + DefTendencies->SurfaceTracerRestoring.Enabled = false; + + // Compute tendencies + int ThickTimeLevel = 0; + int VelTimeLevel = 0; + int TracerTimeLevel = 0; + TimeInstant Time; + TimeInterval Interval(1., TimeUnits::Seconds); + + // because vertical advection tendencies are always on, we need to compute a + // baseline first. the actual test is whether the total tendencies change + // with the flag toggling. + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + HostArray3DReal TracerTendBaseH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendBaseH, DefTendencies->TracerTend); + const Real BaselineTempTend = TracerTendBaseH(TempIndex, ICellTest, KTop); + const Real BaselineSaltTend = TracerTendBaseH(SaltIndex, ICellTest, KTop); + + // Now enable SfcTracerForcing and compute again + DefTendencies->SfcTracerForcing.Enabled = true; + + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + // Build a reference expectations for temperature tendency: + // using TEOS-10 freezing CT (expected to pass under strict tolerance). + + HostArray2DReal PressureMidH = createHostMirrorCopy(VCoord->PressureMid); + deepCopy(PressureMidH, VCoord->PressureMid); + const Real PTopDb = PressureMidH(ICellTest, KTop) * Pa2Db; + const Real CtFrz = + Eos::calcCtFreezing(EosInst->EosChoice, SaTopValue, PTopDb, 0.0_Real); + const Real ExpectedTempTend = + (TestSensibleHeat + TestRain * Cp0Sw * CtTopValue + + TestSnow * (Cp0Sw * CtFrz - LatIce)) * + HFluxFac; + + // SaltTend = SeaIceSaltFlux * SFluxFac + const Real ExpectedSaltTend = TestSeaIceSaltFlux * SFluxFac; + + HostArray3DReal TracerTendH = + createHostMirrorCopy(DefTendencies->TracerTend); + deepCopy(TracerTendH, DefTendencies->TracerTend); + const Real ComputedTempTend = + TracerTendH(TempIndex, ICellTest, KTop) - BaselineTempTend; + const Real ComputedSaltTend = + TracerTendH(SaltIndex, ICellTest, KTop) - BaselineSaltTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; // flux precision is ~e-15 + + // Expected-pass check with TEOS freezing CT reference. + if (!isApprox(ComputedTempTend, ExpectedTempTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing temp tendency FAIL"); + LOG_ERROR(" with TEOS-CtFrz Expected: {}, Computed: {}, Diff: {}", + ExpectedTempTend, ComputedTempTend, + Kokkos::abs(ComputedTempTend - ExpectedTempTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing temp tendency PASS"); + } + + // Check salinity tendency + if (!isApprox(ComputedSaltTend, ExpectedSaltTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcTracerForcing salt tendency FAIL"); + LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedSaltTend, + ComputedSaltTend, + Kokkos::abs(ComputedSaltTend - ExpectedSaltTend)); + } else { + LOG_INFO("TendenciesTest: SfcTracerForcing salt tendency PASS"); + } + + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; + DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; + DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; + DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; + DefTendencies->KEGrad.Enabled = OrigKEGrad; + DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; + DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; + DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; + DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; + DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; + DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; + + return Err; +} + +int testSfcThicknessForcing() { + int Err = 0; + + auto *VCoord = VertCoord::getDefault(); + auto *DefTendencies = Tendencies::getDefault(); + auto *State = OceanState::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + auto *DefForcing = Forcing::getDefault(); + + Array3DReal TracerArray = Tracers::getAll(0); + + deepCopy(DefTendencies->PseudoThicknessTend, 0._Real); + + // Set up single test cell at top layer + const I4 ICellTest = 0; + const I4 KTop = VCoord->MinLayerCellH(ICellTest); + + if (KTop > VCoord->MaxLayerCellH(ICellTest)) { + LOG_ERROR("TendenciesTest: Test cell has no layers for thickness test"); + return -1; + } + + // Retrieve forcing field views for thickness + auto &SnowFlux = DefForcing->TracerForcing.SnowFluxCell; + auto &RainFlux = DefForcing->TracerForcing.RainFluxCell; + auto &EvaporationFlux = DefForcing->TracerForcing.EvaporationFluxCell; + auto &SeaIceFreshWater = DefForcing->TracerForcing.SeaIceFreshWaterFluxCell; + auto &IceRunoffFlux = DefForcing->TracerForcing.IceRunoffFluxCell; + auto &RiverRunoffFlux = DefForcing->TracerForcing.RiverRunoffFluxCell; + auto &SeaIceSaltFlux = DefForcing->TracerForcing.SeaIceSaltFluxCell; + + // Initialize all fluxes to zero + deepCopy(SnowFlux, 0._Real); + deepCopy(RainFlux, 0._Real); + deepCopy(EvaporationFlux, 0._Real); + deepCopy(SeaIceFreshWater, 0._Real); + deepCopy(IceRunoffFlux, 0._Real); + deepCopy(RiverRunoffFlux, 0._Real); + deepCopy(SeaIceSaltFlux, 0._Real); + + // Set test freshwater flux values + // Rain: 1e-8 kg/m²/s + const Real TestRain = 1.0e-8_Real; + // Snow: 5e-9 kg/m²/s + const Real TestSnow = 5.0e-9_Real; + // Ice runoff: 2e-9 kg/m²/s + const Real TestIceRunoff = 2.0e-9_Real; + // River runoff: 3e-9 kg/m²/s + const Real TestRiverRunoff = 3.0e-9_Real; + // Sea ice freshwater: 1e-9 kg/m²/s + const Real TestSeaIceFreshWater = 1.0e-9_Real; + // Sea ice salt flux: 1e-4 kg/m²/s (affects thickness via salt) + const Real TestSeaIceSaltFlux = 1.0e-4_Real; + + OMEGA_SCOPE(LocSnowFlux, SnowFlux); + OMEGA_SCOPE(LocRainFlux, RainFlux); + OMEGA_SCOPE(LocIceRunoffFlux, IceRunoffFlux); + OMEGA_SCOPE(LocRiverRunoffFlux, RiverRunoffFlux); + OMEGA_SCOPE(LocSeaIceFreshWater, SeaIceFreshWater); + OMEGA_SCOPE(LocSeaIceSaltFlux, SeaIceSaltFlux); + Kokkos::parallel_for( + "SetTestForcingThickness", 1, KOKKOS_LAMBDA(int i) { + LocRainFlux(ICellTest) = TestRain; + LocSnowFlux(ICellTest) = TestSnow; + LocIceRunoffFlux(ICellTest) = TestIceRunoff; + LocRiverRunoffFlux(ICellTest) = TestRiverRunoff; + LocSeaIceFreshWater(ICellTest) = TestSeaIceFreshWater; + LocSeaIceSaltFlux(ICellTest) = TestSeaIceSaltFlux; + }); + DefForcing->computeAll(); + + const bool OrigSfcStressEnabled = DefTendencies->SfcStressForcing.Enabled; + const bool OrigSfcThicknessEnabled = + DefTendencies->SfcThicknessForcing.Enabled; + const bool OrigSfcTracerEnabled = DefTendencies->SfcTracerForcing.Enabled; + const bool OrigPseudoThicknessDiv = + DefTendencies->PseudoThicknessFluxDiv.Enabled; + const bool OrigPotentialVortHAdv = DefTendencies->PotentialVortHAdv.Enabled; + const bool OrigKEGrad = DefTendencies->KEGrad.Enabled; + const bool OrigVelocityDiffusion = DefTendencies->VelocityDiffusion.Enabled; + const bool OrigVelocityHyperDiff = DefTendencies->VelocityHyperDiff.Enabled; + const bool OrigTracerHorzAdv = DefTendencies->TracerHorzAdv.Enabled; + const bool OrigTracerDiffusion = DefTendencies->TracerDiffusion.Enabled; + const bool OrigTracerHyperDiff = DefTendencies->TracerHyperDiff.Enabled; + const bool OrigSurfaceTracerRestoring = + DefTendencies->SurfaceTracerRestoring.Enabled; + + DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->SfcThicknessForcing.Enabled = false; + DefTendencies->SfcTracerForcing.Enabled = false; + DefTendencies->PseudoThicknessFluxDiv.Enabled = false; + DefTendencies->PotentialVortHAdv.Enabled = false; + DefTendencies->KEGrad.Enabled = false; + DefTendencies->VelocityDiffusion.Enabled = false; + DefTendencies->VelocityHyperDiff.Enabled = false; + DefTendencies->TracerHorzAdv.Enabled = false; + DefTendencies->TracerDiffusion.Enabled = false; + DefTendencies->TracerHyperDiff.Enabled = false; + DefTendencies->SurfaceTracerRestoring.Enabled = false; + + // Compute baseline tendencies (vertical advection is always on) + int ThickTimeLevel = 0; + int VelTimeLevel = 0; + TimeInstant Time; + DefTendencies->computePseudoThicknessTendenciesOnly( + State, AuxState, ThickTimeLevel, VelTimeLevel, Time); + + HostArray2DReal PseudoThicknessTendBaseH = + createHostMirrorCopy(DefTendencies->PseudoThicknessTend); + deepCopy(PseudoThicknessTendBaseH, DefTendencies->PseudoThicknessTend); + const Real BaselineThickTend = PseudoThicknessTendBaseH(ICellTest, KTop); + + // Now enable SfcThicknessForcing and compute again + DefTendencies->SfcThicknessForcing.Enabled = true; + DefTendencies->computePseudoThicknessTendenciesOnly( + State, AuxState, ThickTimeLevel, VelTimeLevel, Time); + + // Calculate expected thickness tendency + // ThickTend = (Rain + Snow + IceRunoff + RiverRunoff + SeaIceFreshWater + + // SeaIceSaltFlux) / RhoSw + const Real ExpectedThickTend = + (TestRain + TestSnow + TestIceRunoff + TestRiverRunoff + + TestSeaIceFreshWater + TestSeaIceSaltFlux) / + RhoSw; + + HostArray2DReal PseudoThicknessTendH = + createHostMirrorCopy(DefTendencies->PseudoThicknessTend); + deepCopy(PseudoThicknessTendH, DefTendencies->PseudoThicknessTend); + const Real ComputedThickTend = + PseudoThicknessTendH(ICellTest, KTop) - BaselineThickTend; + + constexpr Real RelTol = 1.0e-10_Real; + constexpr Real AbsTol = 1.0e-12_Real; + + // Check thickness tendency + if (!isApprox(ComputedThickTend, ExpectedThickTend, RelTol, AbsTol)) { + Err++; + LOG_ERROR("TendenciesTest: SfcThicknessForcing thickness tendency FAIL"); + LOG_INFO(" Expected: {}, Computed: {}, Diff: {}", ExpectedThickTend, + ComputedThickTend, + Kokkos::abs(ComputedThickTend - ExpectedThickTend)); + } else { + LOG_INFO("TendenciesTest: SfcThicknessForcing thickness tendency PASS"); + } + + DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + DefTendencies->SfcThicknessForcing.Enabled = OrigSfcThicknessEnabled; + DefTendencies->SfcTracerForcing.Enabled = OrigSfcTracerEnabled; + DefTendencies->PseudoThicknessFluxDiv.Enabled = OrigPseudoThicknessDiv; + DefTendencies->PotentialVortHAdv.Enabled = OrigPotentialVortHAdv; + DefTendencies->KEGrad.Enabled = OrigKEGrad; + DefTendencies->VelocityDiffusion.Enabled = OrigVelocityDiffusion; + DefTendencies->VelocityHyperDiff.Enabled = OrigVelocityHyperDiff; + DefTendencies->TracerHorzAdv.Enabled = OrigTracerHorzAdv; + DefTendencies->TracerDiffusion.Enabled = OrigTracerDiffusion; + DefTendencies->TracerHyperDiff.Enabled = OrigTracerHyperDiff; + DefTendencies->SurfaceTracerRestoring.Enabled = OrigSurfaceTracerRestoring; return Err; } diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index af84424217ff..d64b461896a5 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -255,6 +255,8 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerDiffusion.Enabled = false; TestTendencies->TracerHyperDiff.Enabled = false; TestTendencies->SfcStressForcing.Enabled = false; + TestTendencies->SfcTracerForcing.Enabled = false; + TestTendencies->SfcThicknessForcing.Enabled = false; TestTendencies->SurfaceTracerRestoring.Enabled = false; TestTendencies->ExplicitBottomDrag.Enabled = false; DefVAdv->ThickVertAdvEnabled = false; diff --git a/driver-mct/cime_config/config_component_e3sm.xml b/driver-mct/cime_config/config_component_e3sm.xml index 28ec8bb61c4f..7a148d9e6e70 100755 --- a/driver-mct/cime_config/config_component_e3sm.xml +++ b/driver-mct/cime_config/config_component_e3sm.xml @@ -617,6 +617,7 @@ 6 $ATM_NCPL $ATM_NCPL + $ATM_NCPL $ATM_NCPL $ATM_NCPL $ATM_NCPL diff --git a/driver-moab/cime_config/config_component_e3sm.xml b/driver-moab/cime_config/config_component_e3sm.xml index 755cf13e5650..e7c59d6897de 100644 --- a/driver-moab/cime_config/config_component_e3sm.xml +++ b/driver-moab/cime_config/config_component_e3sm.xml @@ -626,6 +626,7 @@ 6 $ATM_NCPL $ATM_NCPL + $ATM_NCPL $ATM_NCPL $ATM_NCPL $ATM_NCPL