From 1ed83707770930bb17733c6e371656d977dd9457 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Mon, 20 Oct 2025 15:37:17 +0200 Subject: [PATCH 01/25] read unstructured elmer-xios grid --- src/ebfm/INIT.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index d5f8686f..02fcd05a 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -19,6 +19,8 @@ import logging logger = logging.getLogger(__name__) +import sys + def init_config(args: Namespace): """ Set model parameters, specify grid parameters, model time period, I/O, and physics settings. @@ -182,6 +184,8 @@ def init_grid(grid, io, args: Namespace): grid['input_type'] = GridInputType.CUSTOM elif args.elmer_mesh: grid['input_type'] = GridInputType.ELMER + elif args.elmer_xios_mesh: + grid['input_type'] = GridInputType.ELMER_XIOS else: logger.error(f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. Please please refer to the documentation for correct configuration.") raise Exception("Invalid grid configuration.") @@ -319,10 +323,28 @@ def init_grid(grid, io, args: Namespace): grid['slope_gamma'][(grid['slope_x'] > 0) & (grid['slope_y'] == 0)] = np.pi / 2 grid['slope_gamma'][(grid['slope_x'] < 0) & (grid['slope_y'] == 0)] = -np.pi / 2 grid['slope_gamma'] = -grid['slope_gamma'] + elif grid['input_type'] is GridInputType.ELMER_XIOS: + grid = read_elmer_xios_grid(grid=grid, gridfile=args.elmer_xios_mesh) + grid['gpsum'] = grid['z'].shape[0] + grid['mask'] = (grid['h'] > 1.0) * 1.0 + grid['x'] = np.zeros_like(grid['z']) + grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! + grid["slope_gamma"] = np.zeros_like(grid["x"]) # test values! + return grid + +def read_elmer_xios_grid(grid, gridfile: Path): + import netCDF4 as nc + print(gridfile) + with nc.Dataset(gridfile) as file: + grid['lat'] = file['mesh2D_node_x'][:] + grid['lon'] = file['mesh2D_node_y'][:] + grid['z'] = np.squeeze(file['zs'][:].data) + grid['h'] = np.squeeze(file['h'][:].data) return grid + def read_MATLAB_grid(gridfile: Path): """ Provides grid information by reading from a .mat file or allowing user input. From 1e6dff88e2393fd2bb1e28e4ae1c3f43b5f2a6aa Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Mon, 20 Oct 2025 15:37:49 +0200 Subject: [PATCH 02/25] add elmer_xios type --- src/ebfm/grid.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index aa20e453..2e1bdc84 100644 --- a/src/ebfm/grid.py +++ b/src/ebfm/grid.py @@ -4,3 +4,4 @@ class GridInputType(Enum): MATLAB = 'matlab' CUSTOM = 'custom' ELMER = 'elmer' + ELMER_XIOS = 'elmer_xios' From 7890ea85228b36def8a5b5d3409e47bd2be7fadc Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Mon, 20 Oct 2025 15:38:54 +0200 Subject: [PATCH 03/25] add elmer-xios-mesh option --- src/dummy.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/dummy.py b/src/dummy.py index d316589e..7b890a23 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -15,8 +15,8 @@ log_levels = { 'file':logging.DEBUG, # log level for logging to file - 0: logging.INFO, # log level for rank 0 - # 1: logging.DEBUG, # to log other ranks to console define log level here + # 0: logging.INFO, # log level for rank 0 + 0: logging.DEBUG, # to log other ranks to console define log level here } setup_logging(log_levels=log_levels) @@ -45,6 +45,10 @@ def main(): type=Path, help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required.") + parser.add_argument("--elmer-xios-mesh", + type=Path, + help="Path to the Elmer-xios mesh file. One of the following options is required: --elmer-mesh, --matlab-mesh, --elmer-xios-mesh.") + parser.add_argument("--matlab-mesh", type=Path, help="Path to the MATLAB mesh file. Either --elmer-mesh or --matlab-mesh is required.") From 31384f8983eb6cf83c142ef0d42c435675730725 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Mon, 20 Oct 2025 16:34:01 +0200 Subject: [PATCH 04/25] allow output writing for unstructured grids --- src/dummy.py | 11 ++- src/ebfm/LOOP_write_to_file.py | 133 +++++++++++++++++++++++---------- 2 files changed, 100 insertions(+), 44 deletions(-) diff --git a/src/dummy.py b/src/dummy.py index 7b890a23..46686ed3 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -159,9 +159,14 @@ def main(): # Write output to files (only in uncoupled run and for unpartitioned grid) if not grid['is_partitioned'] and not coupler.has_coupling: - assert (grid['input_type'] is GridInputType.MATLAB), \ - "Output writing currently only implemented for MATLAB input grids." - io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time2, C) + if grid['input_type'] is GridInputType.ELMER_XIOS: + gridtype = 'unstructured' + elif grid['input_type'] is GridInputType.MATLAB: + gridtype = 'structured' + # assert (grid['input_type'] is GridInputType.MATLAB) or GridInputType.ELMER_XIOS, \ + # "Output writing currently only implemented for MATLAB input grids." + io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, + time2, C, gridtype) pass elif grid['is_partitioned'] or coupler.has_coupling: logger.warning('Skipping writing output to file for coupled or partitioned runs.') diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index d9fcce0a..a6dfad04 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -2,9 +2,10 @@ import os import numpy as np from netCDF4 import Dataset, date2num +import sys -def main(OUTFILE, io, OUT, grid, t, time, C): +def main(OUTFILE, io, OUT, grid, t, time, C, gridtype): # Specify variables to be written if t == 1: OUTFILE["varsout"] = [ @@ -117,7 +118,7 @@ def save_binary_files(): return True - def save_netCDF_file(): + def save_netCDF_file(gridtype): """ Save model output stored in OUTFILE to a NetCDF file. Converts 1D data to 2D grids. @@ -136,6 +137,7 @@ def save_netCDF_file(): time_units = "days since 1970-01-01 00:00:00" time_calendar = "gregorian" + # Initialize NetCDF file at the first time step if t == 1: if not os.path.exists(io["outdir"]): @@ -144,12 +146,17 @@ def save_netCDF_file(): # Create NetCDF file nc_filepath = os.path.join(io["outdir"], "model_output.nc") io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4") - - # Define dimensions - io["nc_file"].createDimension("time", None) # Unlimited time dimension - io["nc_file"].createDimension("y", grid["x_2D"].shape[0]) # 2D grid rows - io["nc_file"].createDimension("x", grid["x_2D"].shape[1]) # 2D grid columns - io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables + if gridtype == 'unstructured': + # Define dimensions + io["nc_file"].createDimension("time", None) # Unlimited time dimension + io["nc_file"].createDimension("y", grid["lat"].shape[0]) # 2D grid rows + io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables + elif gridtype == 'structured': + # Define dimensions + io["nc_file"].createDimension("time", None) # Unlimited time dimension + io["nc_file"].createDimension("y", grid["x_2D"].shape[0]) # 2D grid rows + io["nc_file"].createDimension("x", grid["x_2D"].shape[1]) # 2D grid columns + io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables # Define standard output variables for entry in OUTFILE["varsout"]: @@ -159,27 +166,51 @@ def save_netCDF_file(): # Check if variable is a `sub` variable if varname.startswith("sub"): - # Define variable as 4D: (time, y, x, nl) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y", "x", "nl"), - zlib=True, - complevel=4, - fill_value=-9999.0, # Fill missing values - chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]), - ) + if gridtype == "unstructured": + # Define variable as 3D: (time, y, nl) + nc_var = io["nc_file"].createVariable( + varname, + np.float32, + ("time", "y","nl"), + zlib=True, + complevel=4, + fill_value=-9999.0, # Fill missing values + chunksizes=(1, grid["lat"].shape[0], grid["nl"]), + ) + elif gridtype == "structured": + # Define variable as 4D: (time, y, x, nl) + nc_var = io["nc_file"].createVariable( + varname, + np.float32, + ("time", "y", "x", "nl"), + zlib=True, + complevel=4, + fill_value=-9999.0, # Fill missing values + chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]), + ) else: - # Define variable as 3D: (time, y, x) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y", "x"), - zlib=True, - complevel=4, - fill_value=-9999.0, - chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1]), - ) + if gridtype == "unstructured": + # Define variable as 2D: (time, y) + nc_var = io["nc_file"].createVariable( + varname, + np.float32, + ("time", "y"), + zlib=True, + complevel=4, + fill_value=-9999.0, + chunksizes=(1, grid["lat"].shape[0]), + ) + elif gridtype == "structured": + # Define variable as 3D: (time, y, x) + nc_var = io["nc_file"].createVariable( + varname, + np.float32, + ("time", "y", "x"), + zlib=True, + complevel=4, + fill_value=-9999.0, + chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1]), + ) # Assign metadata nc_var.units = var_units @@ -208,18 +239,38 @@ def save_netCDF_file(): varname = entry[0] var_1D = OUTFILE["TEMP"][varname] # 1D data - # Handle `sub` variables (4D: time, y, x, nl) - if varname.startswith("sub"): - var_3D = np.full((grid["x_2D"].size, grid["nl"]), -9999.0) - var_3D[grid["ind"], :] = var_1D - var_4D = var_3D.reshape(-1, grid["nl"]).reshape( - *grid["x_2D"].shape, grid["nl"] - ) - io["nc_file"][varname][time_index, :, :, :] = var_4D - else: - var_2D = np.full(grid["x_2D"].shape, -9999.0) - var_2D.flat[grid["ind"]] = var_1D - io["nc_file"][varname][time_index, :, :] = var_2D + if gridtype == 'unstructured': + # Handle `sub` variables (3D: time, y, nl) + if varname.startswith("sub"): + print("Doing 3D variables") + print(varname) + print(var_1D.shape) + # sys.exit() + # var_3D = np.full((grid["lat"].size, grid["nl"]), -9999.0) + # var_3D[grid["ind"], :] = var_1D + # var_4D = var_3D.reshape(-1, grid["nl"]).reshape( + # *grid["x_2D"].shape, grid["nl"] + # ) + io["nc_file"][varname][time_index, :, :] = var_1D + else: + print(varname) + print(var_1D.shape) + # var_2D = np.full(grid["x_2D"].shape, -9999.0) + # var_2D.flat[grid["ind"]] = var_1D + io["nc_file"][varname][time_index, :] = var_1D + elif gridtype == 'structured': + # Handle `sub` variables (4D: time, y, x, nl) + if varname.startswith("sub"): + var_3D = np.full((grid["x_2D"].size, grid["nl"]), -9999.0) + var_3D[grid["ind"], :] = var_1D + var_4D = var_3D.reshape(-1, grid["nl"]).reshape( + *grid["x_2D"].shape, grid["nl"] + ) + io["nc_file"][varname][time_index, :, :, :] = var_4D + else: + var_2D = np.full(grid["x_2D"].shape, -9999.0) + var_2D.flat[grid["ind"]] = var_1D + io["nc_file"][varname][time_index, :, :] = var_2D # Close the NetCDF file at the final time step if t == time["tn"]: @@ -231,7 +282,7 @@ def save_netCDF_file(): if io["output_type"] == 1: save_binary_files() elif io["output_type"] == 2: - save_netCDF_file() + save_netCDF_file(gridtype=gridtype) else: print("Invalid output type. Please choose 1 for binary files or 2 for NetCDF.") From cbb2c5c71a9d55df26387631be755e1858e503d6 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Thu, 23 Oct 2025 14:54:51 +0200 Subject: [PATCH 05/25] restructure code for reading surface geometry from elmer-xios output --- src/dummy.py | 24 ++++++++++++------------ src/ebfm/INIT.py | 21 ++++++++++++--------- src/ebfm/grid.py | 1 - src/reader.py | 21 ++++++++++++++++++++- 4 files changed, 44 insertions(+), 23 deletions(-) diff --git a/src/dummy.py b/src/dummy.py index 46686ed3..2651a8ac 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -45,10 +45,6 @@ def main(): type=Path, help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required.") - parser.add_argument("--elmer-xios-mesh", - type=Path, - help="Path to the Elmer-xios mesh file. One of the following options is required: --elmer-mesh, --matlab-mesh, --elmer-xios-mesh.") - parser.add_argument("--matlab-mesh", type=Path, help="Path to the MATLAB mesh file. Either --elmer-mesh or --matlab-mesh is required.") @@ -57,6 +53,10 @@ def main(): type=Path, help="Path to the NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.") + parser.add_argument("--netcdf-mesh-unstructured", + type=Path, + help="Path to the unstructured NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.") + parser.add_argument("--is-partitioned-elmer-mesh", action="store_true", help="Indicate if the provided Elmer mesh is partitioned for parallel runs.") @@ -159,15 +159,15 @@ def main(): # Write output to files (only in uncoupled run and for unpartitioned grid) if not grid['is_partitioned'] and not coupler.has_coupling: - if grid['input_type'] is GridInputType.ELMER_XIOS: - gridtype = 'unstructured' - elif grid['input_type'] is GridInputType.MATLAB: + if grid['input_type'] is GridInputType.MATLAB: gridtype = 'structured' - # assert (grid['input_type'] is GridInputType.MATLAB) or GridInputType.ELMER_XIOS, \ - # "Output writing currently only implemented for MATLAB input grids." - io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, - time2, C, gridtype) - pass + # if grid['input_type'] is GridInputType.MATLAB or grid['input_type'] is GridInputType.ELMER_XIOS: + # if grid['input_type'] is GridInputType.MATLAB: + # assert (grid['input_type'] is GridInputType.MATLAB), \ + # "Output writing currently only implemented for MATLAB and ELMER_XIOS input grids." + io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time2, C, gridtype=gridtype) + else: + logger.warning('Skipping writing output to file for Elmer input grids.') elif grid['is_partitioned'] or coupler.has_coupling: logger.warning('Skipping writing output to file for coupled or partitioned runs.') else: diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 02fcd05a..dcbe53aa 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -10,7 +10,7 @@ from netCDF4 import Dataset, num2date from pathlib import Path -from reader import read_elmer_mesh, read_dem +from reader import read_elmer_mesh, read_dem, read_dem_xios from elmer.mesh import Mesh from ebfm.grid import GridInputType @@ -180,12 +180,10 @@ def init_grid(grid, io, args: Namespace): if args.matlab_mesh: grid['input_type'] = GridInputType.MATLAB - elif args.netcdf_mesh and args.elmer_mesh: + elif (args.netcdf_mesh or args.netcdf_mesh_unstructured) and args.elmer_mesh: grid['input_type'] = GridInputType.CUSTOM elif args.elmer_mesh: grid['input_type'] = GridInputType.ELMER - elif args.elmer_xios_mesh: - grid['input_type'] = GridInputType.ELMER_XIOS else: logger.error(f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. Please please refer to the documentation for correct configuration.") raise Exception("Invalid grid configuration.") @@ -199,15 +197,18 @@ def init_grid(grid, io, args: Namespace): mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh) grid['x'], grid['y'] = mesh.x_vertices, mesh.y_vertices - grid['z'] = read_dem(args.netcdf_mesh, grid['x'], grid['y']) + if args.netcdf_mesh: + grid['z'] = read_dem(args.netcdf_mesh, grid['x'], grid['y']) + grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! + grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! + if args.netcdf_mesh_unstructured: + grid = read_dem_xios(args.netcdf_mesh_unstructured, grid) + grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier + grid["gpsum"] = np.sum(grid['mask'] == 1) # number of modelled grid cells grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! - grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! - grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! grid["slope_gamma"] = np.zeros_like(grid["x"]) # test values! - grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier - grid["gpsum"] = np.sum(grid['mask'] == 1) # number of modelled grid cells grid["mesh"] = mesh # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy @@ -339,6 +340,8 @@ def read_elmer_xios_grid(grid, gridfile: Path): with nc.Dataset(gridfile) as file: grid['lat'] = file['mesh2D_node_x'][:] grid['lon'] = file['mesh2D_node_y'][:] + grid['x'] = file['x'][:] + grid['y'] = file['y'][:] grid['z'] = np.squeeze(file['zs'][:].data) grid['h'] = np.squeeze(file['h'][:].data) return grid diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index 2e1bdc84..aa20e453 100644 --- a/src/ebfm/grid.py +++ b/src/ebfm/grid.py @@ -4,4 +4,3 @@ class GridInputType(Enum): MATLAB = 'matlab' CUSTOM = 'custom' ELMER = 'elmer' - ELMER_XIOS = 'elmer_xios' diff --git a/src/reader.py b/src/reader.py index 4d2b7032..494deafe 100644 --- a/src/reader.py +++ b/src/reader.py @@ -7,6 +7,7 @@ from elmer.mesh import TriangleMesh import elmer.parser +import sys def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: int = -1) -> TriangleMesh: @@ -65,6 +66,24 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: cell_ids=global_cell_ids ) +def read_dem_xios(dem_file: Path, grid: dict): + """Read digital elevation model (DEM) file. + Args: + dem_file (Path): Path to the DEM NetCDF file. + grid (dict): grid-related parameters + Returns: + grid (dict): dictionary containing grid-related parameters + """ + assert dem_file.is_file(), f"DEM file {dem_file} does not exist." + + import netCDF4 + nc = netCDF4.Dataset(dem_file) + assert (np.squeeze(nc['x'][:]).shape == grid['x'].shape), "Surface mesh and Elmer mesh do not have the same number of vertices" + grid['z'] = np.squeeze(nc['zs'][:]).data + grid['lat'] = np.squeeze(nc['mesh2D_node_x'][:]).data + grid['lon'] = np.squeeze(nc['mesh2D_node_y'][:]).data + return grid + def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]): """Read digital elevation model (DEM) file. Args: @@ -174,4 +193,4 @@ def fortran_style_sci(x, precision=15): y = mesh.y_vertices h = read_dem(args.dem, x, y) - write_dem_as_elmer(mesh, h, args.elmer_mesh / "MESH" / "mesh.nodes", allow_overwrite=True) \ No newline at end of file + write_dem_as_elmer(mesh, h, args.elmer_mesh / "MESH" / "mesh.nodes", allow_overwrite=True) From e6a2fc257cefc527317cc2dc7e377442c9943f08 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Thu, 23 Oct 2025 14:55:23 +0200 Subject: [PATCH 06/25] clean up print statements --- src/ebfm/LOOP_write_to_file.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index a6dfad04..a155b195 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -242,21 +242,8 @@ def save_netCDF_file(gridtype): if gridtype == 'unstructured': # Handle `sub` variables (3D: time, y, nl) if varname.startswith("sub"): - print("Doing 3D variables") - print(varname) - print(var_1D.shape) - # sys.exit() - # var_3D = np.full((grid["lat"].size, grid["nl"]), -9999.0) - # var_3D[grid["ind"], :] = var_1D - # var_4D = var_3D.reshape(-1, grid["nl"]).reshape( - # *grid["x_2D"].shape, grid["nl"] - # ) io["nc_file"][varname][time_index, :, :] = var_1D else: - print(varname) - print(var_1D.shape) - # var_2D = np.full(grid["x_2D"].shape, -9999.0) - # var_2D.flat[grid["ind"]] = var_1D io["nc_file"][varname][time_index, :] = var_1D elif gridtype == 'structured': # Handle `sub` variables (4D: time, y, x, nl) From 18331317b98d0e6b859bff85c0fbdb7d62d32a8c Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Thu, 6 Nov 2025 15:39:14 +0100 Subject: [PATCH 07/25] Fix formatting --- src/dummy.py | 20 +++++++++++--------- src/ebfm/INIT.py | 28 ++++++++++++++-------------- src/ebfm/LOOP_write_to_file.py | 7 +++---- src/reader.py | 16 ++++++++++------ 4 files changed, 38 insertions(+), 33 deletions(-) diff --git a/src/dummy.py b/src/dummy.py index 479a400d..4eb88998 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -77,7 +77,8 @@ def main(): parser.add_argument( "--netcdf-mesh-unstructured", type=Path, - help="Path to the unstructured NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.", + help="Path to the unstructured NetCDF mesh file. Optional if using --elmer-mesh." + " If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.", ) parser.add_argument( @@ -188,16 +189,17 @@ def main(): # IN['dhdy'] = data_from_elmer('dhdy') # Write output to files (only in uncoupled run and for unpartitioned grid) - if not grid['is_partitioned'] and not coupler.has_coupling: - if grid['input_type'] is GridInputType.MATLAB: - gridtype = 'structured' - # if grid['input_type'] is GridInputType.MATLAB or grid['input_type'] is GridInputType.ELMER_XIOS: - # if grid['input_type'] is GridInputType.MATLAB: - # assert (grid['input_type'] is GridInputType.MATLAB), \ - # "Output writing currently only implemented for MATLAB and ELMER_XIOS input grids." + if not grid["is_partitioned"] and not coupler.has_coupling: + if grid["input_type"] is GridInputType.MATLAB: + gridtype = "structured" + # if grid['input_type'] is GridInputType.MATLAB or grid['input_type'] is GridInputType.ELMER_XIOS: + # ... + # if grid['input_type'] is GridInputType.MATLAB: + # assert (grid['input_type'] is GridInputType.MATLAB), \ + # "Output writing currently only implemented for MATLAB and ELMER_XIOS input grids." io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time2, C, gridtype=gridtype) else: - logger.warning('Skipping writing output to file for Elmer input grids.') + logger.warning("Skipping writing output to file for Elmer input grids.") elif grid["is_partitioned"] or coupler.has_coupling: logger.warning("Skipping writing output to file for coupled or partitioned runs.") else: diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 98b60d49..2eb5bcdd 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -22,7 +22,6 @@ logger = logging.getLogger(__name__) -import sys def init_config(args: Namespace): """ @@ -208,13 +207,13 @@ def init_grid(grid, io, args: Namespace): grid["x"], grid["y"] = mesh.x_vertices, mesh.y_vertices if args.netcdf_mesh: - grid['z'] = read_dem(args.netcdf_mesh, grid['x'], grid['y']) + grid["z"] = read_dem(args.netcdf_mesh, grid["x"], grid["y"]) grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! if args.netcdf_mesh_unstructured: grid = read_dem_xios(args.netcdf_mesh_unstructured, grid) - grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier - grid["gpsum"] = np.sum(grid['mask'] == 1) # number of modelled grid cells + grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier + grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! @@ -342,11 +341,11 @@ def init_grid(grid, io, args: Namespace): grid["slope_gamma"][(grid["slope_x"] > 0) & (grid["slope_y"] == 0)] = np.pi / 2 grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2 grid["slope_gamma"] = -grid["slope_gamma"] - elif grid['input_type'] is GridInputType.ELMER_XIOS: + elif grid["input_type"] is GridInputType.ELMER_XIOS: grid = read_elmer_xios_grid(grid=grid, gridfile=args.elmer_xios_mesh) - grid['gpsum'] = grid['z'].shape[0] - grid['mask'] = (grid['h'] > 1.0) * 1.0 - grid['x'] = np.zeros_like(grid['z']) + grid["gpsum"] = grid["z"].shape[0] + grid["mask"] = (grid["h"] > 1.0) * 1.0 + grid["x"] = np.zeros_like(grid["z"]) grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! grid["slope_gamma"] = np.zeros_like(grid["x"]) # test values! return grid @@ -354,14 +353,15 @@ def init_grid(grid, io, args: Namespace): def read_elmer_xios_grid(grid, gridfile: Path): import netCDF4 as nc + print(gridfile) with nc.Dataset(gridfile) as file: - grid['lat'] = file['mesh2D_node_x'][:] - grid['lon'] = file['mesh2D_node_y'][:] - grid['x'] = file['x'][:] - grid['y'] = file['y'][:] - grid['z'] = np.squeeze(file['zs'][:].data) - grid['h'] = np.squeeze(file['h'][:].data) + grid["lat"] = file["mesh2D_node_x"][:] + grid["lon"] = file["mesh2D_node_y"][:] + grid["x"] = file["x"][:] + grid["y"] = file["y"][:] + grid["z"] = np.squeeze(file["zs"][:].data) + grid["h"] = np.squeeze(file["h"][:].data) return grid diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 68e6061d..3f2c4b1f 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -5,7 +5,6 @@ import os import numpy as np from netCDF4 import Dataset, date2num -import sys def main(OUTFILE, io, OUT, grid, t, time, C, gridtype): @@ -142,12 +141,12 @@ def save_netCDF_file(gridtype): # Create NetCDF file nc_filepath = os.path.join(io["outdir"], "model_output.nc") io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4") - if gridtype == 'unstructured': + if gridtype == "unstructured": # Define dimensions io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["lat"].shape[0]) # 2D grid rows io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables - elif gridtype == 'structured': + elif gridtype == "structured": # Define dimensions io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["x_2D"].shape[0]) # 2D grid rows @@ -167,7 +166,7 @@ def save_netCDF_file(gridtype): nc_var = io["nc_file"].createVariable( varname, np.float32, - ("time", "y","nl"), + ("time", "y", "nl"), zlib=True, complevel=4, fill_value=-9999.0, # Fill missing values diff --git a/src/reader.py b/src/reader.py index 86ac1c00..bdae6c42 100644 --- a/src/reader.py +++ b/src/reader.py @@ -11,7 +11,6 @@ from elmer.mesh import TriangleMesh import elmer.parser -import sys def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: int = -1) -> TriangleMesh: @@ -89,7 +88,8 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: cell_ids=global_cell_ids, ) -def read_dem_xios(dem_file: Path, grid: dict): + +def read_dem_xios(dem_file: Path, grid: dict): """Read digital elevation model (DEM) file. Args: dem_file (Path): Path to the DEM NetCDF file. @@ -100,13 +100,17 @@ def read_dem_xios(dem_file: Path, grid: dict): assert dem_file.is_file(), f"DEM file {dem_file} does not exist." import netCDF4 + nc = netCDF4.Dataset(dem_file) - assert (np.squeeze(nc['x'][:]).shape == grid['x'].shape), "Surface mesh and Elmer mesh do not have the same number of vertices" - grid['z'] = np.squeeze(nc['zs'][:]).data - grid['lat'] = np.squeeze(nc['mesh2D_node_x'][:]).data - grid['lon'] = np.squeeze(nc['mesh2D_node_y'][:]).data + assert ( + np.squeeze(nc["x"][:]).shape == grid["x"].shape + ), "Surface mesh and Elmer mesh do not have the same number of vertices" + grid["z"] = np.squeeze(nc["zs"][:]).data + grid["lat"] = np.squeeze(nc["mesh2D_node_x"][:]).data + grid["lon"] = np.squeeze(nc["mesh2D_node_y"][:]).data return grid + def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]): """Read digital elevation model (DEM) file. Args: From 3cb21e9344a3417d41ae58ae33058da2cce39cc7 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Tue, 11 Nov 2025 07:50:09 +0100 Subject: [PATCH 08/25] avoid using surface topography from elmer --- src/coupling.py | 2 ++ src/dummy.py | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/coupling.py b/src/coupling.py index 09aa8bbc..54dce22e 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -305,10 +305,12 @@ def exchange_elmer_ice(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra received_data = {} for field_name, field in self.target_fields.items(): if field_name in elmer_fields: + logger.debug(f"EBFM: {field.component_name}, {field.grid_name}, {field.name}") role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) assert role == yac.ExchangeType.TARGET, f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." # Only get the first element, since we only have collection size of 1 received_data[field_name] = self._get(field_name)[0] + logger.debug(received_data) return received_data diff --git a/src/dummy.py b/src/dummy.py index 2651a8ac..e1d84d62 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -151,8 +151,8 @@ def main(): logger.debug('Done.') logger.debug('Received the following data from Elmer/Ice:', data_from_elmer) - IN['h'] = data_from_elmer['h'] - grid['z'] = IN['h'][0].ravel() + # IN['h'] = data_from_elmer['h'] + # grid['z'] = IN['h'][0].ravel() # TODO add gradient field later # IN['dhdx'] = data_from_elmer('dhdx') # IN['dhdy'] = data_from_elmer('dhdy') From 4be61a84fe7a87563df5758b3221ef77328b10c7 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Tue, 11 Nov 2025 07:52:28 +0100 Subject: [PATCH 09/25] adjust runtime that it matches with elmer --- src/ebfm/INIT.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index dcbe53aa..5b211904 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -36,7 +36,7 @@ def init_config(args: Namespace): # --------------------------------------------------------------------- time = {} time['ts'] = datetime.strptime('1-Jan-1979 00:00', '%d-%b-%Y %H:%M') # Start date and time - time['te'] = datetime.strptime('2-Jan-1979 09:00', '%d-%b-%Y %H:%M') # End date and time + time['te'] = datetime.strptime('2-Jan-1979 00:00', '%d-%b-%Y %H:%M') # End date and time time['dt'] = 0.125 # Time step in days # Calculate the number of time steps From 29350ade58c1375d1f3f870bd3eedb6ef10b64b2 Mon Sep 17 00:00:00 2001 From: Clemens Schannwell Date: Tue, 11 Nov 2025 09:31:38 +0100 Subject: [PATCH 10/25] fix formatting --- src/dummy.py | 3 ++- src/ebfm/INIT.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/dummy.py b/src/dummy.py index e8ea95e0..17651b9e 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -191,7 +191,8 @@ def main(): if not grid["is_partitioned"] and not coupler.has_coupling: if grid["input_type"] is GridInputType.MATLAB: gridtype = "structured" - # if grid['input_type'] is GridInputType.MATLAB or grid['input_type'] is GridInputType.ELMER_XIOS: + # if grid['input_type'] is GridInputType.MATLAB or + # grid['input_type'] is GridInputType.ELMER_XIOS: # ... # if grid['input_type'] is GridInputType.MATLAB: # assert (grid['input_type'] is GridInputType.MATLAB), \ diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index e4a863bd..1b4e2d28 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -37,9 +37,9 @@ def init_config(args: Namespace): # Time parameters # --------------------------------------------------------------------- time = {} - time['ts'] = datetime.strptime('1-Jan-1979 00:00', '%d-%b-%Y %H:%M') # Start date and time - time['te'] = datetime.strptime('2-Jan-1979 00:00', '%d-%b-%Y %H:%M') # End date and time - time['dt'] = 0.125 # Time step in days + time["ts"] = datetime.strptime("1-Jan-1979 00:00", "%d-%b-%Y %H:%M") # Start date and time + time["te"] = datetime.strptime("2-Jan-1979 00:00", "%d-%b-%Y %H:%M") # End date and time + time["dt"] = 0.125 # Time step in days # Calculate the number of time steps time["tn"] = int(round((time["te"] - time["ts"]).total_seconds() / 86400 / time["dt"])) + 1 From dcea3f7115a1c4504cb741b7df939f2a06bbbd3e Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Fri, 14 Nov 2025 12:13:29 +0100 Subject: [PATCH 11/25] Rename (for safer merge). --- src/{dummy.py => main.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{dummy.py => main.py} (100%) diff --git a/src/dummy.py b/src/main.py similarity index 100% rename from src/dummy.py rename to src/main.py From 1428593c0ed27cf8100b8cfe2e69956d7c15ca24 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 25 Nov 2025 12:37:03 +0100 Subject: [PATCH 12/25] Add minimal EBFMCouplingConfig class. --- src/coupling.py | 15 ++++++--------- src/ebfm/config.py | 28 ++++++++++++++++++++++++++++ src/main.py | 16 +++++++++++----- 3 files changed, 45 insertions(+), 14 deletions(-) create mode 100644 src/ebfm/config.py diff --git a/src/coupling.py b/src/coupling.py index ded4df10..b6adb5f2 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -9,6 +9,7 @@ from typing import Dict from elmer.mesh import Mesh as Grid # for now use an alias +from ebfm.config import EBFMCouplingConfig # from ebfm.geometry import Grid # TODO: consider introducing a new data structure native to EBFM? @@ -372,26 +373,22 @@ def finalize(self): del self.interface -def init( - coupler_config: Path, ebfm_coupling_config: Path, couple_with_icon_atmo: bool, couple_with_elmer_ice: bool -) -> Coupler: +def init(coupler_config: Path, ebfm_coupling_config: EBFMCouplingConfig) -> Coupler: """Create interface to the coupler and register component @param[in] path to global Coupler configuration file - @param[in] ebfm_coupling_config path to local configuration of coupling for EBFM - @param[in] couple_with_icon_atmo whether to couple with ICON atmosphere - @param[in] couple_with_elmer_ice whether to couple with Elmer/Ice + @param[in] ebfm_coupling_config local configuration of coupling for EBFM @returns Coupler object """ # TODO: use ebfm_coupling_config for EBFM specific configuration? - component_name = "ebfm" # TODO: get from ebfm_coupling_config? + component_name = ebfm_coupling_config.component_name coupler = YACCoupler(component_name, coupler_config) - coupler.couple_to_icon_atmo = couple_with_icon_atmo - coupler.couple_to_elmer_ice = couple_with_elmer_ice + coupler.couple_to_icon_atmo = ebfm_coupling_config.couple_with_icon_atmo + coupler.couple_to_elmer_ice = ebfm_coupling_config.couple_with_elmer_ice # TODO: yac_cget_comp_comm(comp_id, elmer_comm) needed? diff --git a/src/ebfm/config.py b/src/ebfm/config.py new file mode 100644 index 00000000..9f278626 --- /dev/null +++ b/src/ebfm/config.py @@ -0,0 +1,28 @@ +# SPDX-FileCopyrightText: 2025 EBFM Authors +# +# SPDX-License-Identifier: BSD-3-Clause +"""Configuration for EBFM coupling options. + +This file exposes a `EBFMCouplingConfig` dataclass. +""" + +from dataclasses import dataclass + + +@dataclass +class EBFMCouplingConfig: + """Coupling configuration. + + Keep this class intentionally small; add fields later as needed. + """ + + component_name: str # Name of this component + couple_to_icon_atmo: bool # Whether to couple this component to ICON atmosphere + couple_to_elmer_ice: bool # Whether to couple this component to Elmer/Ice + + def defines_coupling(self) -> bool: + """Check if any coupling is defined in this configuration. + + @returns True if coupling to any component is enabled, False otherwise + """ + return self.couple_to_icon_atmo or self.couple_to_elmer_ice diff --git a/src/main.py b/src/main.py index 253506d9..63d5138c 100644 --- a/src/main.py +++ b/src/main.py @@ -16,6 +16,7 @@ ) from ebfm import LOOP_write_to_file, FINAL_create_restart_file from ebfm.grid import GridInputType +from ebfm.config import EBFMCouplingConfig from mpi4py import MPI from utils import setup_logging @@ -167,6 +168,13 @@ def main(): logger.info(f"Starting EBFM version {ebfm.get_version()}...") + # TODO consider introducing an ebfm_adapter_config.yaml to be parsed alternatively/additionally to command line args + ebfm_coupling_config = EBFMCouplingConfig( + component_name="ebfm", + couple_to_icon_atmo=args.couple_to_icon_atmo, + couple_to_elmer_ice=args.couple_to_elmer_ice, + ) + logger.info("Done parsing command line arguments.") logger.debug("Parsed the following command line arguments:") for arg, val in vars(args).items(): @@ -179,18 +187,16 @@ def main(): OUT, IN, OUTFILE = INIT.init_initial_conditions(C, grid, io, time2) - if args.couple_to_icon_atmo or args.couple_to_elmer_ice: + if ebfm_coupling_config.defines_coupling(): # TODO: introduce minimal stub implementation # TODO consider introducing an ebfm_adapter_config.yaml coupler = coupling.init( coupler_config=args.coupler_config, - ebfm_coupling_config=Path("dummies") / "EBFM" / "ebfm-config.yaml", - couple_with_icon_atmo=args.couple_to_icon_atmo, - couple_with_elmer_ice=args.couple_to_elmer_ice, + component_coupling_config=ebfm_coupling_config, ) coupling.setup(coupler, grid["mesh"], time2) else: - coupler = NoCoupler("ebfm") + coupler = NoCoupler(component_name=ebfm_coupling_config.component_name) # Time-loop logger.info("Entering time loop...") From e163957d119f5a71dd762a4db1722082d0d71103 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 25 Nov 2025 19:47:59 +0100 Subject: [PATCH 13/25] Introduce GridConfig and general refactoring. --- src/coupling.py | 19 ++++----- src/ebfm/INIT.py | 57 +++++++-------------------- src/ebfm/config.py | 97 ++++++++++++++++++++++++++++++++++++++++++---- src/main.py | 32 +++++++-------- 4 files changed, 128 insertions(+), 77 deletions(-) diff --git a/src/coupling.py b/src/coupling.py index b6adb5f2..fbbd4ce9 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -9,7 +9,7 @@ from typing import Dict from elmer.mesh import Mesh as Grid # for now use an alias -from ebfm.config import EBFMCouplingConfig +from ebfm.config import CouplingConfig # from ebfm.geometry import Grid # TODO: consider introducing a new data structure native to EBFM? @@ -321,6 +321,7 @@ def exchange_icon_atmo(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra logger.debug(f"Receiving field {field_name} from ICON atmosphere...") # Only get the first element, since we only have collection size of 1 received_data[field_name] = self._get(field_name)[0] + logger.debug(f"Received field '{field_name}' from ICON. Data: {received_data[field_name]}") else: logger.debug(f"Skipping field {field_name} as it is not part of ICON atmosphere exchange.") @@ -365,6 +366,7 @@ def exchange_elmer_ice(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra ), f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." # Only get the first element, since we only have collection size of 1 received_data[field_name] = self._get(field_name)[0] + logger.debug(f"Received field '{field_name}' from Elmer/Ice. Data: {received_data[field_name]}") return received_data @@ -373,22 +375,21 @@ def finalize(self): del self.interface -def init(coupler_config: Path, ebfm_coupling_config: EBFMCouplingConfig) -> Coupler: +def init(coupling_config: CouplingConfig) -> Coupler: """Create interface to the coupler and register component @param[in] path to global Coupler configuration file - @param[in] ebfm_coupling_config local configuration of coupling for EBFM @returns Coupler object """ - # TODO: use ebfm_coupling_config for EBFM specific configuration? + # TODO: use coupling_config for EBFM specific configuration? - component_name = ebfm_coupling_config.component_name - coupler = YACCoupler(component_name, coupler_config) + component_name = coupling_config.component_name + coupler = YACCoupler(component_name, coupling_config.coupler_config) - coupler.couple_to_icon_atmo = ebfm_coupling_config.couple_with_icon_atmo - coupler.couple_to_elmer_ice = ebfm_coupling_config.couple_with_elmer_ice + coupler.couple_to_icon_atmo = coupling_config.couple_to_icon_atmo + coupler.couple_to_elmer_ice = coupling_config.couple_to_elmer_ice # TODO: yac_cget_comp_comm(comp_id, elmer_comm) needed? @@ -403,7 +404,7 @@ def setup(coupler: Coupler, grid: Grid, time: Dict[str, float]): @param[in] grid Grid used by EBFM where coupling happens @param[in] time dictionary with time parameters, e.g. {'tn': 12, 'dt': 0.125} """ - grid_name = "ebfm_grid" # TODO: get from ebfm_coupling_config? + grid_name = "ebfm_grid" # TODO: get from coupling_config? coupler.add_grid(grid_name, grid) coupler.add_couples(time) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 782af163..76af2b3d 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -4,7 +4,6 @@ import os from datetime import datetime -from argparse import Namespace from typing import Any import numpy as np from numpy import ndarray, dtype @@ -16,6 +15,7 @@ from reader import read_elmer_mesh, read_dem from elmer.mesh import Mesh +from ebfm.config import GridConfig from ebfm.grid import GridInputType import logging @@ -23,7 +23,7 @@ logger = logging.getLogger(__name__) -def init_config(args: Namespace): +def init_config(): """ Set model parameters, specify grid parameters, model time period, I/O, and physics settings. Returns: @@ -160,53 +160,22 @@ def init_constants(): return C -def init_grid(grid, io, args: Namespace): - - # Check configuration - - if args.elmer_mesh and args.matlab_mesh: - logger.error("Please provide either --elmer-mesh or --matlab-mesh, not both.") - raise Exception("Invalid grid configuration.") - - if args.is_partitioned_elmer_mesh and not args.elmer_mesh: - logger.error("--is-partitioned-elmer-mesh requires --elmer-mesh.") - raise Exception("Invalid grid configuration.") - - grid["is_partitioned"] = args.is_partitioned_elmer_mesh - if grid["is_partitioned"]: - assert args.netcdf_mesh, ( - "--is-partitioned-elmer-mesh requires --netcdf-mesh. " - "(Without --netcdf-mesh should also work but is untested.)" - ) - logger.info("Using partitioned grid...") - else: - logger.info("Using non-partitioned grid...") - - if args.matlab_mesh: - grid["input_type"] = GridInputType.MATLAB - elif args.netcdf_mesh and args.elmer_mesh: - grid["input_type"] = GridInputType.CUSTOM - elif args.elmer_mesh: - grid["input_type"] = GridInputType.ELMER - else: - logger.error( - f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. " - "Please please refer to the documentation for correct configuration." - ) - raise Exception("Invalid grid configuration.") +def init_grid(grid, io, config: GridConfig): + grid["input_type"] = config.grid_type + grid["is_partitioned"] = config.is_partitioned if grid["input_type"] is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine - if grid["is_partitioned"]: + if config.is_partitioned: mesh: Mesh = read_elmer_mesh( - mesh_root=args.elmer_mesh, - is_partitioned=True, - partition_id=args.use_part, + mesh_root=config.mesh_file, + is_partitioned=config.is_partitioned, + partition_id=config.partition_id, ) else: - mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh) + mesh: Mesh = read_elmer_mesh(mesh_root=config.mesh_file) grid["x"], grid["y"] = mesh.x_vertices, mesh.y_vertices - grid["z"] = read_dem(args.netcdf_mesh, grid["x"], grid["y"]) + grid["z"] = read_dem(config.dem_file, grid["x"], grid["y"]) grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! @@ -219,7 +188,7 @@ def init_grid(grid, io, args: Namespace): # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy elif grid["input_type"] is GridInputType.ELMER: # Read grid and elevations from Elmer - mesh: Mesh = read_elmer_mesh(args.elmer_mesh) + mesh: Mesh = read_elmer_mesh(config.mesh_file) # assuming mesh/MESH/mesh.nodes contains DEM data in the z component # see mesh/README.md for the required preprocessing steps. @@ -246,7 +215,7 @@ def init_grid(grid, io, args: Namespace): # Read and process grid information # --------------------------------------------------------------------- # Read grid data - input_data = read_MATLAB_grid(args.matlab_mesh) + input_data = read_MATLAB_grid(config.mesh_file) grid["x_2D"] = input_data["x"][0][0] grid["y_2D"] = input_data["y"][0][0] grid["z_2D"] = input_data["z"][0][0] diff --git a/src/ebfm/config.py b/src/ebfm/config.py index 9f278626..18ece779 100644 --- a/src/ebfm/config.py +++ b/src/ebfm/config.py @@ -1,24 +1,48 @@ # SPDX-FileCopyrightText: 2025 EBFM Authors # # SPDX-License-Identifier: BSD-3-Clause -"""Configuration for EBFM coupling options. +"""Configuration for EBFM. -This file exposes a `EBFMCouplingConfig` dataclass. +This file exposes a some configuration dataclasses for EBFM components. """ -from dataclasses import dataclass +from argparse import Namespace +from pathlib import Path +from ebfm.grid import GridInputType -@dataclass -class EBFMCouplingConfig: - """Coupling configuration. +import logging - Keep this class intentionally small; add fields later as needed. +logger = logging.getLogger(__name__) + + +class CouplingConfig: + """ + Coupling configuration. """ component_name: str # Name of this component couple_to_icon_atmo: bool # Whether to couple this component to ICON atmosphere couple_to_elmer_ice: bool # Whether to couple this component to Elmer/Ice + coupler_config: Path # Path to the coupler configuration file + + def __init__(self, args: Namespace, component_name: str): + """ + Initialize coupling configuration from command line arguments. + + @param[in] args command line arguments + @param[in] component_name name of this component + """ + + self.component_name = component_name + self.couple_to_icon_atmo = args.couple_to_icon_atmo + self.couple_to_elmer_ice = args.couple_to_elmer_ice + + if self.defines_coupling() and not args.coupler_config: + logger.error("Coupling enabled but no coupler configuration file provided (--coupler-config).") + raise Exception("Missing coupler configuration file.") + + self.coupler_config = args.coupler_config def defines_coupling(self) -> bool: """Check if any coupling is defined in this configuration. @@ -26,3 +50,62 @@ def defines_coupling(self) -> bool: @returns True if coupling to any component is enabled, False otherwise """ return self.couple_to_icon_atmo or self.couple_to_elmer_ice + + +class GridConfig: + """ + Grid configuration. + """ + + grid_type: GridInputType # Name of the grid used in coupling + mesh_file: Path # Path to the grid file + dem_file: Path # Path to the DEM file (only relevant for CUSTOM grid type) + is_partitioned: bool # Whether the grid is partitioned + partition_id: int # Partition ID (only relevant if is_partitioned is True) + + def __init__(self, args: Namespace): + """ + Initialize grid configuration from command line arguments. + + @param[in] args command line arguments + """ + if not (args.elmer_mesh or args.matlab_mesh): + logger.error("Grid needed. Please provide either --elmer-mesh or --matlab-mesh.") + raise Exception("Missing grid.") + + if args.elmer_mesh and args.matlab_mesh: + logger.error("Please provide either --elmer-mesh or --matlab-mesh, not both.") + raise Exception("Invalid grid configuration.") + + if args.is_partitioned_elmer_mesh and not args.elmer_mesh: + logger.error("--is-partitioned-elmer-mesh requires --elmer-mesh.") + raise Exception("Invalid grid configuration.") + + self.is_partitioned = args.is_partitioned_elmer_mesh + if self.is_partitioned: + assert args.netcdf_mesh, ( + "--is-partitioned-elmer-mesh requires --netcdf-mesh. " + "(Without --netcdf-mesh should also work but is untested.)" + ) + logger.info("Using partitioned grid...") + self.partition_id = args.use_part + logger.info(f"{self.partition_id=}") + else: + logger.info("Using non-partitioned grid...") + + if args.matlab_mesh: + self.grid_type = GridInputType.MATLAB + self.mesh_file = args.matlab_mesh + elif args.netcdf_mesh and args.elmer_mesh: + self.grid_type = GridInputType.CUSTOM + self.mesh_file = args.elmer_mesh + self.dem_file = args.netcdf_mesh + elif args.elmer_mesh: + self.grid_type = GridInputType.ELMER + self.mesh_file = args.elmer_mesh + else: + logger.error( + f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. " + "Please refer to the documentation for correct configuration." + ) + raise Exception("Invalid grid configuration.") diff --git a/src/main.py b/src/main.py index 63d5138c..3a5deb0d 100644 --- a/src/main.py +++ b/src/main.py @@ -16,7 +16,7 @@ ) from ebfm import LOOP_write_to_file, FINAL_create_restart_file from ebfm.grid import GridInputType -from ebfm.config import EBFMCouplingConfig +from ebfm.config import CouplingConfig, GridConfig from mpi4py import MPI from utils import setup_logging @@ -159,7 +159,7 @@ def main(): {coupling_import_error} -Hint: If you are missing 'yac', please install YAC and the python binds as described under +Hint: If you are missing 'yac', please install YAC and the python bindings as described under https://dkrz-sw.gitlab-pages.dkrz.de/yac/d1/d9f/installing_yac.html" """ ) @@ -168,35 +168,33 @@ def main(): logger.info(f"Starting EBFM version {ebfm.get_version()}...") - # TODO consider introducing an ebfm_adapter_config.yaml to be parsed alternatively/additionally to command line args - ebfm_coupling_config = EBFMCouplingConfig( - component_name="ebfm", - couple_to_icon_atmo=args.couple_to_icon_atmo, - couple_to_elmer_ice=args.couple_to_elmer_ice, - ) - logger.info("Done parsing command line arguments.") logger.debug("Parsed the following command line arguments:") for arg, val in vars(args).items(): logger.debug(f" {arg}: {val}") + logger.debug("Reading configuration and checking for consistency.") + + # TODO consider introducing an ebfm_adapter_config.yaml to be parsed alternatively/additionally to command line args + coupling_config = CouplingConfig(args, component_name="ebfm") + grid_config = GridConfig(args) + + logger.debug("Successfully completed consistency checks.") + # Model setup & initialization - grid, time2, io, phys = INIT.init_config(args) + grid, time2, io, phys = INIT.init_config() C = INIT.init_constants() - grid = INIT.init_grid(grid, io, args) + grid = INIT.init_grid(grid, io, grid_config) OUT, IN, OUTFILE = INIT.init_initial_conditions(C, grid, io, time2) - if ebfm_coupling_config.defines_coupling(): + if coupling_config.defines_coupling(): # TODO: introduce minimal stub implementation # TODO consider introducing an ebfm_adapter_config.yaml - coupler = coupling.init( - coupler_config=args.coupler_config, - component_coupling_config=ebfm_coupling_config, - ) + coupler = coupling.init(coupling_config=coupling_config) coupling.setup(coupler, grid["mesh"], time2) else: - coupler = NoCoupler(component_name=ebfm_coupling_config.component_name) + coupler = NoCoupler(component_name=coupling_config.component_name) # Time-loop logger.info("Entering time loop...") From b6c1b0b697eb73c570e7d221251ff9f9292b7b74 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 25 Nov 2025 20:21:14 +0100 Subject: [PATCH 14/25] Simplify and check for shading. --- .gitignore | 2 ++ src/ebfm/INIT.py | 12 ++++++++---- src/ebfm/LOOP_EBM_insolation.py | 4 +--- src/main.py | 9 ++++++++- 4 files changed, 19 insertions(+), 8 deletions(-) diff --git a/.gitignore b/.gitignore index ba6a9644..ccae9d66 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ **__pycache__ **.egg-info build +Output/*.nc +Reboot/*.nc diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 76af2b3d..afdadcc3 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -161,10 +161,9 @@ def init_constants(): def init_grid(grid, io, config: GridConfig): - grid["input_type"] = config.grid_type grid["is_partitioned"] = config.is_partitioned - if grid["input_type"] is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine + if config.grid_type is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine if config.is_partitioned: mesh: Mesh = read_elmer_mesh( mesh_root=config.mesh_file, @@ -185,9 +184,10 @@ def init_grid(grid, io, config: GridConfig): grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells grid["mesh"] = mesh + grid["has_shading"] = False # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy - elif grid["input_type"] is GridInputType.ELMER: # Read grid and elevations from Elmer + elif config.grid_type is GridInputType.ELMER: # Read grid and elevations from Elmer mesh: Mesh = read_elmer_mesh(config.mesh_file) # assuming mesh/MESH/mesh.nodes contains DEM data in the z component @@ -210,7 +210,8 @@ def init_grid(grid, io, config: GridConfig): # TODO later add slope # grid["slope_x"], grid["slope_y"] = mesh.dzdy, mesh.dzdy grid["mesh"] = mesh - elif grid["input_type"] is GridInputType.MATLAB: # Read grid and elevations from example MATLAB file + grid["has_shading"] = False + elif config.grid_type is GridInputType.MATLAB: # Read grid and elevations from example MATLAB file # --------------------------------------------------------------------- # Read and process grid information # --------------------------------------------------------------------- @@ -222,6 +223,7 @@ def init_grid(grid, io, config: GridConfig): grid["mask_2D"] = input_data["mask"][0][0] # Determine domain extent + grid["has_shading"] = True grid["Lx"], grid["Ly"] = grid["x_2D"].shape # Flip grid E-W or N-S when needed @@ -307,6 +309,8 @@ def init_grid(grid, io, config: GridConfig): grid["slope_gamma"][(grid["slope_x"] > 0) & (grid["slope_y"] == 0)] = np.pi / 2 grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2 grid["slope_gamma"] = -grid["slope_gamma"] + else: + raise ValueError(f"Unsupported grid input type {config.grid_type} specified in configuration.") return grid diff --git a/src/ebfm/LOOP_EBM_insolation.py b/src/ebfm/LOOP_EBM_insolation.py index 94ba5bc0..27d3a6c4 100644 --- a/src/ebfm/LOOP_EBM_insolation.py +++ b/src/ebfm/LOOP_EBM_insolation.py @@ -5,8 +5,6 @@ import numpy as np import datetime -from ebfm.grid import GridInputType - def main(grid, time2, OUT): """ @@ -71,7 +69,7 @@ def main(grid, time2, OUT): np.sin(lat_rad) * np.sin(d_rad) + np.cos(lat_rad) * np.cos(d_rad) * np.cos(h_rad) ) # SOURCE: Iqbal (1983) - if grid["input_type"] is GridInputType.MATLAB: # only use shading routine in uncoupled runs on a regular grid + if grid["has_shading"]: yl = grid["Ly"] xl = grid["Lx"] diff --git a/src/main.py b/src/main.py index 3a5deb0d..314450b9 100644 --- a/src/main.py +++ b/src/main.py @@ -186,6 +186,13 @@ def main(): C = INIT.init_constants() grid = INIT.init_grid(grid, io, grid_config) + # Ensure shading routine is only used in uncoupled runs on unpartitioned MATLAB grids; + # see https://github.com/EBFMorg/EBFM/issues/11 for details. + if grid["has_shading"]: + assert grid_config.is_partitioned is False, "Shading routine only implemented for unpartitioned grids." + assert grid_config.grid_type is GridInputType.MATLAB, "Shading routine only implemented for MATLAB input grids." + assert coupling_config.defines_coupling() is False, "Shading routine not implemented for coupled runs." + OUT, IN, OUTFILE = INIT.init_initial_conditions(C, grid, io, time2) if coupling_config.defines_coupling(): @@ -265,7 +272,7 @@ def main(): # Write output to files (only in uncoupled run and for unpartitioned grid) if not grid["is_partitioned"] and not coupler.has_coupling: assert ( - grid["input_type"] is GridInputType.MATLAB + grid_config.grid_type is GridInputType.MATLAB ), "Output writing currently only implemented for MATLAB input grids." io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time2, C) pass From 2cf7b25dea379618180d1820896d26740a5b3a35 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 25 Nov 2025 20:24:10 +0100 Subject: [PATCH 15/25] Add pointers to #11. --- src/ebfm/INIT.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index afdadcc3..2d8e13a0 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -184,7 +184,7 @@ def init_grid(grid, io, config: GridConfig): grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells grid["mesh"] = mesh - grid["has_shading"] = False + grid["has_shading"] = False # TODO: see https://github.com/EBFMorg/EBFM/issues/11 # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy elif config.grid_type is GridInputType.ELMER: # Read grid and elevations from Elmer @@ -210,7 +210,7 @@ def init_grid(grid, io, config: GridConfig): # TODO later add slope # grid["slope_x"], grid["slope_y"] = mesh.dzdy, mesh.dzdy grid["mesh"] = mesh - grid["has_shading"] = False + grid["has_shading"] = False # TODO: see https://github.com/EBFMorg/EBFM/issues/11 elif config.grid_type is GridInputType.MATLAB: # Read grid and elevations from example MATLAB file # --------------------------------------------------------------------- # Read and process grid information From 9b044b9ec37ad37484978bb8832c26a983e40cec Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 25 Nov 2025 21:00:57 +0100 Subject: [PATCH 16/25] Apply some changes suggested in review. --- src/ebfm/INIT.py | 10 +++++----- src/ebfm/grid.py | 14 +++++++++++--- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 691573a0..b4cc0188 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -175,12 +175,12 @@ def init_grid(grid, io, config: GridConfig): mesh: Mesh = read_elmer_mesh(mesh_root=config.mesh_file) grid["x"], grid["y"] = mesh.x_vertices, mesh.y_vertices - if args.netcdf_mesh: - grid["z"] = read_dem(args.netcdf_mesh, grid["x"], grid["y"]) + if config.grid_type is GridInputType.CUSTOM: + grid["z"] = read_dem(config.dem_file, grid["x"], grid["y"]) grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! - if args.netcdf_mesh_unstructured: - grid = read_dem_xios(args.netcdf_mesh_unstructured, grid) + if config.grid_type is GridInputType.XIOS_CUSTOM: + grid = read_dem_xios(config.dem_file, grid) grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells grid["slope_x"] = np.zeros_like(grid["x"]) # test values! @@ -314,7 +314,7 @@ def init_grid(grid, io, config: GridConfig): grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2 grid["slope_gamma"] = -grid["slope_gamma"] elif grid["input_type"] is GridInputType.ELMER_XIOS: - grid = read_elmer_xios_grid(grid=grid, gridfile=args.elmer_xios_mesh) + grid = read_elmer_xios_grid(grid=grid, gridfile=config.dem_file) grid["gpsum"] = grid["z"].shape[0] grid["mask"] = (grid["h"] > 1.0) * 1.0 grid["x"] = np.zeros_like(grid["z"]) diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index 5428b33c..ae0581ce 100644 --- a/src/ebfm/grid.py +++ b/src/ebfm/grid.py @@ -6,6 +6,14 @@ class GridInputType(Enum): - MATLAB = "matlab" # .mat grid file with elevation - CUSTOM = "custom" # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file - ELMER = "elmer" # Elmer/Ice mesh file with elevation in z-coordinate + # .mat grid file with elevation + MATLAB = "matlab" + + # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file + CUSTOM = "custom" + + # Elmer/Ice mesh file for xy-coordinates and separate unstructured NetCDF elevation file obtained from XIOS + XIOS_CUSTOM = "xios_custom" + + # Elmer/Ice mesh file with elevation in z-coordinate + ELMER = "elmer" From f222ff976ca36a62771f4c2a9e4c0052425e6442 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 26 Nov 2025 17:03:41 +0100 Subject: [PATCH 17/25] Check for property dem_file instead of grid type. --- src/ebfm/INIT.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index b4cc0188..23c5e369 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -164,7 +164,11 @@ def init_grid(grid, io, config: GridConfig): grid["is_partitioned"] = config.is_partitioned # Read grid from Elmer, elevations from BedMachine - if config.grid_type is GridInputType.CUSTOM or config.grid_type is GridInputType.XIOS_CUSTOM: + if config.dem_file: + assert config.grid_type in { + GridInputType.CUSTOM, + GridInputType.XIOS_CUSTOM, + }, "DEM file can only be specified for CUSTOM or XIOS_CUSTOM grid types." if config.is_partitioned: mesh: Mesh = read_elmer_mesh( mesh_root=config.mesh_file, From a8596994dca119de85e65859648f6ba8a5bbe8b6 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 26 Nov 2025 17:09:00 +0100 Subject: [PATCH 18/25] Set dem_file=None by default. --- src/ebfm/config.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ebfm/config.py b/src/ebfm/config.py index 27a553b5..f159c2d8 100644 --- a/src/ebfm/config.py +++ b/src/ebfm/config.py @@ -59,7 +59,7 @@ class GridConfig: grid_type: GridInputType # Name of the grid used in coupling mesh_file: Path # Path to the grid file - dem_file: Path # Path to the DEM file (only relevant for CUSTOM grid type) + dem_file: Path = None # Path to the DEM file (only relevant for CUSTOM grid type) is_partitioned: bool # Whether the grid is partitioned partition_id: int # Partition ID (only relevant if is_partitioned is True) From c495e6e049656354a0ccc20d50bb2566bf26d330 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 11:28:25 +0100 Subject: [PATCH 19/25] Rename XIOS_CUSTOM to ELMERXIOS. --- src/ebfm/INIT.py | 10 +++++----- src/ebfm/config.py | 2 +- src/ebfm/grid.py | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 23c5e369..8051ca9d 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -165,10 +165,10 @@ def init_grid(grid, io, config: GridConfig): # Read grid from Elmer, elevations from BedMachine if config.dem_file: - assert config.grid_type in { - GridInputType.CUSTOM, - GridInputType.XIOS_CUSTOM, - }, "DEM file can only be specified for CUSTOM or XIOS_CUSTOM grid types." + grid_input_type_supporting_dem = [GridInputType.CUSTOM, GridInputType.ELMERXIOS] + assert ( + config.grid_type in grid_input_type_supporting_dem + ), f"DEM file can only be specified for {grid_input_type_supporting_dem}." if config.is_partitioned: mesh: Mesh = read_elmer_mesh( mesh_root=config.mesh_file, @@ -183,7 +183,7 @@ def init_grid(grid, io, config: GridConfig): grid["z"] = read_dem(config.dem_file, grid["x"], grid["y"]) grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! - if config.grid_type is GridInputType.XIOS_CUSTOM: + if config.grid_type is GridInputType.ELMERXIOS: grid = read_dem_xios(config.dem_file, grid) grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells diff --git a/src/ebfm/config.py b/src/ebfm/config.py index f159c2d8..ea365071 100644 --- a/src/ebfm/config.py +++ b/src/ebfm/config.py @@ -101,7 +101,7 @@ def __init__(self, args: Namespace): self.mesh_file = args.elmer_mesh self.dem_file = args.netcdf_mesh elif args.netcdf_mesh_unstructured and args.elmer_mesh: - self.grid_type = GridInputType.XIOS_CUSTOM + self.grid_type = GridInputType.ELMERXIOS self.mesh_file = args.elmer_mesh self.dem_file = args.netcdf_mesh_unstructured elif args.elmer_mesh: diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index ae0581ce..9f39dfd6 100644 --- a/src/ebfm/grid.py +++ b/src/ebfm/grid.py @@ -13,7 +13,7 @@ class GridInputType(Enum): CUSTOM = "custom" # Elmer/Ice mesh file for xy-coordinates and separate unstructured NetCDF elevation file obtained from XIOS - XIOS_CUSTOM = "xios_custom" + ELMERXIOS = "elmerxios" # Elmer/Ice mesh file with elevation in z-coordinate ELMER = "elmer" From 2101e637362d0056c65ee025a6e7295c008e09e4 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 20:14:40 +0100 Subject: [PATCH 20/25] Change default value. --- src/main.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main.py b/src/main.py index 79554678..e527c265 100644 --- a/src/main.py +++ b/src/main.py @@ -97,7 +97,7 @@ def extract_active_coupling_features(args: argparse.Namespace) -> List[str]: def main(): - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument( "--version", @@ -138,15 +138,15 @@ def main(): time_group.add_argument( "--start-time", type=str, - help="Start time of the simulation in format 'DD-Mon-YYYY HH:MM', e.g., '01-Jan-1979 00:00'.", + help="Start time of the simulation in format 'DD-Mon-YYYY HH:MM'", default="1-Jan-1979 00:00", ) time_group.add_argument( "--end-time", type=str, - help="End time of the simulation in format 'DD-Mon-YYYY HH:MM', e.g., '02-Jan-1979 09:00'.", - default="2-Jan-1979 09:00", + help="End time of the simulation in format 'DD-Mon-YYYY HH:MM'", + default="2-Jan-1979 00:00", ) time_group.add_argument( From 01229e7f0a5d5cd3c216a437b14706abcf921726 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 20:27:20 +0100 Subject: [PATCH 21/25] Store structured/unstructured in grid. --- src/ebfm/INIT.py | 1 + src/ebfm/LOOP_write_to_file.py | 18 +++++++++--------- src/ebfm/config.py | 2 ++ src/main.py | 3 +-- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index e1e83d2d..49d68f24 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -153,6 +153,7 @@ def init_constants(): def init_grid(grid, io, config: GridConfig): grid["is_partitioned"] = config.is_partitioned + grid["is_unstructured"] = config.is_unstructured # Read grid from Elmer, elevations from BedMachine if config.dem_file: diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 3f2c4b1f..0ad898d4 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -7,7 +7,7 @@ from netCDF4 import Dataset, date2num -def main(OUTFILE, io, OUT, grid, t, time, C, gridtype): +def main(OUTFILE, io, OUT, grid, t, time, C): # Specify variables to be written if t == 1: OUTFILE["varsout"] = [ @@ -114,7 +114,7 @@ def save_binary_files(): return True - def save_netCDF_file(gridtype): + def save_netCDF_file(): """ Save model output stored in OUTFILE to a NetCDF file. Converts 1D data to 2D grids. @@ -141,12 +141,12 @@ def save_netCDF_file(gridtype): # Create NetCDF file nc_filepath = os.path.join(io["outdir"], "model_output.nc") io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4") - if gridtype == "unstructured": + if grid["is_unstructured"]: # Define dimensions io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["lat"].shape[0]) # 2D grid rows io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables - elif gridtype == "structured": + else: # structured grid # Define dimensions io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["x_2D"].shape[0]) # 2D grid rows @@ -161,7 +161,7 @@ def save_netCDF_file(gridtype): # Check if variable is a `sub` variable if varname.startswith("sub"): - if gridtype == "unstructured": + if grid["is_unstructured"]: # Define variable as 3D: (time, y, nl) nc_var = io["nc_file"].createVariable( varname, @@ -172,7 +172,7 @@ def save_netCDF_file(gridtype): fill_value=-9999.0, # Fill missing values chunksizes=(1, grid["lat"].shape[0], grid["nl"]), ) - elif gridtype == "structured": + else: # structured grid # Define variable as 4D: (time, y, x, nl) nc_var = io["nc_file"].createVariable( varname, @@ -184,7 +184,7 @@ def save_netCDF_file(gridtype): chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]), ) else: - if gridtype == "unstructured": + if grid["is_unstructured"]: # Define variable as 2D: (time, y) nc_var = io["nc_file"].createVariable( varname, @@ -195,7 +195,7 @@ def save_netCDF_file(gridtype): fill_value=-9999.0, chunksizes=(1, grid["lat"].shape[0]), ) - elif gridtype == "structured": + else: # structured grid # Define variable as 3D: (time, y, x) nc_var = io["nc_file"].createVariable( varname, @@ -253,7 +253,7 @@ def save_netCDF_file(gridtype): if io["output_type"] == 1: save_binary_files() elif io["output_type"] == 2: - save_netCDF_file(gridtype=gridtype) + save_netCDF_file() else: print("Invalid output type. Please choose 1 for binary files or 2 for NetCDF.") diff --git a/src/ebfm/config.py b/src/ebfm/config.py index 9f38b1f3..4a9aa45f 100644 --- a/src/ebfm/config.py +++ b/src/ebfm/config.py @@ -64,6 +64,7 @@ class GridConfig: mesh_file: Path # Path to the grid file dem_file: Path = None # Path to the DEM file (only relevant for CUSTOM grid type) is_partitioned: bool # Whether the grid is partitioned + is_unstructured: bool = False # Whether the grid is unstructured partition_id: int # Partition ID (only relevant if is_partitioned is True) def __init__(self, args: Namespace): @@ -107,6 +108,7 @@ def __init__(self, args: Namespace): self.grid_type = GridInputType.ELMERXIOS self.mesh_file = args.elmer_mesh self.dem_file = args.netcdf_mesh_unstructured + self.is_unstructured = True elif args.elmer_mesh: self.grid_type = GridInputType.ELMER self.mesh_file = args.elmer_mesh diff --git a/src/main.py b/src/main.py index e527c265..74b898fe 100644 --- a/src/main.py +++ b/src/main.py @@ -304,14 +304,13 @@ def main(): # Write output to files (only in uncoupled run and for unpartitioned grid) if not grid["is_partitioned"] and not coupler.has_coupling: if grid_config.grid_type is GridInputType.MATLAB: - gridtype = "structured" # if grid_config.grid_type is GridInputType.MATLAB or # grid['input_type'] is GridInputType.ELMER_XIOS: # ... # if grid_config.grid_type is GridInputType.MATLAB: # assert (grid['input_type'] is GridInputType.MATLAB), \ # "Output writing currently only implemented for MATLAB and ELMER_XIOS input grids." - io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time, C, gridtype=gridtype) + io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time, C) else: logger.warning("Skipping writing output to file for Elmer input grids.") elif grid["is_partitioned"] or coupler.has_coupling: From e7d2807a0811d0efb2d4ebbe0b647eb35c5d05c3 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 20:41:54 +0100 Subject: [PATCH 22/25] Reduce code duplication. --- src/ebfm/LOOP_write_to_file.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 0ad898d4..64db829d 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -141,17 +141,17 @@ def save_netCDF_file(): # 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"].createDimension("time", None) # Unlimited time dimension + if grid["is_unstructured"]: # Define dimensions - io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["lat"].shape[0]) # 2D grid rows - io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables else: # structured grid # Define dimensions - io["nc_file"].createDimension("time", None) # Unlimited time dimension io["nc_file"].createDimension("y", grid["x_2D"].shape[0]) # 2D grid rows io["nc_file"].createDimension("x", grid["x_2D"].shape[1]) # 2D grid columns - io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables + + io["nc_file"].createDimension("nl", grid["nl"]) # Vertical layers for `sub` variables # Define standard output variables for entry in OUTFILE["varsout"]: From db2c65267d620e5b5a100dbaf3b5e89e68055319 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 20:59:21 +0100 Subject: [PATCH 23/25] Reduce code duplication and restructure. --- src/ebfm/LOOP_write_to_file.py | 81 ++++++++++++++-------------------- 1 file changed, 32 insertions(+), 49 deletions(-) diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 64db829d..1d6d1e18 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -156,60 +156,43 @@ def save_netCDF_file(): # Define standard output variables for entry in OUTFILE["varsout"]: varname = entry[0] - var_units = entry[1] - var_desc = entry[3] - # Check if variable is a `sub` variable - if varname.startswith("sub"): - if grid["is_unstructured"]: + dimensions: tuple + chunksizes: tuple + + if grid["is_unstructured"]: + # Check if variable is a `sub` variable + if varname.startswith("sub"): # Define variable as 3D: (time, y, nl) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y", "nl"), - zlib=True, - complevel=4, - fill_value=-9999.0, # Fill missing values - chunksizes=(1, grid["lat"].shape[0], grid["nl"]), - ) - else: # structured grid - # Define variable as 4D: (time, y, x, nl) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y", "x", "nl"), - zlib=True, - complevel=4, - fill_value=-9999.0, # Fill missing values - chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]), - ) - else: - if grid["is_unstructured"]: + dimensions = ("time", "y", "nl") + chunksizes = (1, grid["lat"].shape[0], grid["nl"]) + else: # Define variable as 2D: (time, y) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y"), - zlib=True, - complevel=4, - fill_value=-9999.0, - chunksizes=(1, grid["lat"].shape[0]), - ) - else: # structured grid + dimensions = ("time", "y") + chunksizes = (1, grid["lat"].shape[0]) + else: + # Check if variable is a `sub` variable + if varname.startswith("sub"): + # Define variable as 4D: (time, y, x, nl) + dimensions = ("time", "y", "x", "nl") + chunksizes = (1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]) + else: # Define variable as 3D: (time, y, x) - nc_var = io["nc_file"].createVariable( - varname, - np.float32, - ("time", "y", "x"), - zlib=True, - complevel=4, - fill_value=-9999.0, - chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1]), - ) - + dimensions = ("time", "y", "x") + chunksizes = (1, grid["x_2D"].shape[0], grid["x_2D"].shape[1]) + + nc_var = io["nc_file"].createVariable( + varname=varname, + datatype=np.float32, + dimensions=dimensions, + zlib=True, + complevel=4, + fill_value=-9999.0, # Fill missing values + chunksizes=chunksizes, + ) # Assign metadata - nc_var.units = var_units - nc_var.description = var_desc + nc_var.units = entry[1] + nc_var.description = entry[3] # Define a time variable to track simulation steps nc_time = io["nc_file"].createVariable("time", np.float64, ("time",), zlib=True, fill_value=-9999.0) From 48f8008b5a320059622a2f5d10cc6885068b0b59 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Tue, 2 Dec 2025 21:25:07 +0100 Subject: [PATCH 24/25] Cleanup. --- src/ebfm/LOOP_write_to_file.py | 6 ++++-- src/main.py | 15 +++++---------- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 1d6d1e18..f7e35df7 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -156,6 +156,8 @@ def save_netCDF_file(): # Define standard output variables for entry in OUTFILE["varsout"]: varname = entry[0] + var_units = entry[1] + var_desc = entry[3] dimensions: tuple chunksizes: tuple @@ -191,8 +193,8 @@ def save_netCDF_file(): chunksizes=chunksizes, ) # Assign metadata - nc_var.units = entry[1] - nc_var.description = entry[3] + nc_var.units = var_units + nc_var.description = var_desc # Define a time variable to track simulation steps nc_time = io["nc_file"].createVariable("time", np.float64, ("time",), zlib=True, fill_value=-9999.0) diff --git a/src/main.py b/src/main.py index 74b898fe..379ae278 100644 --- a/src/main.py +++ b/src/main.py @@ -34,8 +34,8 @@ log_levels = { "file": logging.DEBUG, # log level for logging to file - # 0: logging.INFO, # log level for rank 0 - 0: logging.DEBUG, # to log other ranks to console define log level here + 0: logging.INFO, # log level for rank 0 + # 1: logging.DEBUG, # to log other ranks to console define log level here } setup_logging(log_levels=log_levels) @@ -296,20 +296,15 @@ def main(): logger.debug("Done.") logger.debug("Received the following data from Elmer/Ice:", data_from_elmer) - # IN['h'] = data_from_elmer['h'] - # grid['z'] = IN['h'][0].ravel() + IN["h"] = data_from_elmer["h"] + grid["z"] = IN["h"][0].ravel() + # TODO add gradient field later # IN['dhdx'] = data_from_elmer('dhdx') # IN['dhdy'] = data_from_elmer('dhdy') # Write output to files (only in uncoupled run and for unpartitioned grid) if not grid["is_partitioned"] and not coupler.has_coupling: if grid_config.grid_type is GridInputType.MATLAB: - # if grid_config.grid_type is GridInputType.MATLAB or - # grid['input_type'] is GridInputType.ELMER_XIOS: - # ... - # if grid_config.grid_type is GridInputType.MATLAB: - # assert (grid['input_type'] is GridInputType.MATLAB), \ - # "Output writing currently only implemented for MATLAB and ELMER_XIOS input grids." io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time, C) else: logger.warning("Skipping writing output to file for Elmer input grids.") From 67f59f505f329a465c0e7f66cd3fe1756d1801ff Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 3 Dec 2025 11:57:36 +0100 Subject: [PATCH 25/25] Add Changelog entry. --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f0271ae9..519cd8e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,8 +6,9 @@ SPDX-License-Identifier: CC-BY-4.0 # develop +* Support new input mesh format. EBFM now accepts Elmer/Ice mesh file for xy-coordinates and separate unstructured NetCDF elevation file obtained from XIOS. To use this feature please provide `--elmer-mesh` together with the new option `--netcdf-mesh-unstructured`. See https://github.com/EBFMorg/EBFM/pull/12. * Use `setuptools_scm` as backend for `--version` information. See https://github.com/EBFMorg/EBFM/pull/46. -* Remove `pathlib` from requirements, because this can lead to a bug. https://github.com/EBFMorg/EBFM/pull/48 +* Remove `pathlib` from requirements, because this can lead to a bug. https://github.com/EBFMorg/EBFM/pull/48. # v0.1.0