From 8f6552618eddca4ea8a4a2c98cfd3e3a003dff06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 6 Dec 2024 06:51:42 +0100 Subject: [PATCH] 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. --- .../src/include/BrittleCGDynamicsKernel.hpp | 32 ++++++------------- dynamics/src/include/CGDynamicsKernel.hpp | 19 +++++++++-- .../src/include/FreeDriftDynamicsKernel.hpp | 32 ++++++------------- dynamics/src/include/VPCGDynamicsKernel.hpp | 32 ++++++------------- 4 files changed, 47 insertions(+), 68 deletions(-) diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index d366dc1c7..bedd45270 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -163,32 +163,20 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern } } - CGVector getIceOceanStress(const std::string& name) const override + double getIceOceanStressElement(const std::string& name, const int i) const override { - CGVector 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 index 37abdf3b4..a3bf5670b 100644 --- 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 */ @@ -47,7 +47,22 @@ class CGDynamicsKernel : public DynamicsKernel { void applyBoundaries() override; void prepareAdvection() override; - virtual CGVector getIceOceanStress(const std::string& name) const = 0; + virtual inline double getIceOceanStressElement(const std::string& name, const int i) const = 0; + CGVector 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 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 index b7972fa09..973f21395 100644 --- a/dynamics/src/include/FreeDriftDynamicsKernel.hpp +++ b/dynamics/src/include/FreeDriftDynamicsKernel.hpp @@ -70,32 +70,20 @@ template class FreeDriftDynamicsKernel : public CGDynamicsKern } } - CGVector getIceOceanStress(const std::string& name) const override + double getIceOceanStressElement(const std::string& name, const int i) const override { - CGVector 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 index 2db4ac630..b7c51d8bc 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -89,32 +89,20 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::update(tst); } - CGVector getIceOceanStress(const std::string& name) const override + double getIceOceanStressElement(const std::string& name, const int i) const override { - CGVector 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: