Hotfix for issue #720

dynamics/src/CGDynamicsKernel.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file CGDynamicsKernel.cpp
  *
- * @date 07 Oct 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar ÓLason <einar.olason@nersc.no>
  */
@@ -105,12 +105,12 @@ ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) co
     } else if (name == uIOStressName) {
         ModelArray data(ModelArray::Type::U);
         DGVector<DGadvection> utmp(*smesh);
-        Nextsim::Interpolations::CG2DG(*smesh, utmp, getIceOceanStress().first);
+        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().second);
+        Nextsim::Interpolations::CG2DG(*smesh, vtmp, getIceOceanStress(name));
         return DGModelArray::dg2ma(vtmp, data);
     } else {
         return DynamicsKernel<DGadvection, DGstressComp>::getDG0Data(name);
@@ -443,28 +443,6 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::applyBoundaries()
     // TODO Periodic boundary conditions.
 }
 
-template <int DGadvection>
-std::pair<CGVector<CGdegree>, CGVector<CGdegree>>
-CGDynamicsKernel<DGadvection>::getIceOceanStress() const
-{
-    CGVector<CGdegree> taux, tauy;
-    taux.resizeLike(uStress);
-    tauy.resizeLike(vStress);
-
-#pragma omp parallel for
-    for (int i = 0; i < taux.rows(); ++i) {
-        const double uOceanRel = uOcean(i) - uStress(i);
-        const double vOceanRel = vOcean(i) - vStress(i);
-        const double cPrime = DynamicsKernel<DGadvection, DGstressComp>::getParams().F_ocean
-            * std::hypot(uOceanRel, vOceanRel);
-
-        taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
-        tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
-    }
-
-    return { taux, tauy };
-}
-
 // Instantiate the templates for all (1, 3, 6) degrees of DGadvection
 template class CGDynamicsKernel<1>;
 template class CGDynamicsKernel<3>;

dynamics/src/include/BrittleCGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file BrittleCGDynamicsKernel.hpp
  *
- * @date 07 Oct 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -38,8 +38,6 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
     using DynamicsKernel<DGadvection, DGstressComp>::applyBoundaries;
     using DynamicsKernel<DGadvection, DGstressComp>::advectionAndLimits;
     using DynamicsKernel<DGadvection, DGstressComp>::dgtransport;
-    using DynamicsKernel<DGadvection, DGstressComp>::cosOceanAngle;
-    using DynamicsKernel<DGadvection, DGstressComp>::sinOceanAngle;
 
     using CGDynamicsKernel<DGadvection>::u;
     using CGDynamicsKernel<DGadvection>::v;
@@ -58,11 +56,12 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
     using CGDynamicsKernel<DGadvection>::dStressY;
     using CGDynamicsKernel<DGadvection>::pmap;
 
+    double cosOceanAngle, sinOceanAngle;
+
 public:
     BrittleCGDynamicsKernel(StressUpdateStep<DGadvection, DGstressComp>& stressStepIn,
         const DynamicsParameters& paramsIn)
-        : CGDynamicsKernel<DGadvection>(cos(radians * paramsIn.ocean_turning_angle),
-              sin(radians * paramsIn.ocean_turning_angle), paramsIn, avgU, avgV)
+        : CGDynamicsKernel<DGadvection>()
         , stressStep(stressStepIn)
         , params(reinterpret_cast<const MEBParameters&>(paramsIn))
         , stresstransport(nullptr)
@@ -81,6 +80,9 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
         damage.resize_by_mesh(*smesh);
         avgU.resize_by_mesh(*smesh);
         avgV.resize_by_mesh(*smesh);
+
+        cosOceanAngle = cos(radians * params.ocean_turning_angle);
+        sinOceanAngle = sin(radians * params.ocean_turning_angle);
     }
 
     // The brittle rheologies use avgU and avgV to do the advection, not u and v, like mEVP
@@ -163,12 +165,38 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
         }
     }
 
