diff --git a/core/src/ParaGridIO_Xios.cpp b/core/src/ParaGridIO_Xios.cpp index dffd08878..7f8b63292 100644 --- a/core/src/ParaGridIO_Xios.cpp +++ b/core/src/ParaGridIO_Xios.cpp @@ -70,9 +70,13 @@ ModelState ParaGridIO::getModelState(const std::string& filePath) DGField field(ModelArray::Type::DG); field.resize(); state.merge(ModelState { { { fieldId, field } }, {} }); + } else if (type == ModelArray::Type::DGSTRESS) { + DGSField field(ModelArray::Type::DGSTRESS); + field.resize(); + state.merge(ModelState { { { fieldId, field } }, {} }); } else { - throw std::runtime_error( - "ParaGridIO::getModelState: field type for field" + fieldId + " is not supported."); + throw std::runtime_error("ParaGridIO::getModelState: field type for field " + fieldId + + " is not supported."); } } diff --git a/core/src/Xios.cpp b/core/src/Xios.cpp index f33ea7716..f99c8e21a 100644 --- a/core/src/Xios.cpp +++ b/core/src/Xios.cpp @@ -781,21 +781,27 @@ void Xios::affixModelMetadata() size_t counter = 0; for (ModelArray::Dimension dim : ModelArray::typeDimensions[type]) { const std::string domainName = ModelArray::definedDimensions[dim].name; - // Add 1 in the case of vertex-based field types - ModelArray::Base base = ModelArray::baseTypes[type]; - size_t extension = (base == ModelArray::Base::Vertex) ? 1 : 0; if (counter == 0) { - cxios_set_domain_ni_glo(domain, (int)metadata.getGlobalExtentX() + extension); + if (dim == ModelArray::Dimension::X) { + cxios_set_domain_ni_glo(domain, metadata.getGlobalExtentX()); + cxios_set_domain_ni(domain, metadata.getLocalExtentX()); + } else if (dim == ModelArray::Dimension::XVERTEX) { + cxios_set_domain_ni_glo(domain, metadata.getGlobalExtentX() + 1); + cxios_set_domain_ni(domain, metadata.getLocalExtentX() + 1); + } else { + throw std::runtime_error( + "Xios: Could not set domain extents based on dimension '" + + ModelArray::definedDimensions.at(dim).name + "'"); + } if (!cxios_is_defined_domain_ni_glo(domain)) { throw std::runtime_error( "Xios: Failed to set global x-size for domain '" + domainId + "'"); } - cxios_set_domain_ni(domain, (int)metadata.getLocalExtentX() + extension); if (!cxios_is_defined_domain_ni(domain)) { throw std::runtime_error( "Xios: Failed to set local x-size for domain '" + domainId + "'"); } - cxios_set_domain_ibegin(domain, (int)metadata.getLocalCornerX()); + cxios_set_domain_ibegin(domain, metadata.getLocalCornerX()); if (!cxios_is_defined_domain_ibegin(domain)) { throw std::runtime_error( "Xios: Failed to set local starting x-index for domain '" + domainId + "'"); @@ -811,17 +817,26 @@ void Xios::affixModelMetadata() "Xios: Failed to set longitude name for domain '" + domainId + "'"); } } else if (counter == 1) { - cxios_set_domain_nj_glo(domain, (int)metadata.getGlobalExtentY() + extension); + if (dim == ModelArray::Dimension::Y) { + cxios_set_domain_nj_glo(domain, metadata.getGlobalExtentY()); + cxios_set_domain_nj(domain, metadata.getLocalExtentY()); + } else if (dim == ModelArray::Dimension::YVERTEX) { + cxios_set_domain_nj_glo(domain, metadata.getGlobalExtentY() + 1); + cxios_set_domain_nj(domain, metadata.getLocalExtentY() + 1); + } else { + throw std::runtime_error( + "Xios: Could not set domain extents based on dimension '" + + ModelArray::definedDimensions.at(dim).name + "'"); + } if (!cxios_is_defined_domain_nj_glo(domain)) { throw std::runtime_error( "Xios: Failed to set global y-size for domain '" + domainId + "'"); } - cxios_set_domain_nj(domain, (int)metadata.getLocalExtentY() + extension); if (!cxios_is_defined_domain_nj(domain)) { throw std::runtime_error( "Xios: Failed to set local y-size for domain '" + domainId + "'"); } - cxios_set_domain_jbegin(domain, (int)metadata.getLocalCornerY()); + cxios_set_domain_jbegin(domain, metadata.getLocalCornerY()); if (!cxios_is_defined_domain_jbegin(domain)) { throw std::runtime_error( "Xios: Failed to set local starting y-index for domain '" + domainId + "'"); @@ -1655,8 +1670,12 @@ void Xios::write(const std::string fieldId, ModelArray& modelarray) } else if (type == ModelArray::Type::DG) { cxios_write_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], dims[1], ModelArray::size(ModelArray::Dimension::DG), -1); + } else if (type == ModelArray::Type::DGSTRESS) { + cxios_write_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], + dims[1], ModelArray::size(ModelArray::Dimension::DGSTRESS), -1); } else { - throw std::invalid_argument("Only HFields, VertexFields, and DGFields are supported"); + throw std::invalid_argument( + "Only HFields, VertexFields, DGFields, and DGSFields are supported"); } } @@ -1688,8 +1707,12 @@ void Xios::read(const std::string fieldId, ModelArray& modelarray) } else if (type == ModelArray::Type::DG) { cxios_read_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], dims[1], ModelArray::size(ModelArray::Dimension::DG)); + } else if (type == ModelArray::Type::DGSTRESS) { + cxios_read_data_k83(fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], + dims[1], ModelArray::size(ModelArray::Dimension::DGSTRESS)); } else { - throw std::invalid_argument("Only HFields, VertexFields, and DGFields are supported"); + throw std::invalid_argument( + "Only HFields, VertexFields, DGFields, and DGSFields are supported"); } } } diff --git a/core/src/discontinuousgalerkin/ModelArrayDetails.cpp b/core/src/discontinuousgalerkin/ModelArrayDetails.cpp index adfc1f2f1..bd2ff7ccb 100644 --- a/core/src/discontinuousgalerkin/ModelArrayDetails.cpp +++ b/core/src/discontinuousgalerkin/ModelArrayDetails.cpp @@ -84,7 +84,7 @@ const std::map ModelArray::typeNames = { { ModelArray::Type::U, "UField" }, { ModelArray::Type::V, "VField" }, { ModelArray::Type::DG, "DGField" }, - { ModelArray::Type::DGSTRESS, "DGStressField" }, + { ModelArray::Type::DGSTRESS, "DGSField" }, { ModelArray::Type::CG, "CGField" }, }; diff --git a/core/src/include/Xios.hpp b/core/src/include/Xios.hpp index b7bea1cf6..471b59cbe 100644 --- a/core/src/include/Xios.hpp +++ b/core/src/include/Xios.hpp @@ -180,10 +180,12 @@ class Xios : public Configured { std::map axisIds = { { ModelArray::Type::VERTEX, "VertexAxis" }, { ModelArray::Type::DG, "DGAxis" }, + { ModelArray::Type::DGSTRESS, "DGSAxis" }, }; std::map axisNames = { { "VertexAxis", "ncoords" }, { "DGAxis", "dg_comp" }, + { "DGSAxis", "dgstress_comp" }, }; xios::CAxisGroup* getAxisGroup(); xios::CAxis* getAxis(const std::string axisId); @@ -193,6 +195,7 @@ class Xios : public Configured { { ModelArray::Type::H, "HDomain" }, { ModelArray::Type::VERTEX, "VertexDomain" }, { ModelArray::Type::DG, "HDomain" }, + { ModelArray::Type::DGSTRESS, "HDomain" }, }; xios::CDomainGroup* getDomainGroup(); xios::CDomain* getDomain(std::string domainId); @@ -217,6 +220,7 @@ class Xios : public Configured { { ModelArray::Type::H, "HGrid" }, { ModelArray::Type::VERTEX, "VertexGrid" }, { ModelArray::Type::DG, "DGGrid" }, + { ModelArray::Type::DGSTRESS, "DGSGrid" }, }; /* File */ diff --git a/core/test/XiosRead_test.cpp b/core/test/XiosRead_test.cpp index 8b0d7957d..2e05cf7df 100644 --- a/core/test/XiosRead_test.cpp +++ b/core/test/XiosRead_test.cpp @@ -24,6 +24,7 @@ const std::string restartFilename = testFilesDir + "/xios_test_input.nc"; const std::string forcingFilename = testFilesDir + "/xios_test_forcing.nc"; static const int DG = 3; +static const int DGSTRESSCOMP = 8; namespace Nextsim { @@ -46,7 +47,8 @@ MPI_TEST_CASE("TestXiosRead", 2) config << "restart_period = P0-0T01:30:00" << std::endl; config << "partition_file = xios_test_partition_metadata_2.nc" << std::endl; config << "[XiosInput]" << std::endl; - config << "field_names = " << maskName << "," << coordsName << "," << hiceName << std::endl; + config << "field_names = " << maskName << "," << coordsName << "," << hiceName << "," + << ticeName << std::endl; config << "[XiosForcing]" << std::endl; config << "filename = " << forcingFilename << std::endl; config << "field_names = " << uName << std::endl; @@ -76,8 +78,10 @@ MPI_TEST_CASE("TestXiosRead", 2) // TODO: We could deduce this from the NetCDF file ModelArray::setNComponents(ModelArray::Type::DG, DG); + ModelArray::setNComponents(ModelArray::Type::DGSTRESS, DGSTRESSCOMP); ModelArray::setNComponents(ModelArray::Type::VERTEX, ModelArray::nCoords); REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DG); + REQUIRE(ModelArray::nComponents(ModelArray::Type::DGSTRESS) == DGSTRESSCOMP); REQUIRE(ModelArray::nComponents(ModelArray::Type::VERTEX) == ModelArray::nCoords); // Affix ModelMetadata to Xios handler @@ -119,6 +123,11 @@ MPI_TEST_CASE("TestXiosRead", 2) float expected = 1.0 * (d + DG * (i + nx * j)); REQUIRE(entry.second.components({ i, j })[d] == doctest::Approx(expected)); } + } else if (entry.first == ticeName) { + for (size_t d = 0; d < DGSTRESSCOMP; ++d) { + REQUIRE(entry.second.components({ i, j })[d] + == doctest::Approx(2.0 * (d + DGSTRESSCOMP * (i + nx * j)))); + } } } } @@ -149,5 +158,4 @@ MPI_TEST_CASE("TestXiosRead", 2) xiosHandler.context_finalize(); Finalizer::finalize(); } - } diff --git a/core/test/XiosWrite_test.cpp b/core/test/XiosWrite_test.cpp index 38dd925ab..254c2ff0c 100644 --- a/core/test/XiosWrite_test.cpp +++ b/core/test/XiosWrite_test.cpp @@ -24,6 +24,7 @@ const std::string restartFilename = testFilesDir + "/xios_test_output.nc"; const std::string diagnosticFilename = testFilesDir + "/xios_test_diagnostic.nc"; static const int DG = 3; +static const int DGSTRESSCOMP = 8; namespace Nextsim { @@ -46,7 +47,8 @@ MPI_TEST_CASE("TestXiosWrite", 2) config << "partition_file = xios_test_partition_metadata_2.nc" << std::endl; config << "restart_period = P0-0T01:30:00" << std::endl; config << "[XiosOutput]" << std::endl; - config << "field_names = " << maskName << "," << coordsName << "," << hiceName << std::endl; + config << "field_names = " << maskName << "," << coordsName << "," << hiceName << "," + << ticeName << std::endl; config << "[XiosDiagnostic]" << std::endl; config << "filename = " << diagnosticFilename << std::endl; config << "field_names = " << uName << std::endl; @@ -84,8 +86,10 @@ MPI_TEST_CASE("TestXiosWrite", 2) ModelArray::setDimension(ModelArray::Dimension::XVERTEX, nx_glo + 1, nx + 1, 0); ModelArray::setDimension(ModelArray::Dimension::YVERTEX, ny_glo + 1, ny + 1, 0); ModelArray::setNComponents(ModelArray::Type::DG, DG); + ModelArray::setNComponents(ModelArray::Type::DGSTRESS, DGSTRESSCOMP); ModelArray::setNComponents(ModelArray::Type::VERTEX, ModelArray::nCoords); REQUIRE(ModelArray::nComponents(ModelArray::Type::DG) == DG); + REQUIRE(ModelArray::nComponents(ModelArray::Type::DGSTRESS) == DGSTRESSCOMP); REQUIRE(ModelArray::nComponents(ModelArray::Type::VERTEX) == ModelArray::nCoords); // Affix ModelMetadata to Xios handler @@ -96,6 +100,7 @@ MPI_TEST_CASE("TestXiosWrite", 2) xiosHandler.setFieldType(maskName, ModelArray::Type::H); xiosHandler.setFieldType(coordsName, ModelArray::Type::VERTEX); xiosHandler.setFieldType(hiceName, ModelArray::Type::DG); + xiosHandler.setFieldType(ticeName, ModelArray::Type::DGSTRESS); xiosHandler.setFieldType(uName, ModelArray::Type::H); // Set file split frequency for restarts (but not diagnostics) @@ -130,6 +135,15 @@ MPI_TEST_CASE("TestXiosWrite", 2) } } } + DGSField tice(ModelArray::Type::DGSTRESS); + tice.resize(); + for (size_t j = 0; j < ny; ++j) { + for (size_t i = 0; i < nx; ++i) { + for (size_t d = 0; d < DGSTRESSCOMP; ++d) { + tice.components({ i, j })[d] = 2.0 * (d + DGSTRESSCOMP * (i + nx * j)); + } + } + } HField u(ModelArray::Type::H); u.resize(); @@ -160,6 +174,7 @@ MPI_TEST_CASE("TestXiosWrite", 2) { maskName, mask }, { coordsName, coordinates }, { hiceName, hice }, + { ticeName, tice }, }, {} }; ModelState diagnostics = { { diff --git a/core/test/xios_test_input.cdl b/core/test/xios_test_input.cdl index 1d027ed1e..a5cd6d585 100644 --- a/core/test/xios_test_input.cdl +++ b/core/test/xios_test_input.cdl @@ -36,6 +36,11 @@ variables: DGAxis:standard_name = "DG DoFs" ; DGAxis:long_name = "Discontinuous Galerkin number of degrees of freedom" ; DGAxis:units = "-" ; + float DGSAxis(dgstress_comp) ; + DGSAxis:axis = "DGS" ; + DGSAxis:standard_name = "DGS DoFs" ; + DGSAxis:long_name = "Discontinuous Galerkin stress component number of degrees of freedom" ; + DGSAxis:units = "-" ; float VertexAxis(ncoords) ; VertexAxis:axis = "DIM" ; VertexAxis:standard_name = "Dimension" ; @@ -50,6 +55,9 @@ variables: float mask(y, x) ; mask:online_operation = "once" ; mask:coordinates = "lat lon" ; + float tice(y, x, dgstress_comp) ; + tice:online_operation = "once" ; + tice:coordinates = "lat lon DGSAxis" ; // global attributes: :structure_name = "parametric_rectangular" ; @@ -98,4 +106,14 @@ data: mask = 0, 0, 0, 0, 1, 1, 1, 1 ; + + tice = + 0, 2, 4, 6, 8, 10, 12, 14, + 16, 18, 20, 22, 24, 26, 28, 30, + 0, 2, 4, 6, 8, 10, 12, 14, + 16, 18, 20, 22, 24, 26, 28, 30, + 32, 34, 36, 38, 40, 42, 44, 46, + 48, 50, 52, 54, 56, 58, 60, 62, + 32, 34, 36, 38, 40, 42, 44, 46, + 48, 50, 52, 54, 56, 58, 60, 62 ; } diff --git a/docs/xios.rst b/docs/xios.rst index ffaa30cc7..3aa7084e4 100644 --- a/docs/xios.rst +++ b/docs/xios.rst @@ -78,6 +78,11 @@ As elsewhere in the model, these configuration values are parsed by calling the building with MPI, you will need to pass the MPI communicator to this member function when calling it. +The XIOS I/O implementation supports fields of ``HField``, ``VertexField``, +``DGField``, and ``DGSField``. When reading from file, the field type will +be deduced from its dimension. When writing to file, the field type should be +set using the ``setFieldType`` member function of the ``Xios`` handler class. + Developer notes ^^^^^^^^^^^^^^^