Prepare for setting XIOS up via nextSIM-DG (#728)

# Prepare for setting XIOS up via nextSIM-DG

Partially addresses #633.

This PR reworks the XIOS setup workflow in preparation for farming much
of it off to happen automatically via nextSIM-DG. For example, when
`ModelArray::setDimension` is called with `ModelArray::Dimension::X`, we
can have the XIOS handler set the local and global X sizes
correspondingly.

Main changes:
* Use `xy_domain` and `z_axis` consistently in XIOS read and write
demos.
* Reorder operations in XIOS read and write demos, with TODOs planning
out how the setup steps will be automated via nextSIM-DG.
* Automatically create an XIOS Grid `grid_2D` upon creating an XIOS
Domain with the special name `xy_domain`.
* Automatically create an XIOS Grid `grid_3D` upon creating XIOS Domain
and Axis with the special names `xy_domain` and `z_axis`, respectively.

Other minor changes:
* Finish moving from initialising strings with `{}` to `=`.
* Add a check to ensure local and global sizes match for `Z` dimension
because we only MPI parallelise in the horizontal.

core/src/ModelArray.cpp

@@ -1,7 +1,7 @@
 /*!
- * @file ModelData.cpp
+ * @file   ModelArray.cpp
  *
- * @date Feb 24, 2022
+ * @date   Feb 24, 2022
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -203,6 +203,9 @@ void ModelArray::setNComponents(std::map<Type, size_t> cMap)
 #ifdef USE_MPI
 void ModelArray::setDimension(Dimension dim, size_t globalLength, size_t localLength, size_t start)
 {
+    if ((dim == ModelArray::Dimension::Z) && (globalLength != localLength)) {
+        throw std::invalid_argument("Parallelism in the vertical not supported.");
+    }
     definedDimensions.at(dim).setLengths(globalLength, localLength, start);
 #else
 void ModelArray::setDimension(Dimension dim, size_t globalLength)

core/src/Xios.cpp

@@ -2,7 +2,7 @@
  * @file    Xios.cpp
  * @author  Tom Meltzer <tdm39@cam.ac.uk>
  * @author  Joe Wallwork <jw2423@cam.ac.uk>
- * @date    21 August 2024
+ * @date    01 Nov 2024
  * @brief   XIOS interface implementation
  * @details
  *
@@ -18,7 +18,6 @@
 #include <boost/date_time/posix_time/time_parsers.hpp>
 #if USE_XIOS
 
-#include "include/ModelArray.hpp"
 #include "include/Xios.hpp"
 
 #include <boost/algorithm/string.hpp>
@@ -125,7 +124,7 @@ int Xios::getClientMPIRank() { return mpi_rank; }
  */
 bool Xios::isInitialized()
 {
-    bool init { false };
+    bool init = false;
     cxios_context_is_initialized(contextId.c_str(), contextId.length(), &init);
     return init;
 }
@@ -327,7 +326,7 @@ void Xios::updateCalendar(const int stepNumber) { cxios_update_calendar(stepNumb
  */
 xios::CAxisGroup* Xios::getAxisGroup()
 {
-    const std::string groupId = { "axis_definition" };
+    const std::string groupId = "axis_definition";
     xios::CAxisGroup* group = NULL;
     cxios_axisgroup_handle_create(&group, groupId.c_str(), groupId.length());
     if (!group) {
@@ -358,7 +357,10 @@ xios::CAxis* Xios::getAxis(const std::string axisId)
 }
 
 /*!
- * Create an axis with some ID
+ * Create an axis with some ID.
+ *
+ * If the axis ID is 'z_axis' and a domain called 'xy_domain' exists then a grid called 'grid_3D'
+ * will automatically be created with this axis and that domain.
  *
  * @param the axis ID
  */
@@ -378,6 +380,15 @@ void Xios::createAxis(const std::string axisId)
     if (!exists) {
         throw std::runtime_error("Xios: Failed to create axis '" + axisId + "'");
     }
+    if (axisId == "z_axis") {
+        std::string domainId = "xy_domain";
+        cxios_domain_valid_id(&exists, domainId.c_str(), domainId.length());
+        if (exists) {
+            createGrid("grid_3D");
+            gridAddDomain("grid_3D", "xy_domain");
+            gridAddAxis("grid_3D", "z_axis");
+        }
+    }
 }
 
 /*!
@@ -467,7 +478,7 @@ std::vector<double> Xios::getAxisValues(const std::string axisId)
  */
 xios::CDomainGroup* Xios::getDomainGroup()
 {
-    const std::string groupId = { "domain_definition" };
+    const std::string groupId = "domain_definition";
     xios::CDomainGroup* group = NULL;
     cxios_domaingroup_handle_create(&group, groupId.c_str(), groupId.length());
     if (!group) {
@@ -498,7 +509,11 @@ xios::CDomain* Xios::getDomain(const std::string domainId)
 }
 
 /*!
- * Create a domain with some ID
+ * Create a domain with some ID.
+ *
+ * If the domain ID is 'xy_domain' then a grid called 'grid_2D' will automatically be created with
+ * this domain. If an axis called 'z_axis' also exists then a grid called 'grid_3D' will
+ * automatically be created with this domain and that axis.
  *
  * @param the domain ID
  */
@@ -518,6 +533,17 @@ void Xios::createDomain(const std::string domainId)
     if (!exists) {
         throw std::runtime_error("Xios: Failed to create domain '" + domainId + "'");
     }
+    if (domainId == "xy_domain") {
+        createGrid("grid_2D");
+        gridAddDomain("grid_2D", "xy_domain");
+        std::string axisId = "z_axis";
+        cxios_axis_valid_id(&exists, axisId.c_str(), axisId.length());
+        if (exists) {
+            createGrid("grid_3D");
+            gridAddDomain("grid_3D", "xy_domain");
+            gridAddAxis("grid_3D", "z_axis");
+        }
+    }
 }
 
 /*!
@@ -862,7 +888,7 @@ std::vector<double> Xios::getDomainLocalYValues(const std::string domainId)
  */
 xios::CGridGroup* Xios::getGridGroup()
 {
-    const std::string groupId = { "grid_definition" };
+    const std::string groupId = "grid_definition";
     xios::CGridGroup* group = NULL;
     cxios_gridgroup_handle_create(&group, groupId.c_str(), groupId.length());
     if (!group) {
@@ -1003,7 +1029,7 @@ std::vector<std::string> Xios::gridGetDomainIds(const std::string gridId)
  */
 xios::CFieldGroup* Xios::getFieldGroup()
 {
-    const std::string groupId = { "field_definition" };
+    const std::string groupId = "field_definition";
     xios::CFieldGroup* group = NULL;
     cxios_fieldgroup_handle_create(&group, groupId.c_str(), groupId.length());
     if (!group) {
@@ -1241,7 +1267,7 @@ Duration Xios::getFieldFreqOffset(const std::string fieldId)
  */
 xios::CFileGroup* Xios::getFileGroup()
 {
-    const std::string groupId = { "file_definition" };
+    const std::string groupId = "file_definition";
     xios::CFileGroup* group = NULL;
     cxios_filegroup_handle_create(&group, groupId.c_str(), groupId.length());
     if (!group) {

core/src/include/ModelArray.hpp

@@ -1,7 +1,7 @@
 /*!
- * @file ModelData.hpp
+ * @file   ModelArray.hpp
  *
- * @date Feb 24, 2022
+ * @date   31 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 

core/test/XiosGrid_test.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file    XiosGrid_test.cpp
  * @author  Joe Wallwork <jw2423@cam.ac.uk>
- * @date    5 August 2024
+ * @date    01 Nov 2024
  * @brief   Tests for XIOS axes
  * @details
  * This test is designed to test axis functionality of the C++ interface
@@ -48,18 +48,18 @@ MPI_TEST_CASE("TestXiosGrid", 2)
     xios_handler.setCalendarTimestep(Duration("P0-0T01:30:00"));
 
     // Create a 4x2 horizontal domain with a partition halving the x-extent
-    xios_handler.createDomain("xy_domain");
-    xios_handler.setDomainType("xy_domain", "rectilinear");
-    xios_handler.setDomainGlobalXSize("xy_domain", 4);
-    xios_handler.setDomainGlobalYSize("xy_domain", 2);
-    xios_handler.setDomainLocalXStart("xy_domain", 2 * rank);
-    xios_handler.setDomainLocalYStart("xy_domain", 0);
-    xios_handler.setDomainLocalXValues("xy_domain", { -1.0 + rank, -0.5 + rank });
-    xios_handler.setDomainLocalYValues("xy_domain", { -1.0, 1.0 });
+    xios_handler.createDomain("domain_XY");
+    xios_handler.setDomainType("domain_XY", "rectilinear");
+    xios_handler.setDomainGlobalXSize("domain_XY", 4);
+    xios_handler.setDomainGlobalYSize("domain_XY", 2);
+    xios_handler.setDomainLocalXStart("domain_XY", 2 * rank);
+    xios_handler.setDomainLocalYStart("domain_XY", 0);
+    xios_handler.setDomainLocalXValues("domain_XY", { -1.0 + rank, -0.5 + rank });
+    xios_handler.setDomainLocalYValues("domain_XY", { -1.0, 1.0 });
 
     // Create a vertical axis with 2 points
-    xios_handler.createAxis("z_axis");
-    xios_handler.setAxisValues("z_axis", { 0.0, 1.0 });
+    xios_handler.createAxis("axis_Z");
+    xios_handler.setAxisValues("axis_Z", { 0.0, 1.0 });
 
     // --- Tests for grid API
     const std::string gridId = { "grid_2D" };
@@ -73,15 +73,15 @@ MPI_TEST_CASE("TestXiosGrid", 2)
     xios_handler.setGridName(gridId, gridName);
     REQUIRE(xios_handler.getGridName(gridId) == gridName);
     // Add axis
-    xios_handler.gridAddAxis("grid_2D", "z_axis");
+    xios_handler.gridAddAxis("grid_2D", "axis_Z");
     std::vector<std::string> axisIds = xios_handler.gridGetAxisIds(gridId);
     REQUIRE(axisIds.size() == 1);
-    REQUIRE(axisIds[0] == "z_axis");
+    REQUIRE(axisIds[0] == "axis_Z");
     // Add domain
-    xios_handler.gridAddDomain("grid_2D", "xy_domain");
+    xios_handler.gridAddDomain("grid_2D", "domain_XY");
     std::vector<std::string> domainIds = xios_handler.gridGetDomainIds(gridId);
     REQUIRE(domainIds.size() == 1);
-    REQUIRE(domainIds[0] == "xy_domain");
+    REQUIRE(domainIds[0] == "domain_XY");
 
     xios_handler.close_context_definition();
     xios_handler.context_finalize();

core/test/XiosRead_test.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file    XiosRead_test.cpp
  * @author  Joe Wallwork <jw2423@cam.ac.uk
- * @date    21 August 2024
+ * @date    01 Nov 2024
  * @brief   Tests for XIOS read method
  * @details
  * This test is designed to test the read method of the C++ interface
@@ -14,6 +14,7 @@
 #include "StructureModule/include/ParametricGrid.hpp"
 #include "include/Configurator.hpp"
 #include "include/NextsimModule.hpp"
+#include "include/ParaGridIO.hpp"
 #include "include/Xios.hpp"
 
 #include <filesystem>
@@ -40,7 +41,14 @@ MPI_TEST_CASE("TestXiosRead", 2)
     std::unique_ptr<std::istream> pcstream(new std::stringstream(config.str()));
     Configurator::addStream(std::move(pcstream));
 
+    // Create ParametricGrid and ParaGridIO instances
+    Module::setImplementation<IStructure>("Nextsim::ParametricGrid");
+    ParametricGrid grid;
+    ParaGridIO* pio = new ParaGridIO(grid);
+    grid.setIO(pio);
+
     // Initialize an Xios instance called xios_handler
+    // TODO: Create XIOS handler along with ParaGridIO instance
     Xios xios_handler;
     REQUIRE(xios_handler.isInitialized());
     const size_t size = xios_handler.getClientMPISize();
@@ -53,54 +61,50 @@ MPI_TEST_CASE("TestXiosRead", 2)
     Duration timestep("P0-0T01:30:00");
     xios_handler.setCalendarTimestep(timestep);
 
-    // Axis setup
-    const int n1 = 2;
-    const int n2 = 3;
-    const int n3 = 4;
-    const int n4 = 5;
-    xios_handler.createAxis("axis_A");
-    xios_handler.setAxisSize("axis_A", n1);
-    xios_handler.setAxisValues("axis_A", { 0, 1 });
-    xios_handler.createAxis("axis_B");
-    xios_handler.setAxisSize("axis_B", n2);
-    xios_handler.setAxisValues("axis_B", { 0, 1, 2 });
-    xios_handler.createAxis("axis_C");
-    xios_handler.setAxisSize("axis_C", n3);
-    xios_handler.setAxisValues("axis_C", { 0, 1, 2, 3 });
-    xios_handler.createAxis("axis_D");
-    xios_handler.setAxisSize("axis_D", n4);
-    xios_handler.setAxisValues("axis_D", { 0, 1, 2, 3, 4 });
-
-    // Grid setup
-    xios_handler.createGrid("grid_2D");
-    xios_handler.gridAddAxis("grid_2D", "axis_A");
-    xios_handler.gridAddAxis("grid_2D", "axis_B");
-    xios_handler.createGrid("grid_3D");
-    xios_handler.gridAddAxis("grid_3D", "axis_A");
-    xios_handler.gridAddAxis("grid_3D", "axis_B");
-    xios_handler.gridAddAxis("grid_3D", "axis_C");
-    xios_handler.createGrid("grid_4D");
-    xios_handler.gridAddAxis("grid_4D", "axis_A");
-    xios_handler.gridAddAxis("grid_4D", "axis_B");
-    xios_handler.gridAddAxis("grid_4D", "axis_C");
-    xios_handler.gridAddAxis("grid_4D", "axis_D");
+    // Set ModelArray dimensions
+    const size_t nx_glo = 4;
+    const size_t ny_glo = 2;
+    const size_t nx = 2;
+    const size_t ny = 2;
+    const size_t nz = 2;
+    ModelArray::setDimension(ModelArray::Dimension::X, nx_glo, nx, 0);
+    ModelArray::setDimension(ModelArray::Dimension::Y, ny_glo, ny, 0);
+    ModelArray::setDimension(ModelArray::Dimension::Z, nz, nz, 0);
+
+    // Create a 4x2 horizontal domain with a partition halving the x-extent
+    // TODO: Set local and global domain sizes upon calling ModelArray::setDimension for X and Y
+    xios_handler.createDomain("xy_domain");
+    xios_handler.setDomainType("xy_domain", "rectilinear");
+    xios_handler.setDomainGlobalXSize("xy_domain", nx_glo);
+    xios_handler.setDomainGlobalYSize("xy_domain", ny_glo);
+    xios_handler.setDomainLocalXStart("xy_domain", 2 * rank);
+    xios_handler.setDomainLocalYStart("xy_domain", 0);
+    xios_handler.setDomainLocalXValues("xy_domain", { -1.0 + rank, -0.5 + rank });
+    xios_handler.setDomainLocalYValues("xy_domain", { -1.0, 1.0 });
+
+    // Create a vertical axis with 2 points
+    // TODO: Set axis size upon calling ModelArray::setDimension for Z
+    xios_handler.createAxis("z_axis");
+    xios_handler.setAxisValues("z_axis", { 0.0, 1.0 });
+
+    // Create some fake data to test writing methods
+    HField field_2D(ModelArray::Type::H);
+    field_2D.resize();
+    HField field_3D(ModelArray::Type::Z);
+    field_3D.resize();
 
     // Field setup
+    // TODO: Create field along with HField
     xios_handler.createField("field_2D");
     xios_handler.setFieldOperation("field_2D", "instant");
-    xios_handler.setFieldGridRef("field_2D", "grid_2D");
+    xios_handler.setFieldGridRef("field_2D", "grid_2D"); // NOTE: grid_2D auto-generated
     xios_handler.setFieldReadAccess("field_2D", true);
     xios_handler.setFieldFreqOffset("field_2D", timestep);
     xios_handler.createField("field_3D");
     xios_handler.setFieldOperation("field_3D", "instant");
-    xios_handler.setFieldGridRef("field_3D", "grid_3D");
+    xios_handler.setFieldGridRef("field_3D", "grid_3D"); // NOTE: grid_3D auto-generated
     xios_handler.setFieldReadAccess("field_3D", true);
     xios_handler.setFieldFreqOffset("field_3D", timestep);
-    xios_handler.createField("field_4D");
-    xios_handler.setFieldOperation("field_4D", "instant");
-    xios_handler.setFieldGridRef("field_4D", "grid_4D");
-    xios_handler.setFieldReadAccess("field_4D", true);
-    xios_handler.setFieldFreqOffset("field_4D", timestep);
 
     // File setup
     xios_handler.createFile("xios_test_input");
@@ -110,21 +114,10 @@ MPI_TEST_CASE("TestXiosRead", 2)
     xios_handler.setFileParAccess("xios_test_input", "collective");
     xios_handler.fileAddField("xios_test_input", "field_2D");
     xios_handler.fileAddField("xios_test_input", "field_3D");
-    xios_handler.fileAddField("xios_test_input", "field_4D");
 
     xios_handler.close_context_definition();
 
     // --- Tests for reading to file
-    Module::setImplementation<IStructure>("Nextsim::ParametricGrid");
-    ModelArray::setDimension(ModelArray::Dimension::X, n1, n1, 0);
-    ModelArray::setDimension(ModelArray::Dimension::Y, n2, n2, 0);
-    ModelArray::setDimension(ModelArray::Dimension::Z, n3, n3, 0);
-    // Create some fake data to test writing methods
-    HField field_2D(ModelArray::Type::H);
-    field_2D.resize();
-    HField field_3D(ModelArray::Type::Z);
-    field_3D.resize();
-    // TODO: Implement 4D case
     // Verify calendar step is starting from zero
     REQUIRE(xios_handler.getCalendarStep() == 0);
     // Check the input file exists
@@ -136,24 +129,24 @@ MPI_TEST_CASE("TestXiosRead", 2)
         // Receive data from XIOS that is read from disk
         xios_handler.read("field_2D", field_2D);
         xios_handler.read("field_3D", field_3D);
-        // TODO: Implement 4D case
         // Verify timestep
         REQUIRE(xios_handler.getCalendarStep() == ts);
     }
     // Verify fields have been read in correctly
-    for (size_t j = 0; j < n2; ++j) {
-        for (size_t i = 0; i < n1; ++i) {
-            REQUIRE(field_2D(i, j) == doctest::Approx(1.0 * (i + n1 * j)));
+    for (size_t j = 0; j < ny; ++j) {
+        for (size_t i = 0; i < nx; ++i) {
+            REQUIRE(field_2D(i, j) == doctest::Approx(1.0 * (i + nx * j)));
         }
     }
-    for (size_t k = 0; k < n3; ++k) {
-        for (size_t j = 0; j < n2; ++j) {
-            for (size_t i = 0; i < n1; ++i) {
-                REQUIRE(field_3D(i, j, k) == doctest::Approx(1.0 * (i + n1 * (j + n2 * k))));
+    for (size_t k = 0; k < nz; ++k) {
+        for (size_t j = 0; j < ny; ++j) {
+            for (size_t i = 0; i < nx; ++i) {
+                REQUIRE(field_3D(i, j, k) == doctest::Approx(1.0 * (i + nx * (j + ny * k))));
             }
         }
     }
-
     xios_handler.context_finalize();
+
+    // TODO: Consider adding a 4D test case
 }
 }