From b9ade48ea281adbbd6bf239440718b73986f0b54 Mon Sep 17 00:00:00 2001 From: Toby Searle <14909402+twsearle@users.noreply.github.com> Date: Thu, 25 Jan 2024 17:05:22 +0000 Subject: [PATCH] Read all ORCA field data onto a single processor (#71) * Add nemo_io::ReadServer and const correctness * Refactor and remove unnecessary methods * add Doxygen docs coverage * cpplinting pass --- src/orca-jedi/CMakeLists.txt | 3 + src/orca-jedi/nemo_io/NemoFieldReader.cc | 432 ++---------------- src/orca-jedi/nemo_io/NemoFieldReader.h | 31 +- src/orca-jedi/nemo_io/OrcaIndex.h | 58 +++ src/orca-jedi/nemo_io/ReadServer.cc | 176 +++++++ src/orca-jedi/nemo_io/ReadServer.h | 59 +++ src/orca-jedi/state/StateIOUtils.cc | 17 +- src/tests/orca-jedi/CMakeLists.txt | 4 +- .../orca-jedi/test_nemo_io_field_reader.cc | 76 --- .../orca-jedi/test_nemo_io_field_writer.cc | 21 +- .../test_nemo_io_orca_field_reader.cc | 234 ---------- .../orca-jedi/test_nemo_io_read_server.cc | 158 +++++++ 12 files changed, 520 insertions(+), 749 deletions(-) create mode 100644 src/orca-jedi/nemo_io/OrcaIndex.h create mode 100644 src/orca-jedi/nemo_io/ReadServer.cc create mode 100644 src/orca-jedi/nemo_io/ReadServer.h delete mode 100644 src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc create mode 100644 src/tests/orca-jedi/test_nemo_io_read_server.cc diff --git a/src/orca-jedi/CMakeLists.txt b/src/orca-jedi/CMakeLists.txt index 936aa5d..cbf56c5 100644 --- a/src/orca-jedi/CMakeLists.txt +++ b/src/orca-jedi/CMakeLists.txt @@ -27,6 +27,9 @@ nemo_io/NemoFieldReader.h nemo_io/NemoFieldReader.cc nemo_io/NemoFieldWriter.h nemo_io/NemoFieldWriter.cc +nemo_io/OrcaIndex.h +nemo_io/ReadServer.h +nemo_io/ReadServer.cc utilities/OrcaModelTraits.h variablechanges/VariableChangeParameters.h variablechanges/VariableChange.h diff --git a/src/orca-jedi/nemo_io/NemoFieldReader.cc b/src/orca-jedi/nemo_io/NemoFieldReader.cc index 73bc787..4bd4333 100644 --- a/src/orca-jedi/nemo_io/NemoFieldReader.cc +++ b/src/orca-jedi/nemo_io/NemoFieldReader.cc @@ -7,7 +7,7 @@ #include // Using Lynton Appel's netcdf-cxx4 from // https://github.com/Unidata/netcdf-cxx4 - +// #include #include #include @@ -20,44 +20,19 @@ #include "oops/util/Logger.h" #include "oops/util/Duration.h" -#include "atlas/field.h" #include "atlas-orca/grid/OrcaGrid.h" +#include "orca-jedi/nemo_io/OrcaIndex.h" namespace orcamodel { namespace { -struct IndexGlbArray { - int32_t ix_glb_max; - int32_t iy_glb_max; - int32_t glbarray_offset; - int32_t glbarray_jstride; - int32_t nx_halo_WE; - int32_t ny_halo_NS; - - explicit IndexGlbArray(const atlas::OrcaGrid& orcaGrid) { - iy_glb_max = orcaGrid.ny() + orcaGrid.haloNorth() - 1; - ix_glb_max = orcaGrid.nx() + orcaGrid.haloEast() - 1; - - nx_halo_WE = orcaGrid.nx() + orcaGrid.haloEast() + orcaGrid.haloWest(); - ny_halo_NS = orcaGrid.ny() + orcaGrid.haloNorth() - + orcaGrid.haloSouth(); - - // vector of local indices: necessary for remote indices of ghost nodes - int iy_glb_min = -orcaGrid.haloSouth(); - int ix_glb_min = -orcaGrid.haloWest(); - glbarray_offset = -(nx_halo_WE * iy_glb_min) - ix_glb_min; - glbarray_jstride = nx_halo_WE; - } - - int64_t operator()(int i, int j) { - ATLAS_ASSERT(i <= ix_glb_max, - std::to_string(i) + " > " + std::to_string(ix_glb_max)); - ATLAS_ASSERT(j <= iy_glb_max, - std::to_string(j) + " > " + std::to_string(iy_glb_max)); - return glbarray_offset + j * glbarray_jstride + i; - } -}; +/// \brief Search the netCDF file for a dimension matching a name from a list of possible names. +/// \param ncFile The netCDF file. +/// \param check_dim_for_dimvar Set to true to ensure the dimension has +/// a corresponding dimension variable. +/// \param possible_names A vector of all possible names for the variable. +/// \return The dimension name. std::string find_nc_var_name(const netCDF::NcFile& ncFile, const bool check_dim_for_dimvar, const std::vector& possible_names) { @@ -89,7 +64,7 @@ std::string find_nc_var_name(const netCDF::NcFile& ncFile, } } // namespace -NemoFieldReader::NemoFieldReader(eckit::PathName& filename) +NemoFieldReader::NemoFieldReader(const eckit::PathName& filename) : ncFile(nullptr), datetimes_() { oops::Log::debug() << "orcamodel::NemoFieldReader::NemoFieldReader filename: " << filename.fullName().asString() << std::endl; @@ -117,7 +92,10 @@ NemoFieldReader::NemoFieldReader(eckit::PathName& filename) read_datetimes(); } -size_t NemoFieldReader::read_dim_size(const std::string& name) { +/// \brief Read the dimension size for a given netCDF dimension specified by name +/// \param name The name of the netCDF dimension +/// \return size The size +size_t NemoFieldReader::read_dim_size(const std::string& name) const { auto dim = ncFile->getDim(name); if (dim.isNull()) { std::ostringstream err_stream; @@ -132,7 +110,7 @@ size_t NemoFieldReader::read_dim_size(const std::string& name) { return dim.getSize(); } -/// \brief retrieve the datetime for each time index in file. +/// \brief Update the datetimes_ in the object to contain times for each time index in file. void NemoFieldReader::read_datetimes() { // read time indices from file size_t n_times = read_dim_size(time_dimvar_name_); @@ -201,7 +179,9 @@ void NemoFieldReader::read_datetimes() { /// \brief Read the _FillValue for a variable, defaulting to the minimum value /// for the datatype. -template T NemoFieldReader::read_fillvalue(const std::string& name) { +/// \param name Name of the netCDF variable containing the _FillValue attribute to retrieve. +/// \return The fill value for this netCDF variable +template T NemoFieldReader::read_fillvalue(const std::string& name) const { netCDF::NcVar nc_var = ncFile->getVar(name); if (nc_var.isNull()) { std::ostringstream err_stream; @@ -230,15 +210,17 @@ template T NemoFieldReader::read_fillvalue(const std::string& name) { return fillvalue; } -template int NemoFieldReader::read_fillvalue(const std::string& name); -template float NemoFieldReader::read_fillvalue(const std::string& name); +template int NemoFieldReader::read_fillvalue(const std::string& name) const; +template float NemoFieldReader::read_fillvalue(const std::string& name) const; template double NemoFieldReader::read_fillvalue( - const std::string& name); + const std::string& name) const; /// \brief get the time dimension index corresponding to the nearest datetime /// to a target datetime. +/// \param tgt_datetime Search for the index of the time slice in the file nearest this datetime +/// \return The index of the nearest time slice in the file size_t NemoFieldReader::get_nearest_datetime_index( - const util::DateTime& tgt_datetime) { + const util::DateTime& tgt_datetime) const { int64_t time_diff = INT64_MAX; size_t indx = 0; util::Duration duration; @@ -254,7 +236,9 @@ size_t NemoFieldReader::get_nearest_datetime_index( return indx; } -std::vector NemoFieldReader::read_locs() { +/// \brief Read the latitude longitude locations of all points in a NEMO field file +/// \return A vector of XY points of all nodes in the field. +std::vector NemoFieldReader::read_locs() const { try { size_t nx = read_dim_size("x"); size_t ny = read_dim_size("y"); @@ -302,8 +286,13 @@ std::vector NemoFieldReader::read_locs() { } } +/// \brief Read all data for a given variable at a level and time. +/// \param varname The name of the variable. +/// \param t_indx The time index of the slice. +/// \param z_indx The vertical index of the slice. +/// \return a vector of the variable data. std::vector NemoFieldReader::read_var_slice(const std::string& varname, - const size_t t_indx, const size_t z_indx) { + const size_t t_indx, const size_t z_indx) const { oops::Log::trace() << "orcamodel::NemoFieldReader::read_var_slice" << std::endl; try { @@ -313,7 +302,7 @@ std::vector NemoFieldReader::read_var_slice(const std::string& varname, netCDF::NcVar nc_var = ncFile->getVar(varname); if (nc_var.isNull()) { std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar '" + err_stream << "orcamodel::NemoFieldReader::read_var_slice ncVar '" << varname << "' is not present in NetCDF file"; throw eckit::BadValue(err_stream.str(), Here()); } @@ -329,7 +318,7 @@ std::vector NemoFieldReader::read_var_slice(const std::string& varname, nc_var.getVar({0, 0}, {ny, nx}, var_data.data()); } else { std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar '" + err_stream << "orcamodel::NemoFieldReader::read_var_slice ncVar '" << varname << "' has " << n_dims << " dimensions."; throw eckit::BadValue(err_stream.str(), Here()); } @@ -338,144 +327,30 @@ std::vector NemoFieldReader::read_var_slice(const std::string& varname, } catch(netCDF::exceptions::NcException& e) { std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var varname: " - << varname << " NetCDF exception: " << std::endl << e.what(); - throw eckit::ReadError(err_stream.str(), Here()); - } -} - - void NemoFieldReader::read_volume_var(const std::string& varname, - const size_t t_indx, atlas::array::ArrayView& field_view) { - oops::Log::trace() << "orcamodel::NemoFieldReader::read_volume_var" - << std::endl; - try { - size_t nx = read_dim_size("x"); - size_t ny = read_dim_size("y"); - size_t nz = read_dim_size(z_dimvar_name_); - size_t nlevels = field_view.shape(1); - - if (field_view.shape(0) != nx*ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var field_view 1st" - << " dimension does not match horizontal dimensions" - << " for varname " << varname; - throw eckit::BadValue(err_stream.str(), Here()); - } - - if (nlevels > nz) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var field_view 2nd" - << " dimension " << nlevels << " is larger than NetCDF file" - << " z dimension " << nz << " for varname " << varname; - throw eckit::BadValue(err_stream.str(), Here()); - } - - netCDF::NcVar nc_var = ncFile->getVar(varname); - if (nc_var.isNull()) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var ncVar '" - << varname << "' is not present in NetCDF file"; - throw eckit::BadValue(err_stream.str(), Here()); - } - - std::vector buffer(nx*ny*nlevels); - - size_t n_dims = nc_var.getDimCount(); - std::string first_dim_name = nc_var.getDim(0).getName(); - if (n_dims == 4) { - nc_var.getVar({t_indx, 0, 0, 0}, {1, nlevels, ny, nx}, buffer.data()); - } else if (n_dims == 3 && first_dim_name == z_dimvar_name_) { - nc_var.getVar({0, 0, 0}, {nlevels, ny, nx}, buffer.data()); - } else { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var ncVar '" - << varname << "' has " << n_dims << " dimensions."; - throw eckit::BadValue(err_stream.str(), Here()); - } - - // in atlas fields the levels indices change the fastest, so we need to - // swap the indexing order from the netCDF data. - for (int n = 0; n < nx*ny; ++n) { - for (int k = 0; k < nlevels; ++k) { - field_view(n, k) = buffer[k*nx*ny + n]; - } - } - } catch(netCDF::exceptions::NcException& e) - { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var varname: " - << varname << " NetCDF exception: " << std::endl << e.what(); - throw eckit::ReadError(err_stream.str(), Here()); - } -} - -void NemoFieldReader::read_vertical_var(const std::string& varname, - atlas::array::ArrayView& field_view) { - oops::Log::trace() << "orcamodel::NemoFieldReader::read_vertical_var" - << std::endl; - try { - size_t nx = read_dim_size("x"); - size_t ny = read_dim_size("y"); - size_t nz = read_dim_size(z_dimvar_name_); - size_t nlevels = field_view.shape(1); - - if (nlevels > nz) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_vertical_var field_view " - << "2nd dimension is larger than NetCDF file z dimension "; - throw eckit::BadValue(err_stream.str(), Here()); - } - - netCDF::NcVar nc_var = ncFile->getVar(varname); - if (nc_var.isNull()) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_vertical_var ncVar '" - << varname << "' is not present in NetCDF file"; - throw eckit::BadValue(err_stream.str(), Here()); - } - - std::vector buffer(nlevels); - - size_t n_dims = nc_var.getDimCount(); - std::string first_dim_name = nc_var.getDim(0).getName(); - if (n_dims == 1 && first_dim_name == z_dimvar_name_) { - nc_var.getVar({0}, {nlevels}, buffer.data()); - } else { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_vertical_var ncVar '" - << varname << "' has " << n_dims << " dimensions."; - throw eckit::BadValue(err_stream.str(), Here()); - } - - // Store the data in an atlas 3D field - inefficient but flexible - for (size_t inode = 0; inode < field_view.shape(0); ++inode) { - for (int k = 0; k < nlevels; ++k) { - field_view(inode, k) = buffer[k]; - } - } - } catch(netCDF::exceptions::NcException& e) - { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_vertical_var varname: " + err_stream << "orcamodel::NemoFieldReader::read_var_slice varname: " << varname << " NetCDF exception: " << std::endl << e.what(); throw eckit::ReadError(err_stream.str(), Here()); } } -void NemoFieldReader::read_vertical_var(const std::string& varname, - const atlas::Mesh& mesh, atlas::array::ArrayView& field_view) { +/// \brief Read a 1D variable containing level depth information. +/// \param varname The name of the variable. +/// \param n_levels The number of levels to read (beginning from the surface). +/// \return a vector of the depth values. +std::vector NemoFieldReader::read_vertical_var( + const std::string& varname, + const size_t nlevels) const { oops::Log::trace() << "orcamodel::NemoFieldReader::read_vertical_var" << std::endl; try { size_t nx = read_dim_size("x"); size_t ny = read_dim_size("y"); size_t nz = read_dim_size(z_dimvar_name_); - size_t nlevels = field_view.shape(1); if (nlevels > nz) { std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_vertical_var field_view " - << "2nd dimension is larger than NetCDF file z dimension "; + err_stream << "orcamodel::NemoFieldReader::read_vertical_var num levels " + << "is larger than NetCDF file z dimension "; throw eckit::BadValue(err_stream.str(), Here()); } @@ -500,14 +375,7 @@ void NemoFieldReader::read_vertical_var(const std::string& varname, throw eckit::BadValue(err_stream.str(), Here()); } - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - // Store the data in an atlas 3D field - inefficient but flexible - for (size_t inode = 0; inode < field_view.shape(0); ++inode) { - for (int k = 0; k < nlevels; ++k) { - if (ghost(inode)) continue; - field_view(inode, k) = buffer[k]; - } - } + return buffer; } catch(netCDF::exceptions::NcException& e) { std::ostringstream err_stream; @@ -517,214 +385,4 @@ void NemoFieldReader::read_vertical_var(const std::string& varname, } } -void NemoFieldReader::read_surf_var(const std::string& varname, - const size_t t_indx, atlas::array::ArrayView& field_view) { - oops::Log::trace() << "orcamodel::NemoFieldReader::read_surf_var" - << std::endl; - try { - size_t nx = read_dim_size("x"); - size_t ny = read_dim_size("y"); - - if (field_view.size() != nx*ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var field_view " - << "dimensions do not match dimensions in netCDF file "; - throw eckit::BadValue(err_stream.str(), Here()); - } - - netCDF::NcVar nc_var = ncFile->getVar(varname); - if (nc_var.isNull()) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar '" - << varname << "' is not present in NetCDF file"; - throw eckit::BadValue(err_stream.str(), Here()); - } - - size_t n_dims = nc_var.getDimCount(); - if (n_dims == 4) { - nc_var.getVar({t_indx, 0, 0, 0}, {1, 1, ny, nx}, field_view.data()); - } else if (n_dims == 3) { - nc_var.getVar({t_indx, 0, 0}, {1, ny, nx}, field_view.data()); - } else if (n_dims == 2) { - nc_var.getVar({0, 0}, {ny, nx}, field_view.data()); - } else { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar '" - << varname << "' has " << n_dims << " dimensions."; - throw eckit::BadValue(err_stream.str(), Here()); - } - } catch(netCDF::exceptions::NcException& e) - { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var varname: " - << varname << " NetCDF exception: " << std::endl << e.what(); - throw eckit::ReadError(err_stream.str(), Here()); - } -} - -void NemoFieldReader::read_surf_var(const std::string& varname, - const atlas::Mesh& mesh, const size_t t_indx, - atlas::array::ArrayView& field_view) { - oops::Log::trace() << "orcamodel::NemoFieldReader::read_surf_var" - << std::endl; - try { - size_t nx = read_dim_size("x"); - size_t ny = read_dim_size("y"); - - const atlas::OrcaGrid orcaGrid = atlas::OrcaGrid(mesh.grid()); - if (!orcaGrid) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var" - << " only reads ORCA grid data. " << mesh.grid().name() - << " is not an ORCA grid." << std::endl; - throw eckit::BadValue(err_stream.str(), Here()); - } - - IndexGlbArray index_glbarray(orcaGrid); - if (index_glbarray.nx_halo_WE != nx - || index_glbarray.ny_halo_NS != ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar" - << " grid dimensions don't match file dimensions grid: (" - << index_glbarray.nx_halo_WE << ", " - << index_glbarray.ny_halo_NS << ") file: " << nx << ", " - << ny << ")" << std::endl; - throw eckit::BadValue(err_stream.str()); - } - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - if (mesh.nodes().size() > nx * ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_surf_var ncVar" - << " number of mesh nodes " << "larger than the netCDF file" - << " dimensions, " << nx * ny - << " double check the grids match." << std::endl; - throw eckit::BadValue(err_stream.str()); - } - - netCDF::NcVar nc_var = ncFile->getVar(varname); - if (nc_var.isNull()) { - throw eckit::UserError( - std::string("orcamodel::NemoFieldReader::read_surf_var ncVar ") - + varname + " is not present in NetCDF file", Here()); - } - - std::vector var_data(nx * ny); - size_t n_dims = nc_var.getDimCount(); - if (n_dims == 4) { - nc_var.getVar({t_indx, 0, 0, 0}, {1, 1, ny, nx}, var_data.data()); - } else if (n_dims == 3) { - nc_var.getVar({t_indx, 0, 0}, {1, ny, nx}, var_data.data()); - } else if (n_dims == 2) { - nc_var.getVar({0, 0}, {ny, nx}, var_data.data()); - } else { - throw eckit::UserError( - std::string("orcamodel::NemoFieldReader::read_surf_var ") - + "ncVar " + varname + " has an unreadable number of " - + "dimensions: " + std::to_string(n_dims), Here()); - } - - for (size_t inode = 0; inode < field_view.size(); ++inode) { - if (ghost(inode)) continue; - double data = var_data[index_glbarray(ij(inode, 0), ij(inode, 1))]; - field_view(inode, 0) = data; - } - } - catch (netCDF::exceptions::NcException& e) { - throw eckit::FailedLibraryCall("NetCDF", - "orcamodel::NemoFieldReader::read_surf_var", e.what(), Here()); - } -} - -void NemoFieldReader::read_volume_var(const std::string& varname, - const atlas::Mesh& mesh, const size_t t_indx, - atlas::array::ArrayView& field_view) { - try { - size_t nx = read_dim_size("x"); - size_t ny = read_dim_size("y"); - size_t nz = read_dim_size(z_dimvar_name_); - size_t nlevels = field_view.shape(1); - - const atlas::OrcaGrid orcaGrid = atlas::OrcaGrid(mesh.grid()); - if (!orcaGrid) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var" - << " only reads ORCA grid data. " << mesh.grid().name() - << " is not an ORCA grid." << std::endl; - throw eckit::BadValue(err_stream.str(), Here()); - } - - IndexGlbArray index_glbarray(orcaGrid); - if (index_glbarray.nx_halo_WE != nx || index_glbarray.ny_halo_NS != ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var ncVar" - << " grid dimensions don't match file dimensions grid: (" - << index_glbarray.nx_halo_WE << ", " - << index_glbarray.ny_halo_NS << ") file: " << nx << ", " - << ny << ")" << std::endl; - throw eckit::BadValue(err_stream.str()); - } - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - if (mesh.nodes().size() > nx * ny) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var field_view 1st" - << " dimension does not match horizontal dimensions" - << " for varname " << varname; - throw eckit::BadValue(err_stream.str(), Here()); - } - - if (nlevels > nz) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var field_view 2nd" - << " dimension " << nlevels << " is larger than NetCDF file" - << " z dimension " << nz << " for varname " << varname; - throw eckit::BadValue(err_stream.str(), Here()); - } - - netCDF::NcVar nc_var = ncFile->getVar(varname); - if (nc_var.isNull()) { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var ncVar '" - << varname << "' is not present in NetCDF file"; - throw eckit::BadValue(err_stream.str(), Here()); - } - - std::vector buffer(nx*ny*nlevels); - - size_t n_dims = nc_var.getDimCount(); - std::string first_dim_name = nc_var.getDim(0).getName(); - if (n_dims == 4) { - nc_var.getVar({t_indx, 0, 0, 0}, {1, nlevels, ny, nx}, buffer.data()); - } else if (n_dims == 3 && first_dim_name == z_dimvar_name_) { - nc_var.getVar({0, 0, 0}, {nlevels, ny, nx}, buffer.data()); - } else { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var ncVar '" - << varname << "' has " << n_dims << " dimensions."; - throw eckit::BadValue(err_stream.str(), Here()); - } - - // in atlas fields the levels indices change the fastest, so we need to - // swap the indexing order from the netCDF data. - const size_t numNodes = field_view.shape(0); - atlas_omp_for(int k = 0; k < nlevels; ++k) { - for (size_t inode = 0; inode < numNodes; ++inode) { - if (ghost(inode)) continue; - field_view(inode, k) = - buffer[k*nx*ny + index_glbarray(ij(inode, 0), ij(inode, 1))]; - } - } - } catch(netCDF::exceptions::NcException& e) - { - std::ostringstream err_stream; - err_stream << "orcamodel::NemoFieldReader::read_volume_var varname: " - << varname << " NetCDF exception: " << std::endl << e.what(); - throw eckit::ReadError(err_stream.str(), Here()); - } -} } // namespace orcamodel diff --git a/src/orca-jedi/nemo_io/NemoFieldReader.h b/src/orca-jedi/nemo_io/NemoFieldReader.h index 159e73e..205790c 100644 --- a/src/orca-jedi/nemo_io/NemoFieldReader.h +++ b/src/orca-jedi/nemo_io/NemoFieldReader.h @@ -18,9 +18,6 @@ #include "atlas/runtime/Exception.h" #include "atlas/util/Point.h" -#include "atlas/field.h" -#include "atlas/mesh.h" -#include "atlas/array/ArrayView.h" namespace orcamodel { @@ -29,29 +26,17 @@ class NemoFieldReader : private util::ObjectCounter { public: static const std::string classname() {return "orcamodel::NemoFieldReader";} - explicit NemoFieldReader(eckit::PathName& filename); + explicit NemoFieldReader(const eckit::PathName& filename); - std::vector read_locs(); - size_t read_dim_size(const std::string& name); void read_datetimes(); - size_t get_nearest_datetime_index(const util::DateTime& datetime); - template T read_fillvalue(const std::string& name); + std::vector read_locs() const; + size_t read_dim_size(const std::string& name) const; + size_t get_nearest_datetime_index(const util::DateTime& datetime) const; + template T read_fillvalue(const std::string& name) const; std::vector read_var_slice(const std::string& varname, - const size_t t_indx, const size_t z_indx); - void read_surf_var(const std::string& varname, const size_t t_indx, - atlas::array::ArrayView& field_view); - void read_volume_var(const std::string& varname, const size_t t_indx, - atlas::array::ArrayView& field_view); - void read_vertical_var(const std::string& varname, - atlas::array::ArrayView& field_view); - void read_vertical_var(const std::string& varname, const atlas::Mesh& mesh, - atlas::array::ArrayView& field_view); - - void read_surf_var(const std::string& varname, const atlas::Mesh& mesh, - const size_t t_indx, atlas::array::ArrayView& field_view); - void read_volume_var(const std::string& varname, - const atlas::Mesh& mesh, const size_t t_indx, - atlas::array::ArrayView& field_view); + const size_t t_indx, const size_t z_indx) const; + std::vector read_vertical_var(const std::string& varname, + const size_t nlevels) const; private: NemoFieldReader() : ncFile() {} diff --git a/src/orca-jedi/nemo_io/OrcaIndex.h b/src/orca-jedi/nemo_io/OrcaIndex.h new file mode 100644 index 0000000..7d26590 --- /dev/null +++ b/src/orca-jedi/nemo_io/OrcaIndex.h @@ -0,0 +1,58 @@ +/* + * (C) British Crown Copyright 2024 Met Office + */ + +#pragma once + +#include +#include +#include +#include + +#include "eckit/exception/Exceptions.h" + +#include "oops/util/Logger.h" + +#include "atlas-orca/grid/OrcaGrid.h" + +namespace orcamodel { + +/// \brief Indexer into a 1D array of an ORCA model field, match each point to the corresponding +/// point in the atlas-orca ij space. +struct OrcaIndex { + int32_t ix_glb_max; + int32_t iy_glb_max; + int32_t glbarray_offset; + int32_t glbarray_jstride; + int32_t nx_halo_WE; + int32_t ny_halo_NS; + + explicit OrcaIndex(const atlas::OrcaGrid& orcaGrid) { + iy_glb_max = orcaGrid.ny() + orcaGrid.haloNorth() - 1; + ix_glb_max = orcaGrid.nx() + orcaGrid.haloEast() - 1; + + nx_halo_WE = orcaGrid.nx() + orcaGrid.haloEast() + orcaGrid.haloWest(); + ny_halo_NS = orcaGrid.ny() + orcaGrid.haloNorth() + + orcaGrid.haloSouth(); + + // vector of local indices: necessary for remote indices of ghost nodes + int iy_glb_min = -orcaGrid.haloSouth(); + int ix_glb_min = -orcaGrid.haloWest(); + glbarray_offset = -(nx_halo_WE * iy_glb_min) - ix_glb_min; + glbarray_jstride = nx_halo_WE; + } + + /// \brief Index of a 1D array corresponding to point i, j + /// \param i + /// \param j + /// \return index of a matching 1D array + int64_t operator()(const int i, const int j) const { + ATLAS_ASSERT(i <= ix_glb_max, + std::to_string(i) + " > " + std::to_string(ix_glb_max)); + ATLAS_ASSERT(j <= iy_glb_max, + std::to_string(j) + " > " + std::to_string(iy_glb_max)); + return glbarray_offset + j * glbarray_jstride + i; + } +}; + +} // namespace orcamodel diff --git a/src/orca-jedi/nemo_io/ReadServer.cc b/src/orca-jedi/nemo_io/ReadServer.cc new file mode 100644 index 0000000..767bd42 --- /dev/null +++ b/src/orca-jedi/nemo_io/ReadServer.cc @@ -0,0 +1,176 @@ +/* + * (C) British Crown Copyright 2024 Met Office + */ + + + +#include "orca-jedi/nemo_io/ReadServer.h" +#include +#include +#include + +#include "eckit/exception/Exceptions.h" + +namespace orcamodel { + +ReadServer::ReadServer(const eckit::PathName& file_path, const atlas::Mesh& mesh) : + mesh_(mesh), + index_glbarray_(atlas::OrcaGrid(mesh.grid())) { + if (myrank == mpiroot) { + reader_ = std::make_unique(file_path); + } +} + +/// \brief Read in a 2D horizontal slice of variable data on the root processor only +/// \param var_name +/// \param t_index Index of the time slice. +/// \param z_index Index of the vertical slice. +/// \param buffer Vector to store the data. +void ReadServer::read_var_on_root(const std::string& var_name, + const size_t t_index, + const size_t z_index, + std::vector& buffer) const { + size_t size = index_glbarray_.nx_halo_WE * index_glbarray_.ny_halo_NS; + if (myrank == mpiroot) { + buffer = reader_->read_var_slice(var_name, t_index, z_index); + } else { + buffer.resize(size); + } +} + +/// \brief Read 1D vertical variable data on the root processor only +/// \param var_name NetCDF name of the vertical variable. +/// \param n_levels Number of levels to read from the file +/// \param buffer Vector to store the data. +void ReadServer::read_vertical_var_on_root(const std::string& var_name, + const size_t n_levels, + std::vector& buffer) const { + size_t size = index_glbarray_.nx_halo_WE * index_glbarray_.ny_halo_NS; + if (myrank == mpiroot) { + buffer = reader_->read_vertical_var(var_name, n_levels); + } else { + buffer.resize(size); + } +} + +/// \brief Move data from a buffer into an atlas arrayView. +/// \param buffer Vector of data to read +/// \param z_index Index of the vertical slice. +/// \param field_view View into the atlas field to store the data. +void ReadServer::fill_field(const std::vector& buffer, + const size_t z_index, + atlas::array::ArrayView& field_view) const { + auto ghost = atlas::array::make_view(this->mesh_.nodes().ghost()); + auto ij = atlas::array::make_view(this->mesh_.nodes().field("ij")); + const size_t numNodes = field_view.shape(0); + // "ReadServer buffer size does not equal the number of horizontal nodes in the field_view" + assert(numNodes <= buffer.size()); + for (size_t inode = 0; inode < numNodes; ++inode) { + if (ghost(inode)) continue; + const int64_t ibuf = index_glbarray_(ij(inode, 0), ij(inode, 1)); + field_view(inode, z_index) = buffer[ibuf]; + } +} + +/// \brief Move vertical data from a buffer into an atlas arrayView. +/// \param buffer Vector of data to read +/// \param field_view View into the atlas field to store the data. +void ReadServer::fill_vertical_field(const std::vector& buffer, + atlas::array::ArrayView& field_view) const { + auto ghost = atlas::array::make_view(this->mesh_.nodes().ghost()); + const size_t num_nodes = field_view.shape(0); + const size_t num_levels = field_view.shape(1); + // "ReadServer buffer size does not equal the number of levels in the field_view" + assert(num_levels <= buffer.size()); + + // even for 1D depths, store the data in an atlas 3D field - inefficient but flexible + for (size_t inode = 0; inode < num_nodes; ++inode) { + for (size_t k = 0; k < num_levels; ++k) { + if (ghost(inode)) continue; + field_view(inode, k) = buffer[k]; + } + } +} + +/// \brief Read a NetCDF variable into an atlas field. +/// \param var_name The netCDF name of the variable to read. +/// \param t_index The time index for the data to read. +/// \param field_view View into the atlas field to store the data. +void ReadServer::read_var(const std::string& var_name, + const size_t t_index, + atlas::array::ArrayView& field_view) { + + size_t n_levels = field_view.shape(1); + size_t size = index_glbarray_.nx_halo_WE * index_glbarray_.ny_halo_NS; + + std::vector buffer; + // For each level + for (size_t iLev = 0; iLev < n_levels; iLev++) { + // read the data for that level onto the root processor + this->read_var_on_root(var_name, t_index, iLev, buffer); + + // mpi distribute that data out to all processors + atlas::mpi::comm().broadcast(&buffer.front(), size, mpiroot); + + // each processor fills out its field_view from the buffer + this->fill_field(buffer, iLev, field_view); + } +} + +/// \brief Read a vertical variable into an atlas field. +/// \param var_name The netCDF name of the variable to read. +/// \param field_view View into the atlas field to store the data. +void ReadServer::read_vertical_var(const std::string& var_name, + atlas::array::ArrayView& field_view) { + + size_t n_levels = field_view.shape(1); + size_t size = index_glbarray_.nx_halo_WE * index_glbarray_.ny_halo_NS; + + std::vector buffer; + + // read the data onto the root processor + this->read_vertical_var_on_root(var_name, n_levels, buffer); + + // mpi distribute that data out to all processors + atlas::mpi::comm().broadcast(&buffer.front(), size, mpiroot); + + // each processor fills out its field_view from the buffer + this->fill_vertical_field(buffer, field_view); +} + +/// \brief Find the nearest datetime index to a datetime on the MPI root only. +/// \param datetime Search for the index of the time slice in the file nearest this datetime. +/// \return The index of the nearest time slice in the file. +size_t ReadServer::get_nearest_datetime_index(const util::DateTime& datetime) const { + size_t t_index; + + if (myrank == mpiroot) { + t_index = reader_->get_nearest_datetime_index(datetime); + } + + // mpi distribute that data out to all processors + atlas::mpi::comm().broadcast(t_index, mpiroot); + + return t_index; +} + +/// \brief Read the _FillValue for a variable, defaulting to the minimum value +/// for the datatype. Read on the MPI root process only. +/// \param name Name of the netCDF variable containing the _FillValue attribute to retrieve. +/// \return The fill value for this netCDF variable. +template T ReadServer::read_fillvalue(const std::string& name) const { + T fillvalue; + + if (myrank == mpiroot) { + fillvalue = reader_->read_fillvalue(name); + } + + // mpi distribute that data out to all processors + atlas::mpi::comm().broadcast(fillvalue, mpiroot); + + return fillvalue; +} +template int ReadServer::read_fillvalue(const std::string& name) const; +template float ReadServer::read_fillvalue(const std::string& name) const; +template double ReadServer::read_fillvalue(const std::string& name) const; +} // namespace orcamodel diff --git a/src/orca-jedi/nemo_io/ReadServer.h b/src/orca-jedi/nemo_io/ReadServer.h new file mode 100644 index 0000000..0c97c97 --- /dev/null +++ b/src/orca-jedi/nemo_io/ReadServer.h @@ -0,0 +1,59 @@ +/* + * (C) British Crown Copyright 2024 Met Office + */ + +#pragma once + +#include +#include +#include + +#include "oops/util/DateTime.h" +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/array.h" +#include "atlas/mesh.h" +#include "atlas/grid.h" +#include "atlas-orca/grid/OrcaGrid.h" + +#include "orca-jedi/nemo_io/OrcaIndex.h" +#include "orca-jedi/nemo_io/NemoFieldReader.h" + +#include "eckit/exception/Exceptions.h" + +namespace orcamodel { +class ReadServer { + public: + explicit ReadServer(const eckit::PathName& file_path, + const atlas::Mesh& mesh); + ReadServer(ReadServer &&) = default; + ReadServer(const ReadServer &) = default; + ReadServer &operator=(ReadServer &&) = default; + ReadServer &operator=(const ReadServer &) = default; + void read_var(const std::string& var_name, + const size_t t_index, + atlas::array::ArrayView& field_view); + void read_vertical_var(const std::string& var_name, + atlas::array::ArrayView& field_view); + size_t get_nearest_datetime_index(const util::DateTime& datetime) const; + template T read_fillvalue(const std::string& nemo_var_name) const; + + private: + void read_var_on_root(const std::string& var_name, + const size_t t_index, + const size_t z_index, + std::vector& buffer) const; + void read_vertical_var_on_root(const std::string& var_name, + const size_t n_levels, + std::vector& buffer) const; + void fill_field(const std::vector& buffer, + const size_t z_index, + atlas::array::ArrayView& field_view) const; + void fill_vertical_field(const std::vector& buffer, + atlas::array::ArrayView& field_view) const; + const size_t mpiroot = 0; + const size_t myrank = atlas::mpi::rank(); + const atlas::Mesh& mesh_; + const OrcaIndex index_glbarray_; + std::unique_ptr reader_; +}; +} // namespace orcamodel diff --git a/src/orca-jedi/state/StateIOUtils.cc b/src/orca-jedi/state/StateIOUtils.cc index 0024257..b11abdc 100644 --- a/src/orca-jedi/state/StateIOUtils.cc +++ b/src/orca-jedi/state/StateIOUtils.cc @@ -15,7 +15,7 @@ #include "orca-jedi/geometry/Geometry.h" #include "orca-jedi/state/State.h" -#include "orca-jedi/nemo_io/NemoFieldReader.h" +#include "orca-jedi/nemo_io/ReadServer.h" #include "orca-jedi/nemo_io/NemoFieldWriter.h" #define NEMO_FILL_TOL 1e-6 @@ -43,12 +43,12 @@ void readFieldsFromFile( auto nemo_field_path = eckit::PathName(nemo_file_name); oops::Log::debug() << "orcamodel::readFieldsFromFile:: nemo_field_path " << nemo_field_path << std::endl; - NemoFieldReader nemo_file(nemo_field_path); + ReadServer nemo_reader(nemo_field_path, geom.mesh()); // Read fields from Nemo field file // field names in the atlas fieldset are assumed to match their names in // the field file - const size_t time_indx = nemo_file.get_nearest_datetime_index(valid_date); + const size_t time_indx = nemo_reader.get_nearest_datetime_index(valid_date); oops::Log::debug() << "orcamodel::readFieldsFromFile:: time_indx " << time_indx << std::endl; @@ -71,15 +71,12 @@ void readFieldsFromFile( << std::endl; if (geom.variable_in_variable_type(fieldName, variable_type)) { auto field_view = atlas::array::make_view(field); - if (varCoordTypeMap[fieldName] == "surface") { - nemo_file.read_surf_var(nemoName, geom.mesh(), time_indx, field_view); - } else if (varCoordTypeMap[fieldName] == "vertical") { - nemo_file.read_vertical_var(nemoName, geom.mesh(), field_view); + if (varCoordTypeMap[fieldName] == "vertical") { + nemo_reader.read_vertical_var(nemoName, field_view); } else { - nemo_file.read_volume_var(nemoName, geom.mesh(), time_indx, - field_view); + nemo_reader.read_var(nemoName, time_indx, field_view); } - auto missing_value = nemo_file.read_fillvalue(nemoName); + auto missing_value = nemo_reader.read_fillvalue(nemoName); field.metadata().set("missing_value", missing_value); field.metadata().set("missing_value_type", "approximately-equals"); field.metadata().set("missing_value_epsilon", NEMO_FILL_TOL); diff --git a/src/tests/orca-jedi/CMakeLists.txt b/src/tests/orca-jedi/CMakeLists.txt index 3710833..66bc0e0 100644 --- a/src/tests/orca-jedi/CMakeLists.txt +++ b/src/tests/orca-jedi/CMakeLists.txt @@ -13,8 +13,8 @@ ecbuild_add_test( TARGET test_orcajedi_nemo_io_field_writer.x CONDITION eckit_HAVE_MPI LIBS orcamodel ) -ecbuild_add_test( TARGET test_orcajedi_nemo_io_orca_field_reader.x - SOURCES test_nemo_io_orca_field_reader.cc +ecbuild_add_test( TARGET test_orcajedi_nemo_io_read_server.x + SOURCES test_nemo_io_read_server.cc ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} MPI 2 CONDITION eckit_HAVE_MPI diff --git a/src/tests/orca-jedi/test_nemo_io_field_reader.cc b/src/tests/orca-jedi/test_nemo_io_field_reader.cc index e9f3e29..5c238b0 100644 --- a/src/tests/orca-jedi/test_nemo_io_field_reader.cc +++ b/src/tests/orca-jedi/test_nemo_io_field_reader.cc @@ -93,82 +93,6 @@ CASE("test read_var_slice reads vector") { EXPECT_EQUAL(data[5], 171); } -CASE("test read_surf_var reads field array view") { - eckit::PathName test_data_path("../Data/simple_nemo.nc"); - - auto array_shape = atlas::array::ArrayShape{6, 1}; - auto array_datatype = atlas::array::DataType::real64(); - std::unique_ptr test_array( - atlas::array::Array::create(array_shape) ); - - atlas::array::ArrayView arrayView = - atlas::array::make_view(*test_array); - - NemoFieldReader field_reader(test_data_path); - field_reader.read_surf_var("iiceconc", 2, arrayView); - - EXPECT_EQUAL(arrayView(0, 0), 122); - EXPECT_EQUAL(arrayView(5, 0), 172); -} - -CASE("test read_volume_var reads field array view") { - eckit::PathName test_data_path("../Data/simple_nemo.nc"); - - int nx = 3; int ny = 2; int nz = 2; int nt = 3; - - auto array_shape = atlas::array::ArrayShape{nx*ny, nz}; - auto array_datatype = atlas::array::DataType::real64(); - std::unique_ptr test_array( - atlas::array::Array::create(array_shape) ); - - atlas::array::ArrayView arrayView = - atlas::array::make_view(*test_array); - - NemoFieldReader field_reader(test_data_path); - field_reader.read_volume_var("votemper", 1, arrayView); - - for (int i = 0; i < nx; ++i) { - for (int j = 0; j < ny; ++j) { - for (int k = 0; k < nz; ++k) { - double test_val = 2000 + (k + 1)*100 + (j + 1)*10 + i + 1; - std::cout << "arrayView(" << j << "+" << i << "*" << ny << ", " - << k << ") "; - std::cout << arrayView(i+j*nx, k) << " ~ " << test_val << std::endl; - EXPECT_EQUAL(arrayView(i+j*nx, k), test_val); - } - } - } -} - -CASE("test read_vertical_var reads field array view") { - eckit::PathName test_data_path("../Data/simple_nemo.nc"); - - int nx = 3; int ny = 2; int nz = 2; int nt = 3; - - auto array_shape = atlas::array::ArrayShape{nx*ny, nz}; - auto array_datatype = atlas::array::DataType::real64(); - std::unique_ptr test_array( - atlas::array::Array::create(array_shape) ); - - atlas::array::ArrayView arrayView = - atlas::array::make_view(*test_array); - - NemoFieldReader field_reader(test_data_path); - field_reader.read_vertical_var("nav_lev", arrayView); - - for (int i = 0; i < nx; ++i) { - for (int j = 0; j < ny; ++j) { - for (int k = 0; k < nz; ++k) { - double test_val = k*100; - std::cout << "arrayView(" << j << "+" << i << "*" << ny << ", " - << k << ") "; - std::cout << arrayView(i+j*nx, k) << " ~ " << test_val << std::endl; - EXPECT_EQUAL(arrayView(i+j*nx, k), test_val); - } - } - } -} - } // namespace test } // namespace orcamodel diff --git a/src/tests/orca-jedi/test_nemo_io_field_writer.cc b/src/tests/orca-jedi/test_nemo_io_field_writer.cc index 43c56f3..5e59475 100644 --- a/src/tests/orca-jedi/test_nemo_io_field_writer.cc +++ b/src/tests/orca-jedi/test_nemo_io_field_writer.cc @@ -99,24 +99,11 @@ CASE("test parallel serially distributed write field array views") { SECTION("depth field matches with data in memory") { if (rank == 0) { NemoFieldReader field_reader(test_data_path); - atlas::Field depth_field(funcSpace.createField( - atlas::option::name("depth") | - atlas::option::levels(3))); - auto depth_fv = atlas::array::make_view(temp_field); + auto vert_levels = field_reader.read_vertical_var("z", 3); - field_reader.read_vertical_var("z", depth_fv); - for (atlas::idx_t iNode = 0; iNode < temp_fv.shape(0); ++iNode) { - if (ghost(iNode)) continue; - EXPECT_EQUAL(depth_fv(iNode, 0), 1); - } - for (atlas::idx_t iNode = 0; iNode < temp_fv.shape(0); ++iNode) { - if (ghost(iNode)) continue; - EXPECT_EQUAL(depth_fv(iNode, 1), 2); - } - for (atlas::idx_t iNode = 0; iNode < temp_fv.shape(0); ++iNode) { - if (ghost(iNode)) continue; - EXPECT_EQUAL(depth_fv(iNode, 2), 3); - } + EXPECT_EQUAL(vert_levels[0], 1); + EXPECT_EQUAL(vert_levels[1], 2); + EXPECT_EQUAL(vert_levels[2], 3); } } diff --git a/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc b/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc deleted file mode 100644 index cf55587..0000000 --- a/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc +++ /dev/null @@ -1,234 +0,0 @@ -/* - * (C) British Crown Copyright 2024 Met Office - */ - -#include "eckit/log/Bytes.h" -#include "eckit/testing/Test.h" - -#include "oops/util/DateTime.h" - -#include "atlas/array.h" -#include "atlas/util/Config.h" -#include "atlas/functionspace/NodeColumns.h" -#include "atlas/mesh.h" -#include "atlas/grid.h" -#include "atlas/meshgenerator.h" -#include "atlas/field/Field.h" - -#include "atlas-orca/grid/OrcaGrid.h" - -#include "eckit/exception/Exceptions.h" - -#include "orca-jedi/nemo_io/NemoFieldReader.h" - -#include "tests/orca-jedi/OrcaModelTestEnvironment.h" - -namespace orcamodel { -namespace test { - -//----------------------------------------------------------------------------- - -CASE("test parallel serially distributed reads field array view") { - eckit::PathName test_data_path("../Data/orca2_t_nemo.nc"); - - atlas::OrcaGrid grid("ORCA2_T"); - auto meshgen_config = grid.meshgenerator(); - - atlas::MeshGenerator meshgen(meshgen_config); - auto partitioner_config = grid.partitioner(); - partitioner_config.set("type", "serial"); - auto partitioner = atlas::grid::Partitioner(partitioner_config); - auto mesh = meshgen.generate(grid, partitioner); - auto funcSpace = atlas::functionspace::NodeColumns(mesh); - - NemoFieldReader field_reader(test_data_path); - std::vector raw_data = field_reader.read_var_slice("iiceconc", 0, 0); - - atlas::Field field(funcSpace.createField( - atlas::option::name("iiceconc") | - atlas::option::levels(1))); - auto field_view = atlas::array::make_view(field); - - field_reader.read_surf_var("iiceconc", mesh, 0, field_view); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - for (size_t i = 0; i < field_view.size(); ++i) { - if (ghost(i)) continue; - if (raw_data[i] != field_view(i, 0)) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " 1 proc " << raw_data[i] << " " - << " with mesh " << field_view(i, 0) << std::endl; - } - EXPECT_EQUAL(raw_data[i], field_view(i, 0)); - } - - SECTION("volume field") { - atlas::Field field(funcSpace.createField( - atlas::option::name("votemper") | - atlas::option::levels(3))); - auto field_view = atlas::array::make_view(field); - - field_reader.read_volume_var("votemper", mesh, 0, field_view); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - std::vector raw_data; - for (int k =0; k <3; k++) { - raw_data = field_reader.read_var_slice("votemper", 0, k); - for (int i = 0; i < field_view.shape(0); ++i) { - if (ghost(i)) continue; - if (raw_data[i] != field_view(i, k)) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " 1 proc " << raw_data[i] << " " - << " with mesh " << field_view(i, k) << std::endl; - } - EXPECT_EQUAL(raw_data[i], field_view(i, k)); - } - } - } - - SECTION("depth field") { - atlas::Field field(funcSpace.createField( - atlas::option::name("depth") | - atlas::option::levels(3))); - auto field_view = atlas::array::make_view(field); - - field_reader.read_vertical_var("nav_lev", mesh, field_view); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - std::vector levels{0, 10, 100}; - for (int k =0; k <3; k++) { - for (int i = 0; i < field_view.shape(0); ++i) { - if (ghost(i)) continue; - if (levels[k] != field_view(i, k)) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " 1 proc " << levels[k] << " " - << " with mesh " << field_view(i, k) << std::endl; - } - EXPECT_EQUAL(levels[k], field_view(i, k)); - } - } - } -} - -// Disable this section until jopa-bundle uses a new version of atlas -CASE("test parallel domain distributed read_surf_var reads field array view") { - // NOTE: At this time, atlas-orca is only capable of domain distribution for - // ORCA1 and higher resolution - eckit::PathName test_data_path("../Data/orca1_t_nemo.nc"); - - atlas::OrcaGrid grid("ORCA1_T"); - auto meshgen_config = grid.meshgenerator(); - - atlas::MeshGenerator meshgen(meshgen_config); - auto partitioner_config = grid.partitioner(); - partitioner_config.set("type", "checkerboard"); - auto partitioner = atlas::grid::Partitioner(partitioner_config); - auto mesh = meshgen.generate(grid, partitioner); - auto funcSpace = atlas::functionspace::NodeColumns(mesh); - atlas::Field field(funcSpace.createField( - atlas::option::name("iiceconc") | - atlas::option::levels(1))); - - auto field_view = atlas::array::make_view(field); - - NemoFieldReader field_reader(test_data_path); - field_reader.read_surf_var("iiceconc", mesh, 0, field_view); - - std::vector raw_data = field_reader.read_var_slice("iiceconc", 0, 0); - - int nx_halo_WE = grid.nx() + grid.haloEast() + grid.haloWest(); - int ny_halo_NS = grid.ny() + grid.haloNorth() + grid.haloSouth(); - - int iy_glb_min = -grid.haloSouth(); - int ix_glb_min = -grid.haloWest(); - int glbarray_offset = -(nx_halo_WE * iy_glb_min) - ix_glb_min; - int glbarray_jstride = nx_halo_WE; - - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - for (size_t i = 0; i < mesh.nodes().size(); ++i) { - if (ghost(i)) continue; - int orca_i = ij(i, 0); - int orca_j = ij(i, 1); - size_t raw_idx = glbarray_offset + orca_j * glbarray_jstride + orca_i; - if (field_view(i, 0) != raw_data[raw_idx]) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " checkerboard " << field_view(i, 0) << " " - << " raw " << raw_data[raw_idx] << std::endl; - } - EXPECT_EQUAL(field_view(i, 0), raw_data[raw_idx]); - } - - SECTION("volume field") { - atlas::Field field(funcSpace.createField( - atlas::option::name("votemper") | - atlas::option::levels(3))); - auto field_view = atlas::array::make_view(field); - - field_reader.read_volume_var("votemper", mesh, 0, field_view); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - std::vector raw_data; - for (int k =0; k <3; k++) { - raw_data = field_reader.read_var_slice("votemper", 0, k); - for (int i = 0; i < field_view.shape(0); ++i) { - if (ghost(i)) continue; - int orca_i = ij(i, 0); - int orca_j = ij(i, 1); - size_t raw_idx = glbarray_offset + orca_j * glbarray_jstride + orca_i; - if (field_view(i, k) != raw_data[raw_idx]) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " checkerboard " << field_view(i, k) << " " - << " raw " << raw_data[raw_idx] << std::endl; - } - EXPECT_EQUAL(field_view(i, k), raw_data[raw_idx]); - } - } - } - SECTION("depth field") { - atlas::Field field(funcSpace.createField( - atlas::option::name("depth") | - atlas::option::levels(3))); - auto field_view = atlas::array::make_view(field); - - field_reader.read_vertical_var("z", mesh, field_view); - auto ij = atlas::array::make_view(mesh.nodes().field("ij")); - - auto ghost = atlas::array::make_view(mesh.nodes().ghost()); - std::vector levels{0, 10, 100}; - for (int k =0; k <3; k++) { - for (int i = 0; i < field_view.shape(0); ++i) { - if (ghost(i)) continue; - if (levels[k] != field_view(i, k)) { - std::cout << "mismatch: " - << " ij(" << i << ", 0) " << ij(i, 0) - << " ij(" << i << ", 1) " << ij(i, 1) - << " 1 proc " << levels[k] << " " - << " with mesh " << field_view(i, k) << std::endl; - } - EXPECT_EQUAL(levels[k], field_view(i, k)); - } - } - } -} - -} // namespace test -} // namespace orcamodel - -int main(int argc, char** argv) { - return orcamodel::test::run(argc, argv); -} diff --git a/src/tests/orca-jedi/test_nemo_io_read_server.cc b/src/tests/orca-jedi/test_nemo_io_read_server.cc new file mode 100644 index 0000000..6e6d715 --- /dev/null +++ b/src/tests/orca-jedi/test_nemo_io_read_server.cc @@ -0,0 +1,158 @@ +/* + * (C) British Crown Copyright 2024 Met Office + */ + +#include "eckit/log/Bytes.h" +#include "eckit/testing/Test.h" + +#include "oops/util/DateTime.h" + +#include "atlas/parallel/mpi/mpi.h" +#include "atlas/array.h" +#include "atlas/util/Config.h" +#include "atlas/functionspace/NodeColumns.h" +#include "atlas/mesh.h" +#include "atlas/grid.h" +#include "atlas/meshgenerator.h" +#include "atlas/field/Field.h" + +#include "atlas-orca/grid/OrcaGrid.h" + +#include "eckit/exception/Exceptions.h" + +#include "orca-jedi/nemo_io/ReadServer.h" +#include "orca-jedi/nemo_io/NemoFieldReader.h" +#include "orca-jedi/nemo_io/OrcaIndex.h" + +#include "tests/orca-jedi/OrcaModelTestEnvironment.h" + +namespace orcamodel { +namespace test { + + +//----------------------------------------------------------------------------- + +CASE("test MPI distributed reads field array view") { + eckit::PathName test_data_path("../Data/orca2_t_nemo.nc"); + + auto partitioner_names = std::vector{"serial", "checkerboard"}; + for (std::string partitioner_name : partitioner_names) { + atlas::OrcaGrid grid("ORCA2_T"); + auto meshgen_config = grid.meshgenerator(); + + atlas::MeshGenerator meshgen(meshgen_config); + auto partitioner_config = grid.partitioner(); + partitioner_config.set("type", partitioner_name); + auto partitioner = atlas::grid::Partitioner(partitioner_config); + auto mesh = meshgen.generate(grid, partitioner); + auto funcSpace = atlas::functionspace::NodeColumns(mesh); + auto orca_index = OrcaIndex(mesh.grid()); + + std::vector raw_data; + { + NemoFieldReader field_reader(test_data_path); + raw_data = field_reader.read_var_slice("iiceconc", 0, 0); + } + + atlas::Field field(funcSpace.createField( + atlas::option::name("iiceconc") | + atlas::option::levels(1))); + auto field_view = atlas::array::make_view(field); + + ReadServer read_server(test_data_path, mesh); + read_server.read_var("iiceconc", 0, field_view); + + auto ij = atlas::array::make_view(mesh.nodes().field("ij")); + + auto ghost = atlas::array::make_view(mesh.nodes().ghost()); + for (size_t i = 0; i < field_view.size(); ++i) { + if (ghost(i)) continue; + double raw_value = raw_data[orca_index(ij(i, 0), ij(i, 1))]; + if (raw_value != field_view(i, 0)) { + std::cout << "mismatch: " + << " ij(" << i << ", 0) " << ij(i, 0) + << " ij(" << i << ", 1) " << ij(i, 1) + << " 1 proc " << raw_value << " " + << " with mesh " << field_view(i, 0) << std::endl; + } + EXPECT_EQUAL(raw_value, field_view(i, 0)); + } + + SECTION(partitioner_name + " volume field") { + atlas::Field field(funcSpace.createField( + atlas::option::name("votemper") | + atlas::option::levels(3))); + auto field_view = atlas::array::make_view(field); + + NemoFieldReader field_reader(test_data_path); + read_server.read_var("votemper", 0, field_view); + + auto ij = atlas::array::make_view(mesh.nodes().field("ij")); + + auto ghost = atlas::array::make_view(mesh.nodes().ghost()); + std::vector raw_data; + for (int k =0; k <3; k++) { + raw_data = field_reader.read_var_slice("votemper", 0, k); + for (int i = 0; i < field_view.shape(0); ++i) { + double raw_value = raw_data[orca_index(ij(i, 0), ij(i, 1))]; + if (ghost(i)) continue; + if (raw_value != field_view(i, k)) { + std::cout << "mismatch: " + << " ij(" << i << ", 0) " << ij(i, 0) + << " ij(" << i << ", 1) " << ij(i, 1) + << " 1 proc " << raw_value << " " + << " with mesh " << field_view(i, k) << std::endl; + } + EXPECT_EQUAL(raw_value, field_view(i, k)); + } + } + } + + SECTION(partitioner_name + " depth field") { + atlas::Field field(funcSpace.createField( + atlas::option::name("depth") | + atlas::option::levels(3))); + auto field_view = atlas::array::make_view(field); + + read_server.read_vertical_var("nav_lev", field_view); + + auto ij = atlas::array::make_view(mesh.nodes().field("ij")); + + auto ghost = atlas::array::make_view(mesh.nodes().ghost()); + std::vector levels{0, 10, 100}; + for (int k =0; k <3; k++) { + for (int i = 0; i < field_view.shape(0); ++i) { + if (ghost(i)) continue; + if (levels[k] != field_view(i, k)) { + std::cout << "mismatch: " + << " ij(" << i << ", 0) " << ij(i, 0) + << " ij(" << i << ", 1) " << ij(i, 1) + << " 1 proc " << levels[k] << " " + << " with mesh " << field_view(i, k) << std::endl; + } + EXPECT_EQUAL(levels[k], field_view(i, k)); + } + } + } + SECTION(partitioner_name + " get nearest datetime index") { + util::DateTime test_datetime{"1970-01-01T00:00:00Z"}; + size_t index = read_server.get_nearest_datetime_index(test_datetime); + EXPECT_EQUAL(index, 0); + } + + SECTION(partitioner_name + " get field _FillValue") { + auto missing_value = read_server.read_fillvalue("iiceconc"); + EXPECT_EQUAL(missing_value, -32768.); + + auto default_missing_value = read_server.read_fillvalue("nav_lat"); + EXPECT_EQUAL(default_missing_value, std::numeric_limits::lowest()); + } + } +} + +} // namespace test +} // namespace orcamodel + +int main(int argc, char** argv) { + return orcamodel::test::run(argc, argv); +}