+    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    {
+        CGVector<CGdegree> taux, tauy;
+        taux.resizeLike(avgU);
+        tauy.resizeLike(avgV);
+
+#pragma omp parallel for
+        for (int i = 0; i < taux.rows(); ++i) {
+            const double uOceanRel = uOcean(i) - avgU(i);
+            const double vOceanRel = vOcean(i) - avgV(i);
+            const double cPrime = params.F_ocean * std::hypot(uOceanRel, vOceanRel);
+
+            taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+            tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        }
+
+        if (name == uIOStressName) {
+            return taux;
+        } else if (name == vIOStressName) {
+            return tauy;
+        } else {
+            throw std::logic_error(std::string(__func__) + " called with an unknown argument "
+                + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n");
+        }
+    }
+
 protected:
     CGVector<CGdegree> avgU;
     CGVector<CGdegree> avgV;
 
     StressUpdateStep<DGadvection, DGstressComp>& stressStep;
-    const MEBParameters params;
+    const MEBParameters& params;
 
     Nextsim::DGTransport<DGstressComp>* stresstransport;
 

dynamics/src/include/CGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file CGDynamicsKernel.hpp
  *
- * @date 07 Oct 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -35,22 +35,10 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
     using DynamicsKernel<DGadvection, DGstressComp>::smesh;
     using DynamicsKernel<DGadvection, DGstressComp>::dgtransport;
     using typename DynamicsKernel<DGadvection, DGstressComp>::DataMap;
-    using DynamicsKernel<DGadvection, DGstressComp>::cosOceanAngle;
-    using DynamicsKernel<DGadvection, DGstressComp>::sinOceanAngle;
 
 public:
-    CGDynamicsKernel(double cosOceanAngleIn, double sinOceanAngleIn,
-        const DynamicsParameters& paramsIn, CGVector<CGdegree>& uStressRef,
-        CGVector<CGdegree>& vStressRef)
-        : DynamicsKernel<DGadvection, DGstressComp>(cosOceanAngleIn, sinOceanAngleIn, paramsIn)
-        , pmap(nullptr)
-        , uStress(uStressRef)
-        , vStress(vStressRef)
-    {
-    }
-    CGDynamicsKernel(
-        double cosOceanAngleIn, const DynamicsParameters& paramsIn, double sinOceanAngleIn)
-        : CGDynamicsKernel(cosOceanAngleIn, sinOceanAngleIn, paramsIn, u, v)
+    CGDynamicsKernel()
+        : pmap(nullptr)
     {
     }
     virtual ~CGDynamicsKernel() = default;
@@ -65,7 +53,7 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
     void applyBoundaries() override;
     void prepareAdvection() override;
 
-    virtual std::pair<CGVector<CGdegree>, CGVector<CGdegree>> getIceOceanStress() const;
+    virtual CGVector<CGdegree> getIceOceanStress(const std::string& name) const = 0;
 
 protected:
     void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy);
@@ -94,10 +82,6 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
     CGVector<CGdegree> uAtmos;
     CGVector<CGdegree> vAtmos;
 
-    // Velocities used to calculate ice-ocean stress. The references must be set in derived classes.
-    CGVector<CGdegree>& uStress;
-    CGVector<CGdegree>& vStress;
-
     ParametricMomentumMap<CGdegree, DGadvection>* pmap;
 };
 

dynamics/src/include/DynamicsKernel.hpp

@@ -1,15 +1,14 @@
 /*!
  * @file DynamicsKernel.hpp
  *
- * @date 07 Oct 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
 #ifndef DYNAMICSKERNEL_HPP
 #define DYNAMICSKERNEL_HPP
 
 #include "DGTransport.hpp"
-#include "DynamicsParameters.hpp"
 #include "Interpolations.hpp"
 #include "ParametricMesh.hpp"
 #include "ParametricTools.hpp"
@@ -40,14 +39,7 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
     typedef std::pair<const std::string, const DGVector<DGadvection>&> DataMapping;
     typedef std::map<typename DataMapping::first_type, typename DataMapping::second_type> DataMap;
 
-    DynamicsKernel(
-        double cosOceanAngleIn, double sinOceanAngleIn, const DynamicsParameters& paramsIn)
-        : cosOceanAngle(cosOceanAngleIn)
-        , sinOceanAngle(sinOceanAngleIn)
-        , baseParams(paramsIn)
-        , deltaT(0.)
-    {
-    }
+    DynamicsKernel() = default;
     virtual ~DynamicsKernel() = default;
     virtual void initialise(const ModelArray& coords, bool isSpherical, const ModelArray& mask)
     {
@@ -200,10 +192,6 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
 
     Nextsim::ParametricMesh* smesh;
 
-    // Components of the ocean turning angle
-    const double cosOceanAngle = 1.;
-    const double sinOceanAngle = 0.;
-
     virtual void updateMomentum(const TimestepTime& tst) = 0;
 
     // Pass through functions to the common momentum solver class
@@ -233,18 +221,11 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
      */
     virtual void prepareAdvection() = 0;
 
-    /*!
-     * Returns a const reference to the dynamics parameters
-     */
-    const DynamicsParameters& getParams() const { return baseParams; }
-
 private:
     std::unordered_map<std::string, DGVector<DGadvection>> advectedFields;
 
     // A map from field name to the type of
     std::unordered_map<std::string, ModelArray::Type> fieldType;
