Collect looping logic of getIceOceanStress in base class

CGDynamicsKernel now handles the basic logic of getIceOceanStress and calls getIceOceanStressElement for the physical implementation of each child class. Those are still very similar to each other, but not the same. Maybe one day they will be, and then we need to reconsider this.
nextsimhub · Dec 6, 2024 · 8f65526 · 8f65526
1 parent 9ee8c56
commit 8f65526
Show file tree

Hide file tree

Showing 4 changed files with 47 additions and 68 deletions.
diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp
@@ -163,32 +163,20 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
         }
     }
 
-    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    double getIceOceanStressElement(const std::string& name, const int i) const override
     {
-        CGVector<CGdegree> taux, tauy;
-        taux.resizeLike(avgU);
-        tauy.resizeLike(avgV);
-
         const double FOcean = params.COcean * params.rhoOcean;
 
-#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 = FOcean * std::hypot(uOceanRel, vOceanRel);
-
-            taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
-            tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
-        }
+        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 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");
-        }
+        if (name == uIOStressName)
+            return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+        else if (name == vIOStressName)
+            return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        else
+            return std::nan("");
     }
 
 protected:

diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/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,7 +47,22 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
     void applyBoundaries() override;
     void prepareAdvection() override;
 
-    virtual CGVector<CGdegree> getIceOceanStress(const std::string& name) const = 0;
+    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);

diff --git a/dynamics/src/include/FreeDriftDynamicsKernel.hpp b/dynamics/src/include/FreeDriftDynamicsKernel.hpp
@@ -70,32 +70,20 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
         }
     }
 
-    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    double getIceOceanStressElement(const std::string& name, const int i) const override
     {
-        CGVector<CGdegree> taux, tauy;
-        taux.resizeLike(u);
-        tauy.resizeLike(v);
-
         const double FOcean = params.COcean * params.rhoOcean;
 
-#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 = FOcean * std::hypot(uOceanRel, vOceanRel);
+        const double uOceanRel = uOcean(i) - u(i);
+        const double vOceanRel = vOcean(i) - v(i);
+        const double cPrime = FOcean * 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");
-        }
+        if (name == uIOStressName)
+            return cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
+        else if (name == vIOStressName)
+            return cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
+        else
+            return std::nan("");
     }
 };
 } /* namespace Nextsim */

diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp
@@ -89,32 +89,20 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG
         DynamicsKernel<DGadvection, DGstressComp>::update(tst);
     }
 
-    CGVector<CGdegree> getIceOceanStress(const std::string& name) const override
+    double getIceOceanStressElement(const std::string& name, const int i) const override
     {
-        CGVector<CGdegree> taux, tauy;
-        taux.resizeLike(u);
-        tauy.resizeLike(v);
-
         const double FOcean = params.COcean * params.rhoOcean;
 
-#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));
+        double uOcnRel = u(i) - uOcean(i);
+        double vOcnRel = v(i) - vOcean(i);
+        double absocn = sqrt(SQR(uOcnRel) + SQR(vOcnRel));
 
-            taux(i) = FOcean * absocn * uOcnRel;
-            tauy(i) = FOcean * 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");
-        }
+        if (name == uIOStressName)
+            return FOcean * absocn * uOcnRel;
+        else if (name == vIOStressName)
+            return FOcean * absocn * vOcnRel;
+        else
+            return std::nan("");
     }
 
 protected: