Get the ice-ocean stress into the ModelArrayRef system (#754)

# Get the ice-ocean stress into the ModelArrayRef system
## Fixes #649 

### Task List
- [ ] Defined the tests that specify a complete and functioning change
(*It may help to create a [design specification & test
specification](../../../wiki/Specification-Template)*)
- [x] Implemented the source code change that satisfies the tests
- [ ] Documented the feature by providing worked example
- [ ] Updated the README or other documentation
- [ ] Completed the pre-Request checklist below

---
# Change Description

Introduces a getIceOceanStress function into CGDynamicsKernel, called by
getDG0Data when that function gets the right arguments. All derived
DynamicsKernels (VPCGDynamicsKernel, BrittleCGDynamicsKernel, and
FreeDriftDynamicsKernel) must implement getIceOceanStress.

This is a partial re-implementation of code previously rejected from the
develop branch.

---
# Test Description

No CI test. I just ran the benchmark code and checked that the values
returned are of the right order of magnitude.

---
# Documentation Impact

None

---
# Other Details

None

---
### Pre-Request Checklist

- [x] The requirements of this pull request are fully captured in an
issue or design specification and are linked and summarised in the
description of this PR
- [x] No new warnings are generated
- [x] The documentation has been updated (or an issue has been created
to track the corresponding change)
- [x] Methods and Tests are commented such that they can be understood
without having to obtain additional context
- [x] This PR/Issue is labelled as a bug/feature/enhancement/breaking
change
- [x] File dates have been updated to reflect modification date
- [x] This change conforms to the conventions described in the README

core/src/include/ModelComponent.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file ModelComponent.hpp
  *
- * @date 1 Jul 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -64,6 +64,8 @@ namespace Protected {
     inline constexpr TextTag WIND_V = "WIND_V"; // y(north)-ward component of wind, m s⁻¹
     inline constexpr TextTag ICE_U = "ICE_U"; // x(east)-ward ice velocity, m s⁻¹
     inline constexpr TextTag ICE_V = "ICE_V"; // y(north)-ward ice velocity, m s⁻¹
+    inline constexpr TextTag IO_STRESS_X = "IO_STRESS_X"; // x(east)-ward ice-ocean stress, Pa
+    inline constexpr TextTag IO_STRESS_Y = "IO_STRESS_Y"; // y(north)-ward ice-ocean stress, Pa
     // Slab ocean fields
     inline constexpr TextTag SLAB_SST = "SLAB_SST"; // Slab ocean sea surface temperature, ˚C
     inline constexpr TextTag SLAB_SSS = "SLAB_SSS"; // Slab ocean sea surface salinity, ˚C

core/src/include/gridNames.hpp

@@ -29,6 +29,9 @@ static const std::string vWindName = "vwind";
 static const std::string uOceanName = "uocean";
 static const std::string vOceanName = "vocean";
 
+static const std::string uIOStressName = "uiostress";
+static const std::string vIOStressName = "viostress";
+
 static const std::string coordsName = "coords";
 static const std::string latitudeName = "latitude";
 static const std::string longitudeName = "longitude";

core/src/modules/DynamicsModule/BBMDynamics.cpp

@@ -141,6 +141,9 @@ void BBMDynamics::update(const TimestepTime& tst)
 
     uice = kernel.getDG0Data(uName);
     vice = kernel.getDG0Data(vName);
+
+    taux = kernel.getDG0Data(uIOStressName);
+    tauy = kernel.getDG0Data(vIOStressName);
 }
 
 // All data for prognostic output

core/src/modules/DynamicsModule/MEVPDynamics.cpp

@@ -108,6 +108,9 @@ void MEVPDynamics::update(const TimestepTime& tst)
 
     uice = kernel.getDG0Data(uName);
     vice = kernel.getDG0Data(vName);
+
+    taux = kernel.getDG0Data(uIOStressName);
+    tauy = kernel.getDG0Data(vIOStressName);
 }
 
 ModelState MEVPDynamics::getStateRecursive(const OutputSpec& os) const

core/src/modules/include/IDynamics.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file IDynamics.hpp
  *
- * @date 7 Sep 2023
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -23,6 +23,8 @@ class IDynamics : public ModelComponent {
         : uice(ModelArray::Type::H)
         , vice(ModelArray::Type::H)
         , damage(ModelArray::Type::H)
+        , taux(ModelArray::Type::H)
+        , tauy(ModelArray::Type::H)
         , hice(getStore())
         , cice(getStore())
         , hsnow(getStore())
@@ -34,6 +36,8 @@ class IDynamics : public ModelComponent {
         , m_usesDamage(usesDamageIn)
     {
         getStore().registerArray(Shared::DAMAGE, &damage, RW);
+        getStore().registerArray(Protected::IO_STRESS_X, &taux, RO);
+        getStore().registerArray(Protected::IO_STRESS_Y, &tauy, RO);
     }
     virtual ~IDynamics() = default;
 
@@ -76,12 +80,14 @@ class IDynamics : public ModelComponent {
     HField vice;
     // Updated damage array
     HField damage;
+    // Ice-ocean stress (for the coupler, mostly)
+    HField taux;
+    HField tauy;
     // References to the DG0 finite volume data arrays
     ModelArrayRef<Shared::H_ICE, RW> hice;
     ModelArrayRef<Shared::C_ICE, RW> cice;
     ModelArrayRef<Shared::H_SNOW, RW> hsnow;
     ModelArrayRef<Protected::DAMAGE, RO> damage0;
-    // ModelArrayRef<ModelComponent::SharedArray::D, MARBackingStore, RW> damage;
 
     // References to the forcing velocity arrays
     ModelArrayRef<Protected::WIND_U> uwind;

core/src/modules/include/ProtectedArrayNames.ipp

@@ -40,3 +40,5 @@
     { "sss_slab", "SLAB_SSS" }, // Slab ocean surface salinity PSU
     { "qdw", "SLAB_QDW" }, // Slab ocean temperature nudging heat flux, W m⁻²
     { "fdw", "SLAB_FDW" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+    { "taux", "IO_STRESS_X" }, // Ice-ocean stress x(east) direction, Pa
+    { "tauy", "IO_STRESS_Y" }, // Ice-ocean stress y(north) direction, Pa

dynamics/src/CGDynamicsKernel.cpp

@@ -98,6 +98,16 @@ ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) co
         DGVector<DGadvection> vtmp(*smesh);
         Nextsim::Interpolations::CG2DG(*smesh, vtmp, v);
         return DGModelArray::dg2ma(vtmp, data);
+    } else if (name == uIOStressName) {
+        ModelArray data(ModelArray::Type::U);
+        DGVector<DGadvection> utmp(*smesh);
+        Nextsim::Interpolations::CG2DG(*smesh, utmp, getIceOceanStress(name));
+        return DGModelArray::dg2ma(utmp, data);
+    } else if (name == vIOStressName) {
+        ModelArray data(ModelArray::Type::V);
+        DGVector<DGadvection> vtmp(*smesh);
+        Nextsim::Interpolations::CG2DG(*smesh, vtmp, getIceOceanStress(name));
+        return DGModelArray::dg2ma(vtmp, data);
     } else {
         return DynamicsKernel<DGadvection, DGstressComp>::getDG0Data(name);
     }

dynamics/src/include/BrittleCGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file BrittleCGDynamicsKernel.hpp
  *
- * @date 19 Nov 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -163,6 +163,22 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
         }
     }
 
+    double getIceOceanStressElement(const std::string& name, const int i) const override
+    {
+        const double FOcean = params.COcean * params.rhoOcean;
+
+        const double uOceanRel = uOcean(i) - avgU(i);
+        const double vOceanRel = vOcean(i) - avgV(i);
+        const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel);
+
+        if (name == uIOStressName)
+            return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+        else if (name == vIOStressName)
+            return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        else
+            return std::numeric_limits<double>::quiet_NaN();
+    }
+
 protected:
     CGVector<CGdegree> avgU;
     CGVector<CGdegree> avgV;

dynamics/src/include/CGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file CGDynamicsKernel.hpp
  *
- * @date Jan 31, 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -47,6 +47,23 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
     void applyBoundaries() override;
     void prepareAdvection() override;
 
+    virtual inline double getIceOceanStressElement(const std::string& name, const int i) const = 0;
+    CGVector<CGdegree> getIceOceanStress(const std::string& name) const
+    {
+        if (name != uIOStressName && name != vIOStressName)
+            throw std::logic_error(std::string(__func__) + " called with an unknown argument "
+                + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n");
+
+        CGVector<CGdegree> tau;
+        tau.resizeLike(u);
+
+#pragma omp parallel for
+        for (int i = 0; i < tau.rows(); ++i)
+            tau(i) = getIceOceanStressElement(name, i);
+
+        return tau;
+    }
+
 protected:
     void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy);
     void dirichletZero(CGVector<CGdegree>&) const;

dynamics/src/include/FreeDriftDynamicsKernel.hpp

@@ -4,7 +4,7 @@
  * Implementation of "classic free drift", where we ignore all \rho h terms in the momentum
  * equation. This is equivalent to assuming that the ice is very thin.
  *
- * @date 19 Nov 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -69,6 +69,22 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
                 + NansenNumber * (-uAtmos(i) * sinOceanAngle + vAtmos(i) * cosOceanAngle);
         }
     }
+
+    double getIceOceanStressElement(const std::string& name, const int i) const override
+    {
+        const double FOcean = params.COcean * params.rhoOcean;
+
+        const double uOceanRel = uOcean(i) - u(i);
+        const double vOceanRel = vOcean(i) - v(i);
+        const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel);
+
+        if (name == uIOStressName)
+            return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+        else if (name == vIOStressName)
+            return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        else
+            return std::numeric_limits<double>::quiet_NaN();
+    }
 };
 } /* namespace Nextsim */
 

dynamics/src/include/VPCGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file VPCGDynamicsKernel.hpp
  *
- * @date 19 Nov 2024
+ * @date 06 Dec 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -89,6 +89,22 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG
         DynamicsKernel<DGadvection, DGstressComp>::update(tst);
     }
 
+    double getIceOceanStressElement(const std::string& name, const int i) const override
+    {
+        const double FOcean = params.COcean * params.rhoOcean;
+
+        double uOcnRel = u(i) - uOcean(i);
+        double vOcnRel = v(i) - vOcean(i);
+        double absocn = sqrt(SQR(uOcnRel) + SQR(vOcnRel));
+
+        if (name == uIOStressName)
+            return FOcean * absocn * uOcnRel;
+        else if (name == vIOStressName)
+            return FOcean * absocn * vOcnRel;
+        else
+            return std::numeric_limits<double>::quiet_NaN();
+    }
+
 protected:
     StressUpdateStep<DGadvection, DGstressComp>& stressStep;
     const VPParameters& params;