Merge remote-tracking branch 'origin/develop' into issue459_oasis_oce…

…an_coupling

.github/workflows/clang-compile-tests.yml

@@ -22,7 +22,7 @@ jobs:
         cd core/src
         clang-format --dry-run -Werror *cpp include/*hpp
         cd -
-        for component in physics
+        for component in physics dynamics
         do
           cd $component/src
           clang-format --dry-run -Werror *cpp include/*hpp
@@ -51,7 +51,7 @@ jobs:
         apt install -y wget
         cd build
         (cd core/test && wget "ftp://ftp.nersc.no/nextsim/netCDF/partition_metadata_1.nc")
-        for component in core physics
+        for component in core physics dynamics
         do
             cd $component/test
             for file in $(find test* -maxdepth 0 -type f | grep -vP "_MPI\d*$")
@@ -88,7 +88,7 @@ jobs:
         apt install -y wget
         cd build
         (cd core/test && wget "ftp://ftp.nersc.no/nextsim/netCDF/partition_metadata_1.nc")
-        for component in core physics
+        for component in core physics dynamics
         do
             cd $component/test
             for file in $(find test* -maxdepth 0 -type f | grep -P "_MPI\d+$")
@@ -134,7 +134,7 @@ jobs:
     - name: run serial tests
       run: |
         cd build
-        for component in core physics
+        for component in core physics dynamics
         do
             cd $component/test
             for file in $(find test* -maxdepth 0 -type f | grep -v "_MPI\d*$")

CMakeLists.txt

@@ -90,14 +90,54 @@ set(ScriptDirectory "${PROJECT_SOURCE_DIR}/scripts")
 set(ModuleBuilderScript "${ScriptDirectory}/module_builder.py")
 set(ModuleHeaderScript "${ScriptDirectory}/module_header.py")
 
+# Dynamics type. Defaults to DG1.
+# Available options are:
+# DG1 Discontinuous Galerkin, degree 1 (3 components)
+# DG2 Discontinuous Galerkin, degree 2 (6 components)
+set(isDG FALSE)
+
+# Replace a missing dynamics type with DG1
+if (NOT DynamicsType)
+    set(DynamicsType "DG2")
+    message("Defaulting dynamics to ${DynamicsType}")
+endif()
+
+message("Dynamics is set to ${DynamicsType}")
+
+if (DynamicsType STREQUAL "DG1")
+    set(isDG TRUE)
+    set(DGComp 3)
+    set(CGDegree 1)
+endif()
+
+if (DynamicsType STREQUAL "DG2")
+    set(isDG TRUE)
+    set(DGComp 6)
+    set(CGDegree 2)
+endif()
+
 # Set the components which provide source or object code to the main model
 set(CodeComponents
     "physics"
-    "dynamics"
     )
 
+if (isDG)
+# Add the DG dynamics subdirectory to the list of code components
+    set(CodeComponents
+        "${CodeComponents}"
+        "dynamics"
+        )
+
+# Set the number of DG stress components given the CG degree
+    if (CGDegree EQUAL 1)
+        set(DGStressComp 3)
+    elseif(CGDegree EQUAL 2)
+        set(DGStressComp 8)
+    else()
+        message("Invalid value of CGDegree. Valid values are 1–2")
+    endif()
+
 # The exact selection of dimensions and Types available to ModelArray
-if(True)
     set(ModelArrayStructure "discontinuousgalerkin")
 endif()
 
@@ -124,6 +164,9 @@ target_compile_definitions(nextsimlib
     PRIVATE
     $<$<BOOL:${ENABLE_MPI}>:USE_MPI>
     $<$<BOOL:${ENABLE_XIOS}>:USE_XIOS>
+    DGCOMP=${DGComp}
+    DGSTRESSCOMP=${DGStressComp}
+    CGDEGREE=${CGDegree}
     )
 target_include_directories(nextsimlib
     PRIVATE

core/src/DevStep.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file DevStep.cpp
  *
- * @date Jan 12, 2022
+ * @date 2 Jul 2024
  * @author Tim Spain <[email protected]>
  */
 
@@ -26,7 +26,13 @@ void DevStep::init()
     tryConfigure(ido);
 }
 
-void DevStep::start(const TimePoint& startTime) { lastOutput = startTime; }
+void DevStep::start(const TimePoint& startTime)
+{
+    // Set the last output time for the restart files to the model start time
+    lastOutput = startTime;
+    // Set the model start time for the diagnostic output files
+    Module::getImplementation<IDiagnosticOutput>().setModelStart(startTime);
+}
 
 void DevStep::iterate(const TimestepTime& tst)
 {
@@ -36,11 +42,12 @@ void DevStep::iterate(const TimestepTime& tst)
     mData->incrementTime(tst.step);
     if ((m_restartPeriod.seconds() > 0) && (mData->time() >= lastOutput + m_restartPeriod)) {
         std::string currentFileName = mData->time().format(m_restartFileName);
-        pData->writeRestartFile(currentFileName);
+        pData->writeRestartFile(currentFileName, *mData);
         lastOutput = mData->time();
     }
     // XIOS wants all the fields, every timestep, so I guess that's what everyone gets
-    ModelState overallState = pData->getStateRecursive(true);
+    OutputSpec os; // The default OutputSpec is all fields, but only cell average values
+    ModelState overallState = pData->getStateRecursive(os);
     overallState.merge(ConfiguredModule::getAllModuleConfigurations());
     Module::getImplementation<IDiagnosticOutput>().outputState(*mData);
 }

