diff --git a/src/orca-jedi/geometry/Geometry.cc b/src/orca-jedi/geometry/Geometry.cc index cde12bd..bbd7d83 100644 --- a/src/orca-jedi/geometry/Geometry.cc +++ b/src/orca-jedi/geometry/Geometry.cc @@ -47,14 +47,14 @@ Geometry::Geometry(const eckit::Configuration & config, grid_(config.getString("grid name")) { params_.validateAndDeserialize(config); - int64_t halo = params_.sourceMeshHalo.value().value_or(0); + int64_t halo = params_.sourceMeshHalo.value(); auto meshgen_config = grid_.meshgenerator() | atlas::option::halo(halo); atlas::MeshGenerator meshgen(meshgen_config); auto partitioner_config = grid_.partitioner(); partitioner_config.set("type", - params_.partitioner.value().value_or("serial")); + params_.partitioner.value()); partitioner_ = atlas::grid::Partitioner(partitioner_config); mesh_ = meshgen.generate(grid_, partitioner_); funcSpace_ = atlas::functionspace::NodeColumns( diff --git a/src/orca-jedi/geometry/Geometry.h b/src/orca-jedi/geometry/Geometry.h index b7d1dae..33f2e70 100644 --- a/src/orca-jedi/geometry/Geometry.h +++ b/src/orca-jedi/geometry/Geometry.h @@ -64,7 +64,7 @@ class Geometry : public util::Printable { std::string variable_type) const; bool levelsAreTopDown() const {return true;} std::string distributionType() const { - return params_.partitioner.value().value_or("serial");} + return params_.partitioner.value();} private: void print(std::ostream &) const; diff --git a/src/orca-jedi/geometry/GeometryParameters.h b/src/orca-jedi/geometry/GeometryParameters.h index f65281f..d837618 100644 --- a/src/orca-jedi/geometry/GeometryParameters.h +++ b/src/orca-jedi/geometry/GeometryParameters.h @@ -36,8 +36,17 @@ class OrcaGeometryParameters : public oops::Parameters { oops::RequiredParameter gridName {"grid name", this}; oops::RequiredParameter nLevels {"number levels", this}; - oops::OptionalParameter sourceMeshHalo {"source mesh halo", this}; - oops::OptionalParameter partitioner {"partitioner", this}; + oops::Parameter sourceMeshHalo {"source mesh halo", + "Size of the MPI halo when using a domain-distributed geometry." + " The default is 0 (no MPI halo)", + 0, + this}; + oops::Parameter partitioner { + "partitioner", + "Name of the atlas partitioner to use to MPI distribute the model data" + " The default will not distribute the data ('serial').", + "serial", + this}; }; } // namespace orcamodel diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc index 7a4a785..0533d73 100644 --- a/src/orca-jedi/state/State.cc +++ b/src/orca-jedi/state/State.cc @@ -235,9 +235,11 @@ double State::norm(const std::string & field_name) const { bool has_mv = static_cast(mv); double squares_TP = 0; size_t valid_points_TP = 0; - atlas_omp_for(atlas::idx_t j = 0; j < field_view.shape(0); ++j) { + atlas::idx_t num_h_locs = field_view.shape(0); + atlas::idx_t num_levels = field_view.shape(1); + atlas_omp_for(atlas::idx_t j = 0; j < num_h_locs; ++j) { if (!ghost(j)) { - for (atlas::idx_t k = 0; k < field_view.shape(1); ++k) { + for (atlas::idx_t k = 0; k < num_levels; ++k) { double pointValue = field_view(j, k); if (!has_mv || (has_mv && !mv(pointValue))) { squares_TP += pointValue*pointValue; @@ -251,12 +253,28 @@ double State::norm(const std::string & field_name) const { valid_points += valid_points_TP; } } - // prevent divide by zero when there are no valid model points on the - // partition - if (!valid_points) - return 0; - return sqrt(squares)/valid_points; + // serial distributions have the entire model grid on each MPI rank + if (geom_->distributionType() == "serial") { + double local_norm = 0; + // prevent divide by zero when there are no valid model points on this + // MPI rank + if (valid_points) { + local_norm = sqrt(squares)/valid_points; + } + return local_norm; + } + + + // Accumulate values across MPI ranks. + geom_->getComm().allReduceInPlace(squares, eckit::mpi::sum()); + geom_->getComm().allReduceInPlace(valid_points, eckit::mpi::sum()); + + if (valid_points) { + return sqrt(squares)/valid_points; + } + + return 0; } } // namespace orcamodel diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index f85f206..a2899de 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -28,9 +28,19 @@ ecbuild_add_test( TARGET test_orcamodel_hofx_sst ARGS testinput/hofx_nc_sst.yaml COMMAND orcamodel_hofx.x ) -ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel +ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel_checkerboard OMP 1 MPI 2 + ARGS testinput/hofx_nc_ssh_checkerboard.yaml + COMMAND orcamodel_hofx.x ) + +ecbuild_add_test( TARGET test_orcamodel_hofx_ssh_parallel_serial + OMP 1 + MPI 2 + ARGS testinput/hofx_nc_ssh.yaml + COMMAND orcamodel_hofx.x ) + +ecbuild_add_test( TARGET test_orcamodel_hofx_ssh ARGS testinput/hofx_nc_ssh.yaml COMMAND orcamodel_hofx.x ) diff --git a/src/tests/Data/orca1_t_nemo.nc b/src/tests/Data/orca1_t_nemo.nc index 26a43b3..5a4c5cf 100644 --- a/src/tests/Data/orca1_t_nemo.nc +++ b/src/tests/Data/orca1_t_nemo.nc @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2a3b664a7cc89ee768e0ec9101a25f2439a02df180a369de14850e18aa5d006b -size 4253091 +oid sha256:daa2b06c99148f8f4e08b3bdc7de763984beb7ad56700406446048c63291e180 +size 4253928 diff --git a/src/tests/orca-jedi/test_interpolator_parallel.cc b/src/tests/orca-jedi/test_interpolator_parallel.cc new file mode 100644 index 0000000..687078d --- /dev/null +++ b/src/tests/orca-jedi/test_interpolator_parallel.cc @@ -0,0 +1,333 @@ +/* + * (C) British Crown Copyright 2023 Met Office + */ + +#include +#include + +#include "eckit/log/Bytes.h" +#include "eckit/config/LocalConfiguration.h" +#include "eckit/mpi/Comm.h" +#include "eckit/testing/Test.h" + +#include "oops/base/Variables.h" +#include "oops/util/DateTime.h" +#include "oops/util/missingValues.h" + +#include "atlas/library/Library.h" +#include "atlas/util/function/VortexRollup.h" + +#include "orca-jedi/geometry/Geometry.h" +#include "orca-jedi/state/State.h" +#include "orca-jedi/state/StateParameters.h" +#include "orca-jedi/interpolator/Interpolator.h" +#include "orca-jedi/interpolator/InterpolatorParameters.h" + +#include "tests/orca-jedi/OrcaModelTestEnvironment.h" + +namespace orcamodel { +namespace test { + +const double ATOL = 1e-6; +const double CHECKERBOARD_MED_RANK = 1; + +//----------------------------------------------------------------------------- + +CASE("test serial interpolator") { + int nlevs = 3; + eckit::LocalConfiguration config; + config.set("grid name", "eORCA12_T"); + config.set("number levels", nlevs); + config.set("partitioner", "serial"); + + std::vector nemo_var_mappings(4); + nemo_var_mappings[0].set("name", "sea_surface_height_anomaly") + .set("nemo field name", "sossheig") + .set("model space", "surface"); + nemo_var_mappings[1].set("name", "sea_surface_height_anomaly_error") + .set("nemo field name", "ssh_tot_std") + .set("model space", "surface") + .set("variable type", "background variance"); + nemo_var_mappings[2].set("name", "sea_surface_foundation_temperature") + .set("nemo field name", "votemper") + .set("model space", "surface"); + nemo_var_mappings[3].set("name", "sea_water_potential_temperature") + .set("nemo field name", "votemper") + .set("model space", "volume"); + config.set("nemo variables", nemo_var_mappings); + Geometry geometry(config, eckit::mpi::comm()); + + eckit::LocalConfiguration interp_conf; + interp_conf.set("type", "unstructured-bilinear-lonlat"); + interp_conf.set("non_linear", "missing-if-all-missing"); + eckit::LocalConfiguration interpolator_conf; + interpolator_conf.set("atlas-interpolator", interp_conf); + + OrcaInterpolatorParameters params; + params.validateAndDeserialize(interpolator_conf); + + int nlocs = 0; + std::vector lons, lats; + // just off baja, stick a bunch of points as these look bad in orca12 decomp + int nlats = 10; + int nlons = 10; + // -114 to -94 + double lonStart = -114; + double lonEnd = -94; + // 10 to 15 + double latStart = 10; + double latEnd = 15; + nlocs = nlats*nlons; + lons = std::vector(); + lats = std::vector(); + for (int iLon = 0; iLon < nlons; ++iLon) { + double lon = lonStart + std::floor((lonEnd - lonStart)*iLon/nlons); + for (int iLat = 0; iLat < nlats; ++iLat) { + double lat = latStart + std::floor((latEnd - latStart)*iLat/nlats); + lons.emplace_back(lon); + lats.emplace_back(lat); + } + } + + auto rollup_plus = [](const double lon, const double lat) { + return 1 + atlas::util::function::vortex_rollup(lon, lat, 0.0); + }; + + // create a state from the test data + eckit::LocalConfiguration state_config; + std::vector state_variables { + "sea_surface_height_anomaly", + "sea_surface_foundation_temperature", + "sea_water_potential_temperature"}; + state_config.set("state variables", state_variables); + state_config.set("date", "2021-06-30T00:00:00Z"); + state_config.set("analytic initialisation", true); + state_config.set("nemo field file", "../Data/orca2_t_nemo.nc"); + state_config.set("variance field file", "../Data/orca2_t_bkg_var.nc"); + OrcaStateParameters stateParams; + stateParams.validateAndDeserialize(state_config); + State state(geometry, stateParams); + + SECTION("test interpolator succeeds even with no locations") { + Interpolator interpolator(interpolator_conf, geometry, {}, {}); + } + + Interpolator interpolator(interpolator_conf, geometry, lats, lons); + + SECTION("test interpolator.apply fails missing variable") { + oops::Variables variables({"NOTAVARIABLE"}); + std::vector vals(nlocs); + std::vector mask(nlocs); + EXPECT_THROWS_AS(interpolator.apply(variables, state, mask, vals), + eckit::BadParameter); + } + + SECTION("test interpolator.apply") { + // two variables at n locations + std::vector vals(2*nlocs); + std::vector mask(2*nlocs); + interpolator.apply(oops::Variables({"sea_surface_height_anomaly", + "sea_surface_foundation_temperature"}), state, mask, vals); + + double missing_value = util::missingValue(); + std::vector testvals(2*nlocs, 0); + for (int i = 0; i < nlocs; ++i) { + testvals[i] = rollup_plus(lons[i], lats[i]); + testvals[i+nlocs] = rollup_plus(lons[i], lats[i]); + } + + int errCount = 0; + for (size_t i=0; i < testvals.size(); ++i) { + if (std::abs(vals[i] - testvals[i]) > ATOL) { + std::cout << "[" << eckit::mpi::comm().rank() << "] lons [" << i << "] " << lons[i] + << " lats[" << i << "] " << lats[i] + << " vals[" << i << "] " << std::setprecision(12) << vals[i] + << " testvals[" << i << "] " << testvals[i] << std::endl; + ++errCount; + } + EXPECT(std::abs(vals[i] - testvals[i]) < ATOL); + } + std::cout << "[" << eckit::mpi::comm().rank() << "] errCount: " + << errCount << "/" << 2*nlocs << std::endl; + } + // SECTION("test interpolator.apply multiple levels") { + // std::vector vals(nlevs*nlocs); + // std::vector mask(nlevs*nlocs); + // interpolator.apply(oops::Variables({"sea_water_potential_temperature"}), + // state, mask, vals); + + // double missing_value = util::missingValue(); + // std::vector testvals(3*nlocs,0); + // if (eckit::mpi::comm().rank() == CHECKERBOARD_MED_RANK) { + // //testvals = std::vector({0.69885330268, 0.00010463659, 1.38328833133, + // // 0.69885330268, 0.00010463659, 1.38328833133, + // // 0.69885330268, 0.00010463659, 1.38328833133}); + // for (int i = 0; i < nlocs; ++i) { + // testvals[i] = rollup_plus(lons[i], lats[i]); + // testvals[i+nlocs] = rollup_plus(lons[i], lats[i]); + // testvals[i+2*nlocs] = rollup_plus(lons[i], lats[i]); + // } + // } + + // for (size_t i=0; i < testvals.size(); ++i) { + // //std::cout << "[" << eckit::mpi::comm().rank() << "] lons [" << i << "] " << lons[i] + // // << " lats[" << i << "] " << lats[i] + // // << " vals[" << i << "] " << std::setprecision(12) << vals[i] + // // << " testvals[" << i << "] " << testvals[i] << std::endl; + // EXPECT(std::abs(vals[i] - testvals[i]) < ATOL); + // } + // } +} + +CASE("test checkerboard interpolator") { + int nlevs = 3; + eckit::LocalConfiguration config; + config.set("grid name", "eORCA12_T"); + config.set("number levels", nlevs); + config.set("source mesh halo", 0); + config.set("partitioner", "checkerboard"); + + std::vector nemo_var_mappings(4); + nemo_var_mappings[0].set("name", "sea_surface_height_anomaly") + .set("nemo field name", "sossheig") + .set("model space", "surface"); + nemo_var_mappings[1].set("name", "sea_surface_height_anomaly_error") + .set("nemo field name", "ssh_tot_std") + .set("model space", "surface") + .set("variable type", "background variance"); + nemo_var_mappings[2].set("name", "sea_surface_foundation_temperature") + .set("nemo field name", "votemper") + .set("model space", "surface"); + nemo_var_mappings[3].set("name", "sea_water_potential_temperature") + .set("nemo field name", "votemper") + .set("model space", "volume"); + config.set("nemo variables", nemo_var_mappings); + Geometry geometry(config, eckit::mpi::comm()); + + eckit::LocalConfiguration interp_conf; + interp_conf.set("type", "unstructured-bilinear-lonlat"); + interp_conf.set("non_linear", "missing-if-all-missing"); + eckit::LocalConfiguration interpolator_conf; + interpolator_conf.set("atlas-interpolator", interp_conf); + + OrcaInterpolatorParameters params; + params.validateAndDeserialize(interpolator_conf); + + int nlocs = 0; + std::vector lons, lats; + int nlats = 10; + int nlons = 10; + // -114 to -94 + double lonStart = -114; + double lonEnd = -94; + // 10 to 15 + double latStart = 10; + double latEnd = 15; + nlocs = nlats*nlons; + lons = std::vector(); + lats = std::vector(); + for (int iLon = 0; iLon < nlons; ++iLon) { + double lon = lonStart + std::floor((lonEnd - lonStart)*iLon/nlons); + for (int iLat = 0; iLat < nlats; ++iLat) { + double lat = latStart + std::floor((latEnd - latStart)*iLat/nlats); + lons.emplace_back(lon); + lats.emplace_back(lat); + } + } + + auto rollup_plus = [](const double lon, const double lat) { + return 1 + atlas::util::function::vortex_rollup(lon, lat, 0.0); + }; + + // create a state from the test data + eckit::LocalConfiguration state_config; + std::vector state_variables { + "sea_surface_height_anomaly", + "sea_surface_foundation_temperature", + "sea_water_potential_temperature"}; + state_config.set("state variables", state_variables); + state_config.set("date", "2021-06-30T00:00:00Z"); + state_config.set("analytic initialisation", true); + state_config.set("nemo field file", "../Data/orca2_t_nemo.nc"); + state_config.set("variance field file", "../Data/orca2_t_bkg_var.nc"); + OrcaStateParameters stateParams; + stateParams.validateAndDeserialize(state_config); + State state(geometry, stateParams); + + SECTION("test interpolator succeeds even with no locations") { + Interpolator interpolator(interpolator_conf, geometry, {}, {}); + } + + Interpolator interpolator(interpolator_conf, geometry, lats, lons); + + SECTION("test interpolator.apply fails missing variable") { + oops::Variables variables({"NOTAVARIABLE"}); + std::vector vals(nlocs); + std::vector mask(nlocs); + EXPECT_THROWS_AS(interpolator.apply(variables, state, mask, vals), + eckit::BadParameter); + } + + SECTION("test interpolator.apply") { + // two variables at n locations + std::vector vals(2*nlocs); + std::vector mask(2*nlocs); + interpolator.apply(oops::Variables({"sea_surface_height_anomaly", + "sea_surface_foundation_temperature"}), state, mask, vals); + + double missing_value = util::missingValue(); + std::vector testvals(2*nlocs, 0); + for (int i = 0; i < nlocs; ++i) { + testvals[i] = rollup_plus(lons[i], lats[i]); + testvals[i+nlocs] = rollup_plus(lons[i], lats[i]); + } + + int errCount = 0; + for (size_t i=0; i < testvals.size(); ++i) { + if (std::abs(vals[i] - testvals[i]) > ATOL) { + std::cout << "[" << eckit::mpi::comm().rank() << "] lons [" << i << "] " << lons[i] + << " lats[" << i << "] " << lats[i] + << " vals[" << i << "] " << std::setprecision(12) << vals[i] + << " testvals[" << i << "] " << testvals[i] << std::endl; + ++errCount; + } + EXPECT(std::abs(vals[i] - testvals[i]) < ATOL); + } + std::cout << "[" << eckit::mpi::comm().rank() << "] errCount: " + << errCount << "/" << 2*nlocs << std::endl; + } + // SECTION("test interpolator.apply multiple levels") { + // std::vector vals(nlevs*nlocs); + // std::vector mask(nlevs*nlocs); + // interpolator.apply(oops::Variables({"sea_water_potential_temperature"}), + // state, mask, vals); + + // double missing_value = util::missingValue(); + // std::vector testvals(3*nlocs,0); + // if (eckit::mpi::comm().rank() == CHECKERBOARD_MED_RANK) { + // //testvals = std::vector({0.69885330268, 0.00010463659, 1.38328833133, + // // 0.69885330268, 0.00010463659, 1.38328833133, + // // 0.69885330268, 0.00010463659, 1.38328833133}); + // for (int i = 0; i < nlocs; ++i) { + // testvals[i] = rollup_plus(lons[i], lats[i]); + // testvals[i+nlocs] = rollup_plus(lons[i], lats[i]); + // testvals[i+2*nlocs] = rollup_plus(lons[i], lats[i]); + // } + // } + + // for (size_t i=0; i < testvals.size(); ++i) { + // //std::cout << "[" << eckit::mpi::comm().rank() << "] lons [" << i << "] " << lons[i] + // // << " lats[" << i << "] " << lats[i] + // // << " vals[" << i << "] " << std::setprecision(12) << vals[i] + // // << " testvals[" << i << "] " << testvals[i] << std::endl; + // EXPECT(std::abs(vals[i] - testvals[i]) < ATOL); + // } + // } +} + +} // 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_orca_field_reader.cc b/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc index 3a2df5e..cf55587 100644 --- a/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc +++ b/src/tests/orca-jedi/test_nemo_io_orca_field_reader.cc @@ -120,7 +120,6 @@ CASE("test parallel serially distributed reads field array view") { } // 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 @@ -151,12 +150,12 @@ CASE("test parallel domain distributed read_surf_var reads field array view") { 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_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 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); - field_reader.read_vertical_var("nav_lev", mesh, field_view); + 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}; + std::vector levels{0, 10, 100}; for (int k =0; k <3; k++) { - for (int i = 0; i