Fixes to getIceOceanStress functions

They didn't work properly because of changes in the structure of params.
nextsimhub · Dec 6, 2024 · 9ee8c56 · 9ee8c56
1 parent 62c9395
commit 9ee8c56
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 9 deletions.
diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/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>
  */
@@ -169,11 +169,13 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
         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 = params.F_ocean * std::hypot(uOceanRel, vOceanRel);
+            const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel);
 
             taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
             tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);

diff --git a/dynamics/src/include/FreeDriftDynamicsKernel.hpp b/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>
  */
@@ -76,11 +76,13 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
         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 = params.F_ocean * std::hypot(uOceanRel, vOceanRel);
+            const double cPrime = FOcean * std::hypot(uOceanRel, vOceanRel);
 
             taux(i) = cPrime * (uOceanRel * cosOceanAngle - vOceanRel * sinOceanAngle);
             tauy(i) = cPrime * (vOceanRel * cosOceanAngle + uOceanRel * sinOceanAngle);
@@ -91,8 +93,8 @@ template <int DGadvection> class FreeDriftDynamicsKernel : public CGDynamicsKern
         } 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");
+            throw std::logic_error(std::string(__func__) + " called with an unknown argument "
+                + name + ". Only " + uIOStressName + " and " + vIOStressName + " are supported\n");
         }
     }
 };

diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/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>
  */
 
@@ -95,14 +95,16 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG
         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));
 
-            taux(i) = params.F_ocean * absocn * uOcnRel;
-            tauy(i) = params.F_ocean * absocn * vOcnRel;
+            taux(i) = FOcean * absocn * uOcnRel;
+            tauy(i) = FOcean * absocn * vOcnRel;
         }
 
         if (name == uIOStressName) {