core/src/MissingData.cpp

@@ -1,15 +1,16 @@
 /*!
  * @file MissingData.cpp
  *
- * @date Jun 14, 2022
+ * @date Jul 17, 2024
  * @author Tim Spain <[email protected]>
+ * @author Einar Ólason <[email protected]>
  */
 
 #include "include/MissingData.hpp"
 
 namespace Nextsim {
 
-const double MissingData::defaultValue = -0x1p300;
+const double MissingData::defaultValue = 1.7e38;
 double MissingData::value = MissingData::defaultValue;
 
 } /* namespace Nextsim */

core/src/Model.cpp

@@ -193,17 +193,7 @@ void Model::run() { iterator.run(); }
 void Model::writeRestartFile()
 {
     std::string formattedFileName = m_etadata.time().format(finalFileName);
-
-    Logged::notice(std::string("  Writing state-based restart file: ") + formattedFileName + '\n');
-    // Copy the configuration from the ModelState to the ModelMetadata
-    ConfigMap modelConfig = getConfig();
-    modelConfig.merge(pData.getStateRecursive(true).config);
-    modelConfig.merge(ConfiguredModule::getAllModuleConfigurations());
-    m_etadata.setConfig(modelConfig);
-    // Get the model state from PrognosticData and add the coordinates.
-    ModelState state = pData.getState();
-    m_etadata.affixCoordinates(state);
-    StructureFactory::fileFromState(state, m_etadata, formattedFileName, true);
+    pData.writeRestartFile(formattedFileName, m_etadata);
 }
 
 ModelMetadata& Model::metadata() { return m_etadata; }

core/src/ModelMetadata.cpp

@@ -10,6 +10,13 @@
 #include "include/StructureModule.hpp"
 #include "include/gridNames.hpp"
 
