diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bdeaddc..97f3a5e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,6 +32,7 @@ SPDX-License-Identifier: CC-BY-4.0 * Grids built from `--elmer-mesh` without a DEM now place each column at the latitude and longitude of its mesh cell center, as they already did with `--netcdf-dem-mesh`. Every column previously sat at 75 N 40 W regardless of the mesh. Results of such runs change through the solar geometry, unless the surface energy balance is taken from ICON-Land (`--couple-to-icon-land`). The `lon` written to the output now ranges from -180 to 180 degrees east instead of being 320. https://github.com/EBFMorg/EBFM/pull/193 * Grids built from `--elmer-mesh` alone (no `--netcdf-dem-mesh`) now take each column's elevation from the mesh: the mean of the z coordinates of its cell's vertices, as `reader.py` writes them from a DEM. The elevation was previously random between 0 and 100 m. A mesh without elevation data, such as `examples/greenland_mesh_v0/MESH`, now gives 0 m everywhere. https://github.com/EBFMorg/EBFM/pull/193 * Grids built from `--elmer-mesh` alone (no `--netcdf-dem-mesh`) now take each column's surface slope and orientation from the plane through the vertices of its mesh cell, instead of treating every column as flat. This changes the shortwave radiation on sloping columns. A mesh without elevation data still gives flat columns. When coupled to Elmer/Ice (`--couple-to-elmer-ice`), the columns stay flat and a warning is logged, since the slopes do not follow the surface elevation received from Elmer/Ice yet. https://github.com/EBFMorg/EBFM/pull/193 +* Add options to choose the path (`--output-file`, default `Output/model_output.nc`) and the frequency (`--output-frequency`, an ISO 8601 duration, default `P1D`) of the NetCDF output file. https://github.com/EBFMorg/EBFM/pull/211 # v0.8.0 diff --git a/src/ebfm/core/INIT.py b/src/ebfm/core/INIT.py index e25c394f..63dd25c6 100644 --- a/src/ebfm/core/INIT.py +++ b/src/ebfm/core/INIT.py @@ -11,12 +11,12 @@ from netCDF4 import Dataset, num2date from pathlib import Path -from datetime import datetime +from datetime import datetime, timedelta from ebfm.reader import read_elmer_mesh, read_dem, read_dem_xios from ebfm.elmer.mesh import Mesh -from .config import TimeConfig, GridConfig, ColumnDiscretizationConfig, iso8601 +from .config import TimeConfig, GridConfig, ColumnDiscretizationConfig, OutputConfig, iso8601 from .grid import GridInputType, GridDict, ShadingMethod, domain_mean, number_of_columns, validate_grid from .restart import validate_all_variables_present, validate_variable_shape @@ -41,12 +41,19 @@ def create_restart_file_name(time: datetime) -> str: return restart_file_prefix + iso8601(time) + ".nc" -def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initialize_from_restart_file: bool): +def init_config( + time_config: TimeConfig, + grid_config, + output_config: OutputConfig, + restartdir: Path, + initialize_from_restart_file: bool, +): """ Set model parameters, specify grid parameters, I/O, and physics settings. @param[in] time_config Time configuration object. @param[in] grid_config Grid configuration object. + @param[in] output_config Output file configuration object. @param[in] restartdir Path to the directory containing restart files (if initializing from restart file). @param[in] initialize_from_restart_file Boolean flag indicating whether to initialize from a restart file. @@ -84,7 +91,7 @@ def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initiali io = {} io["homedir"] = Path(os.getcwd()) # Home directory - io["outdir"] = io["homedir"] / "Output" # Output directory + io["outfile"] = output_config.output_file # Output file write_restart_file: bool @@ -106,10 +113,13 @@ def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initiali write_restart_file = False io["writebootfile"] = write_restart_file - io["freqout"] = 8 # OUTPUT: frequency of storing output (every n-th time-step) + io["freqout"] = output_config.output_frequency // time_config.time_step # OUTPUT: store every n-th time-step + assert ( + output_config.output_frequency % time_config.time_step == timedelta() + ), "Output frequency must be a multiple of the time step." # Ensure output and reboot directories exist - os.makedirs(io["outdir"], exist_ok=True) + os.makedirs(io["outfile"].parent, exist_ok=True) # Return the initialized parameters return grid, io, phys, column_config diff --git a/src/ebfm/core/LOOP_write_to_file.py b/src/ebfm/core/LOOP_write_to_file.py index 1f4e2f94..1535a02a 100644 --- a/src/ebfm/core/LOOP_write_to_file.py +++ b/src/ebfm/core/LOOP_write_to_file.py @@ -3,6 +3,7 @@ # SPDX-License-Identifier: BSD-3-Clause import os +from pathlib import Path import numpy as np from netCDF4 import Dataset, date2num @@ -254,12 +255,10 @@ def main(OUTFILE, io, OUT, grid, t, time, column, gather: ColumnGather): # Initialize NetCDF file at the first time step if is_first_time_step(t) and gather.is_root: - if not os.path.exists(io["outdir"]): - os.makedirs(io["outdir"]) + os.makedirs(Path(io["outfile"]).parent, exist_ok=True) # Create NetCDF file - nc_filepath = os.path.join(io["outdir"], "model_output.nc") - io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4") + io["nc_file"] = Dataset(io["outfile"], "w", format="NETCDF4") io["nc_file"].createDimension("time", None) # Unlimited time dimension if io["output_grid"]["is_unstructured"]: diff --git a/src/ebfm/core/cli.py b/src/ebfm/core/cli.py index 284671f7..cdbea58a 100644 --- a/src/ebfm/core/cli.py +++ b/src/ebfm/core/cli.py @@ -14,7 +14,7 @@ import argparse import os from argparse import ArgumentParser, Namespace -from datetime import datetime +from datetime import datetime, timedelta from enum import Enum from pathlib import Path @@ -61,6 +61,8 @@ class CliDefaults(Enum): TIME_STEP_SIZE = "PT3H" # = 3 hours LOG_LEVEL_CONSOLE = "INFO" COMPONENT_NAME = "ebfm" + OUTPUT_FILE = Path("Output") / "model_output.nc" + OUTPUT_FREQUENCY = "P1D" # = daily @classmethod def default_time_step_size_in_hours(cls) -> float: @@ -286,6 +288,36 @@ def parse_cli_args(args: list[str] | None = None) -> Namespace: help="Path to the folder containing NetCDF forcing files with meteorological data.", ) + output_group = parser.add_argument_group("output file") + + output_group.add_argument( + "--output-file", + type=Path, + default=CliDefaults.OUTPUT_FILE.value, + help="Path of the NetCDF output file. Relative paths are relative to the working directory. Missing " + "directories are created and an existing file is overwritten.", + ) + + def parse_output_frequency(value: str) -> timedelta: + try: + interval = isodate.parse_duration(value) + except isodate.ISO8601Error: + raise argparse.ArgumentTypeError(f"must be an ISO 8601 duration such as P1D, got '{value}'") + if not isinstance(interval, timedelta): + raise argparse.ArgumentTypeError(f"must not contain years or months, which vary in length, got '{value}'") + if interval <= timedelta(0): + raise argparse.ArgumentTypeError(f"must be positive, got '{value}'") + return interval + + output_group.add_argument( + "--output-frequency", + type=parse_output_frequency, + default=CliDefaults.OUTPUT_FREQUENCY.value, + help="Interval between two output records in ISO 8601 format, e.g. P1D for daily or PT6H for six-hourly " + "output. Must be a multiple of --time-step. Each record holds the mean, sum, or a sample over this " + "interval, depending on the variable.", + ) + time_group = parser.add_argument_group("time configuration") time_group.add_argument( @@ -613,6 +645,19 @@ def validate_forcing_parallel_compat(args: Namespace, comm) -> None: ) +def validate_output_frequency_time_step_compat(args: Namespace, time_config) -> None: + """Check that --output-frequency is a multiple of the time step. + + @param[in] args parsed CLI arguments + @param[in] time_config resolved TimeConfig + """ + if args.output_frequency % time_config.time_step: + _cli_error( + f"--output-frequency {isodate.duration_isoformat(args.output_frequency)} is not a multiple of the time " + f"step {isodate.duration_isoformat(time_config.time_step)} (--time-step)." + ) + + def validate_restart_parallel_compat(args: Namespace, comm) -> None: """Check that no restart files are read or written in a run on several MPI ranks. diff --git a/src/ebfm/core/config/__init__.py b/src/ebfm/core/config/__init__.py index 13c551ec..da816234 100644 --- a/src/ebfm/core/config/__init__.py +++ b/src/ebfm/core/config/__init__.py @@ -7,3 +7,4 @@ from ebfm.core.config.grid import GridConfig, GridInputType # noqa: F401 from ebfm.core.config.forcing import ForcingConfig, ForcingType # noqa: F401 from ebfm.core.config.column import ColumnDiscretizationConfig # noqa: F401 +from ebfm.core.config.output import OutputConfig # noqa: F401 diff --git a/src/ebfm/core/config/output.py b/src/ebfm/core/config/output.py new file mode 100644 index 00000000..4108a953 --- /dev/null +++ b/src/ebfm/core/config/output.py @@ -0,0 +1,26 @@ +# SPDX-FileCopyrightText: 2026 EBFM Authors +# +# SPDX-License-Identifier: BSD-3-Clause +"""Output configuration for EBFM. + +This file provides the configuration of the NetCDF output file for EBFM. +""" + +from argparse import Namespace +from datetime import timedelta +from pathlib import Path + + +class OutputConfig: + """Configuration of the NetCDF output file.""" + + def __init__(self, args: Namespace): + """ + Initialize output configuration from command line arguments. + + @param[in] args command line arguments + """ + # Relative paths are relative to the working directory + self.output_file: Path = Path.cwd() / args.output_file + # Interval between two output records + self.output_frequency: timedelta = args.output_frequency diff --git a/src/ebfm/main.py b/src/ebfm/main.py index e2ab2941..878769ee 100644 --- a/src/ebfm/main.py +++ b/src/ebfm/main.py @@ -18,7 +18,7 @@ ) from ebfm.core import LOOP_write_to_file, FINAL_create_restart_file from ebfm.core.comm.gather import ColumnGather -from ebfm.core.config import CouplingConfig, GridConfig, TimeConfig, ForcingConfig +from ebfm.core.config import CouplingConfig, GridConfig, TimeConfig, ForcingConfig, OutputConfig from ebfm.core.logger import Logger, setup_logging, log_levels_map, getLogger from ebfm.core.cli import ( extract_active_coupling_features, @@ -30,6 +30,7 @@ validate_forcing_parallel_compat, validate_restart_parallel_compat, validate_shading_parallel_compat, + validate_output_frequency_time_step_compat, ) import ebfm.coupling @@ -222,11 +223,16 @@ def _main_impl(): # Ensure shading routine is only used in runs on a single rank. validate_shading_parallel_compat(grid_config, ebfm_comm) + # Ensure every output record covers a whole number of time steps. + validate_output_frequency_time_step_compat(args, time_config) + logger.debug("Successfully completed consistency checks.") # Model setup & initialization time = time_config.to_dict() - grid, io, phys, column = INIT.init_config(time_config, grid_config, args.restart_dir, args.restart_init) + grid, io, phys, column = INIT.init_config( + time_config, grid_config, OutputConfig(args), args.restart_dir, args.restart_init + ) set_random_seed(logger, args.random_seed, ebfm_comm) diff --git a/tests/core/test_cli.py b/tests/core/test_cli.py index bb3e41f3..c4d15817 100644 --- a/tests/core/test_cli.py +++ b/tests/core/test_cli.py @@ -8,6 +8,7 @@ import io import os from contextlib import redirect_stderr, redirect_stdout +from datetime import timedelta from pathlib import Path from unittest.mock import MagicMock, patch @@ -18,6 +19,7 @@ parse_cli_args, validate_forcing_parallel_compat, validate_grid_parallel_compat, + validate_output_frequency_time_step_compat, validate_restart_parallel_compat, validate_shading_parallel_compat, validate_shading_coupling_compat, @@ -194,6 +196,51 @@ def test_netcdf_dem_mesh_default_none(self): def test_random_seed_default_none(self): self.assertIsNone(self.args.random_seed) + def test_output_defaults(self): + self.assertEqual(self.args.output_file, Path("Output") / "model_output.nc") + self.assertEqual(self.args.output_frequency, timedelta(days=1)) + + +class TestOutputOptions(unittest.TestCase): + """--output-file / --output-frequency parsing and checks.""" + + def test_custom_values_accepted(self): + args = _parse("--output-file", "/scratch/run1.nc", "--output-frequency", "PT6H") + self.assertEqual(args.output_file, Path("/scratch/run1.nc")) + self.assertEqual(args.output_frequency, timedelta(hours=6)) + + def test_frequency_other_than_positive_fixed_duration_rejected(self): + cases = { + "integer": "8", + "float": "1.5", + "word": "daily", + "zero": "PT0S", + "negative": "-P1D", + "months vary in length": "P1M", + "years vary in length": "P1Y", + } + for case, value in cases.items(): + with self.subTest(case, value=value): + with self.assertRaises(SystemExit), redirect_stderr(io.StringIO()): + _parse("--output-frequency", value) + + def test_frequency_must_be_multiple_of_time_step(self): + cases = { + "equal to time step": ("PT3H", "PT3H", True), + "multiple of time step": ("P1D", "PT3H", True), + "not a multiple": ("PT4H", "PT3H", False), + "shorter than time step": ("PT1H", "PT3H", False), + } + for case, (frequency, time_step, accepted) in cases.items(): + with self.subTest(case, frequency=frequency, time_step=time_step): + args = _parse("--output-frequency", frequency, "--time-step", time_step) + time_config = TimeConfig(args) + if accepted: + validate_output_frequency_time_step_compat(args, time_config) + else: + with self.assertRaises(SystemExit), redirect_stderr(io.StringIO()): + validate_output_frequency_time_step_compat(args, time_config) + class TestGpuOptions(unittest.TestCase): """--with-gpu / --gpu-vendor parsing and compatibility checks.""" diff --git a/tests/core/test_init.py b/tests/core/test_init.py index 68a954cf..74ed76ed 100644 --- a/tests/core/test_init.py +++ b/tests/core/test_init.py @@ -20,7 +20,7 @@ from ebfm.core.FINAL_create_restart_file import main as write_restart_file from ebfm.core.cli import parse_cli_args from ebfm.core.comm import defaultComm -from ebfm.core.config import ColumnDiscretizationConfig +from ebfm.core.config import ColumnDiscretizationConfig, OutputConfig from ebfm.core.config.grid import GridConfig from ebfm.core.config.time import TimeConfig from ebfm.core.restart import PER_COLUMN_VARIABLES, PER_LAYER_VARIABLES, RESTART_VARIABLES @@ -336,6 +336,56 @@ def test_the_two_groups_do_not_overlap(self): self.assertEqual(len(RESTART_VARIABLES), len(set(RESTART_VARIABLES))) +class TestOutputConfig(unittest.TestCase): + """The output options reach the io dict that LOOP_write_to_file reads.""" + + def setUp(self): + self.tmpdir = tempfile.TemporaryDirectory() + + def tearDown(self): + self.tmpdir.cleanup() + + def _init_io(self, *extra_args): + """Run init_config with the output file in the temporary directory and return its io dict.""" + outfile = Path(self.tmpdir.name) / "runs" / "a" / "run1.nc" + args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--output-file", str(outfile), *extra_args]) + _, io, _, _ = init_config(TimeConfig(args), GridConfig(args), OutputConfig(args), None, False) + return io + + def test_default_output_file_is_model_output_nc_in_output_of_working_directory(self): + output_config = OutputConfig(parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH)])) + + self.assertEqual(output_config.output_file, Path.cwd() / "Output" / "model_output.nc") + + def test_relative_output_file_is_relative_to_working_directory(self): + output_config = OutputConfig(parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--output-file", "a.nc"])) + + self.assertEqual(output_config.output_file, Path.cwd() / "a.nc") + + def test_output_file_ends_up_in_io_and_its_directory_is_created(self): + io = self._init_io() + + self.assertEqual(io["outfile"], Path(self.tmpdir.name) / "runs" / "a" / "run1.nc") + self.assertTrue(io["outfile"].parent.is_dir()) + + def test_default_output_frequency_is_daily(self): + io = self._init_io() + + self.assertEqual(io["freqout"], 8) # time steps of the default PT3H per day + + def test_output_frequency_is_converted_to_time_steps(self): + cases = { + "daily, hourly steps": ("P1D", "PT1H", 24), + "six-hourly, three-hourly steps": ("PT6H", "PT3H", 2), + "every time step": ("PT3H", "PT3H", 1), + "weekly, daily steps": ("P1W", "P1D", 7), + } + for case, (frequency, time_step, expected) in cases.items(): + with self.subTest(case): + io = self._init_io("--output-frequency", frequency, "--time-step", time_step) + self.assertEqual(io["freqout"], expected) + + class TestMatlabShadingLookupTable(unittest.TestCase): """The shading look-up table must only be pre-computed when shading is enabled. @@ -354,7 +404,7 @@ def tearDown(self): def _init_grid(self, shading_flag): args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), shading_flag]) grid_config = GridConfig(args) - grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False) + grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False) return init_grid(grid, io, grid_config, defaultComm) def test_shading_disabled_skips_the_lookup_table(self): @@ -392,7 +442,7 @@ def tearDown(self): def _init_grid(self, num_ranks, rank): args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--no-shading"]) grid_config = GridConfig(args, num_ranks, rank) - grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False) + grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False) return init_grid(grid, io, grid_config, defaultComm) def test_blocks_of_all_ranks_make_up_the_serial_grid(self): @@ -449,7 +499,7 @@ def _init_grid(self, extra_args=()): ["--elmer-mesh", self.tmpdir.name, "--elmer-mesh-crs-epsg", str(self._CRS_EPSG), *extra_args] ) grid_config = GridConfig(args) - grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False) + grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False) with mock.patch.object(INIT, "read_elmer_mesh", return_value=self.mesh): return init_grid(grid, io, grid_config, defaultComm) diff --git a/tests/core/test_loop_write_to_file.py b/tests/core/test_loop_write_to_file.py index 87499ba7..1e6b6116 100644 --- a/tests/core/test_loop_write_to_file.py +++ b/tests/core/test_loop_write_to_file.py @@ -133,12 +133,13 @@ def _write_single_step(outdir, grid, extra_out=None): out = _fake_out(GPSUM, NL) if extra_out: out.update(extra_out) - io = {"outdir": str(outdir), "freqout": 1} + outfile_path = Path(outdir) / "model_output.nc" + io = {"outfile": outfile_path, "freqout": 1} time = {"TCUR": datetime(2020, 1, 1), "tn": 1} _, outfile = write_output({}, io, out, grid, 0, time, column, ColumnGather(defaultComm, grid["column_id"])) - return Path(outdir) / "model_output.nc", outfile, out + return outfile_path, outfile, out class TestNetCDFOutputSmoke(unittest.TestCase):