-
-    const DynamicsParameters& baseParams;
 };
 
 }

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 07 Oct 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -22,8 +22,6 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
     using DynamicsKernel<DGadvection, DGstressComp>::cice;
     using DynamicsKernel<DGadvection, DGstressComp>::advectionAndLimits;
     using DynamicsKernel<DGadvection, DGstressComp>::dgtransport;
-    using DynamicsKernel<DGadvection, DGstressComp>::cosOceanAngle;
-    using DynamicsKernel<DGadvection, DGstressComp>::sinOceanAngle;
 
     using CGDynamicsKernel<DGadvection>::u;
     using CGDynamicsKernel<DGadvection>::v;
@@ -35,8 +33,7 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
 
 public:
     FreeDriftDynamicsKernel(const DynamicsParameters& paramsIn)
-        : CGDynamicsKernel<DGadvection>(cos(radians * paramsIn.ocean_turning_angle),
-              sin(radians * paramsIn.ocean_turning_angle), paramsIn, u, v)
+        : CGDynamicsKernel<DGadvection>()
         , params(paramsIn)
     {
     }
@@ -53,8 +50,10 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
     };
 
 protected:
-    const DynamicsParameters params;
+    const DynamicsParameters& params;
 
+    const double cosOceanAngle = cos(radians * params.ocean_turning_angle);
+    const double sinOceanAngle = sin(radians * params.ocean_turning_angle);
     const double NansenNumber = sqrt(params.F_atm / params.F_ocean);
 
     void updateMomentum(const TimestepTime& tst) override
@@ -68,6 +67,32 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
                 + NansenNumber * (-uAtmos(i) * sinOceanAngle + vAtmos(i) * cosOceanAngle);
         }
     }
+
+    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    {
+        CGVector<CGdegree> taux, tauy;
+        taux.resizeLike(u);
+        tauy.resizeLike(v);
+
+#pragma omp parallel for
+        for (int i = 0; i < taux.rows(); ++i) {
+            const double uOceanRel = uOcean(i) - u(i);
+            const double vOceanRel = vOcean(i) - v(i);
+            const double cPrime = params.F_ocean * std::hypot(uOceanRel, vOceanRel);
+
+            taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+            tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        }
+
+        if (name == uIOStressName) {
+            return taux;
+        } else if (name == vIOStressName) {
+            return tauy;
+        } else {
+            throw std::logic_error(std::string(__func__) + " called with an unknown argument "
+                + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n");
+        }
+    }
 };
 } /* namespace Nextsim */
 

dynamics/src/include/VPCGDynamicsKernel.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file VPCGDynamicsKernel.hpp
  *
- * @date Aug 23, 2024
+ * @date 22 Oct 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -54,13 +54,11 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG
 public:
     VPCGDynamicsKernel(StressUpdateStep<DGadvection, DGstressComp>& stressStepIn,
         const DynamicsParameters& paramsIn)
-        // Note that the ocean turning angle is explicitly zero at present
-        : CGDynamicsKernel<DGadvection>(1., 0., paramsIn, u, v)
+        : CGDynamicsKernel<DGadvection>()
         , stressStep(stressStepIn)
         , params(reinterpret_cast<const VPParameters&>(paramsIn))
     {
     }
-
     virtual ~VPCGDynamicsKernel() = default;
     void update(const TimestepTime& tst) override
     {
@@ -95,9 +93,35 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG
         DynamicsKernel<DGadvection, DGstressComp>::update(tst);
     }
 
+    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    {
+        CGVector<CGdegree> taux, tauy;
+        taux.resizeLike(u);
+        tauy.resizeLike(v);
+
+#pragma omp parallel for
+        for (int i = 0; i < taux.rows(); ++i) {
+            double uOcnRel = u(i) - uOcean(i);
+            double vOcnRel = v(i) - vOcean(i);
+            double absocn = sqrt(SQR(uOcnRel) + SQR(vOcnRel));
+
+            taux(i) = params.F_ocean * absocn * uOcnRel;
+            tauy(i) = params.F_ocean * absocn * vOcnRel;
+        }
+
+        if (name == uIOStressName) {
+            return taux;
+        } else if (name == vIOStressName) {
+            return tauy;
+        } else {
+            throw std::logic_error(std::string(__func__) + " called with an unknown argument "
+                + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n");
+        }
+    }
+
 protected:
     StressUpdateStep<DGadvection, DGstressComp>& stressStep;
-    const VPParameters params;
+    const VPParameters& params;
     const double alpha = 1500.;
     const double beta = 1500.;