+#ifdef USE_MPI
+#include <ncDim.h>
+#include <ncFile.h>
+#include <ncGroup.h>
+#include <ncVar.h>
+#endif
+
 namespace Nextsim {
 
 const std::string& ModelMetadata::structureName() const
@@ -33,6 +40,7 @@ void ModelMetadata::setMpiMetadata(MPI_Comm comm)
 
 void ModelMetadata::getPartitionMetadata(std::string partitionFile)
 {
+    // TODO: Move the reading of the partition file to its own class
     netCDF::NcFile ncFile(partitionFile, netCDF::NcFile::read);
     int sizes = ncFile.getDim("L").getSize();
     int nBoxes = ncFile.getDim("P").getSize();

core/src/PDWriter.cpp

@@ -13,14 +13,19 @@
 #include "include/StructureFactory.hpp"
 
 namespace Nextsim {
-void PrognosticData::writeRestartFile(const std::string& filePath) const
+void PrognosticData::writeRestartFile(
+    const std::string& filePath, const ModelMetadata& metadata) const
 {
-    ConfigMap modelConfig;
-    modelConfig.merge(ModelConfig::getConfig());
+    Logged::notice(std::string("  Writing state-based restart file: ") + filePath + '\n');
+
+    ConfigMap modelConfig = ModelConfig::getConfig();
     modelConfig.merge(getStateRecursive(true).config);
     modelConfig.merge(ConfiguredModule::getAllModuleConfigurations());
-    ModelMetadata meta;
+    ModelMetadata meta(metadata);
     meta.setConfig(modelConfig);
-    StructureFactory::fileFromState(getState(), meta, filePath);
+    ModelState state = getState();
+    meta.affixCoordinates(state);
+
+    StructureFactory::fileFromState(state, meta, filePath, true);
 }
 }

core/src/ParaGridIO.cpp

@@ -27,11 +27,16 @@
 
 namespace Nextsim {
 
+// Accept both post-May 2024 (xdim, ydim, zdim) dimension names and pre-May 2024 (x, y, z)
 const std::map<std::string, ModelArray::Type> ParaGridIO::dimensionKeys = {
     { "yx", ModelArray::Type::H },
+    { "ydimxdim", ModelArray::Type::H },
     { "zyx", ModelArray::Type::Z },
+    { "zdimydimxdim", ModelArray::Type::Z },
     { "yxdg_comp", ModelArray::Type::DG },
+    { "ydimxdimdg_comp", ModelArray::Type::DG },
     { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
+    { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
     { "ycgxcg", ModelArray::Type::CG },
     { "yvertexxvertexncoords", ModelArray::Type::VERTEX },
 };
@@ -90,7 +95,18 @@ ModelState ParaGridIO::getModelState(const std::string& filePath)
                 continue;
 
             ModelArray::DimensionSpec& dimensionSpec = entry.second;
+            // Find dimensions in the netCDF file by their name in the ModelArray details
             netCDF::NcDim dim = dataGroup.getDim(dimensionSpec.name);
+            // Also check the old name
+            if (dim.isNull()) {
+                dim = dataGroup.getDim(dimensionSpec.altName);
+            }
+            // If we didn't find a dimension with the dimensions name or altName, throw.
+            if (dim.isNull()) {
+                throw std::out_of_range(
+                    std::string("No netCDF dimension found corresponding to the dimension named ")
+                    + dimensionSpec.name + std::string(" or ") + dimensionSpec.altName);
+            }
             if (entry.first == ModelArray::Dimension::Z) {
                 // A special case, as the number of levels in the file might not be
                 // the number that the selected ice thermodynamics requires.

core/src/PrognosticData.cpp

@@ -1,8 +1,9 @@
 /*!
  * @file PrognosticData.cpp
  *
- * @date 7 Sep 2023
+ * @date 1 Jul 2024
  * @author Tim Spain <[email protected]>
+ * @author Einar Ólason <[email protected]>
  */
 
 #include "include/PrognosticData.hpp"
@@ -19,6 +20,7 @@ PrognosticData::PrognosticData()
     , m_conc(ModelArray::Type::H)
     , m_snow(ModelArray::Type::H)
     , m_tice(ModelArray::Type::Z)
+    , m_damage(ModelArray::Type::H)
     , pAtmBdy(0)
     , pOcnBdy(0)
     , pDynamics(0)
@@ -58,12 +60,12 @@ void PrognosticData::setData(const ModelState::DataMap& ms)
     m_conc = ms.at("cice");
     m_tice = ms.at("tice");
     m_snow = ms.at("hsnow");
-    // Damage is an optional field, and defaults to zero, if absent
+    // Damage is an optional field, and defaults to 1, if absent
     if (ms.count(damageName) > 0) {
         m_damage = ms.at(damageName);
     } else {
         m_damage.resize();
-        m_damage = 0.5;
+        m_damage = 1.;
     }
 
     pAtmBdy->setData(ms);
@@ -113,10 +115,14 @@ void PrognosticData::updatePrognosticFields()
     m_damage.setData(damageUpd);
 }
 
+// Gets all of the prognostic data, including that in the dynamics
 ModelState PrognosticData::getState() const
 {
     ModelArrayRef<Protected::SST> sst(getStore());
     ModelArrayRef<Protected::SSS> sss(getStore());
+
+    // Get the prognostic data from the dynamics, including the full dynamics state
+    ModelState dynamicsState = pDynamics->getState();
     // clang-format off
     ModelState localState = { {
                  { "mask", ModelArray(oceanMask()) }, // make a copy
@@ -129,27 +135,33 @@ ModelState PrognosticData::getState() const
              },
         {} };
     // clang-format on
-    // Get the state from the dynamics (ice velocity). This allows the
-    // dynamics to define its own dimensions for the velocity grid.
-    localState.merge(pDynamics->getState());
-
-    // Merge in the damage field, if the dynamics uses it.
-    if (pDynamics->usesDamage()) {
-        ModelState damageState = { {
-                                       { "damage", mask(m_damage) },
-                                   },
-            {} };
-        localState.merge(damageState);
-    }
-    return localState;
+
+    // Use the dynamics values of any duplicated fields
+    ModelState state(dynamicsState);
+    state.merge(localState);
+
+    return state;
 }
 
+// Recursively gets the data from all subcomponents
 ModelState PrognosticData::getStateRecursive(const OutputSpec& os) const
 {
-    ModelState state(getState());
-    state.merge(pAtmBdy->getStateRecursive(os));
-    state.merge(iceGrowth.getStateRecursive(os));
-    state.merge(pDynamics->getStateRecursive(os));
+    ModelState state;
+    /* If allComponents is set on the OutputSpec, then for any duplicate fields, the subsystems
+     * take priority, otherwise the fields held by PrognosticData itself. Note that std::map::merge
+     * will not overwrite existing keys, so the first one that exists will survive.
+     */
+    if (os.allComponents()) {
+        state.merge(pAtmBdy->getStateRecursive(os));
+        state.merge(iceGrowth.getStateRecursive(os));
+        state.merge(pDynamics->getStateRecursive(os));
+        state.merge(getState());
+    } else {
+        state.merge(getState());
+        state.merge(pAtmBdy->getStateRecursive(os));
+        state.merge(iceGrowth.getStateRecursive(os));
+        state.merge(pDynamics->getStateRecursive(os));
+    }
     // OceanBoundary does not contribute to the output model state
     return os ? state : ModelState();
 }