From cbb2871a27c522e6e32a8fa4d0316d65b61a3cef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Tue, 22 Oct 2024 13:26:06 +0200 Subject: [PATCH] Hotfix for issue #720 I just reverted the offending comment and fixed a few formatting issues. But this should be done properly at some point. --- dynamics/src/CGDynamicsKernel.cpp | 28 ++----------- .../src/include/BrittleCGDynamicsKernel.hpp | 40 ++++++++++++++++--- dynamics/src/include/CGDynamicsKernel.hpp | 24 ++--------- dynamics/src/include/DynamicsKernel.hpp | 23 +---------- .../src/include/FreeDriftDynamicsKernel.hpp | 37 ++++++++++++++--- dynamics/src/include/VPCGDynamicsKernel.hpp | 34 +++++++++++++--- 6 files changed, 103 insertions(+), 83 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 94e07483b..2e3ef7509 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.cpp * - * @date 07 Oct 2024 + * @date 22 Oct 2024 * @author Tim Spain * @author Einar ÓLason */ @@ -105,12 +105,12 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co } else if (name == uIOStressName) { ModelArray data(ModelArray::Type::U); DGVector 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 vtmp(*smesh); - Nextsim::Interpolations::CG2DG(*smesh, vtmp, getIceOceanStress().second); + Nextsim::Interpolations::CG2DG(*smesh, vtmp, getIceOceanStress(name)); return DGModelArray::dg2ma(vtmp, data); } else { return DynamicsKernel::getDG0Data(name); @@ -443,28 +443,6 @@ template void CGDynamicsKernel::applyBoundaries() // TODO Periodic boundary conditions. } -template -std::pair, CGVector> -CGDynamicsKernel::getIceOceanStress() const -{ - CGVector 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::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>; diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 78264ade1..5b0b907d6 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file BrittleCGDynamicsKernel.hpp * - * @date 07 Oct 2024 + * @date 22 Oct 2024 * @author Tim Spain * @author Einar Ólason */ @@ -38,8 +38,6 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using DynamicsKernel::applyBoundaries; using DynamicsKernel::advectionAndLimits; using DynamicsKernel::dgtransport; - using DynamicsKernel::cosOceanAngle; - using DynamicsKernel::sinOceanAngle; using CGDynamicsKernel::u; using CGDynamicsKernel::v; @@ -58,11 +56,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using CGDynamicsKernel::dStressY; using CGDynamicsKernel::pmap; + double cosOceanAngle, sinOceanAngle; + public: BrittleCGDynamicsKernel(StressUpdateStep& stressStepIn, const DynamicsParameters& paramsIn) - : CGDynamicsKernel(cos(radians * paramsIn.ocean_turning_angle), - sin(radians * paramsIn.ocean_turning_angle), paramsIn, avgU, avgV) + : CGDynamicsKernel() , stressStep(stressStepIn) , params(reinterpret_cast(paramsIn)) , stresstransport(nullptr) @@ -81,6 +80,9 @@ template 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 class BrittleCGDynamicsKernel : public CGDynamicsKern } } + CGVector getIceOceanStress(const std::string& name) const override + { + CGVector 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 avgU; CGVector avgV; StressUpdateStep& stressStep; - const MEBParameters params; + const MEBParameters& params; Nextsim::DGTransport* stresstransport; diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index 8802bd3ef..cb1b07f58 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file CGDynamicsKernel.hpp * - * @date 07 Oct 2024 + * @date 22 Oct 2024 * @author Tim Spain * @author Einar Ólason */ @@ -35,22 +35,10 @@ class CGDynamicsKernel : public DynamicsKernel { using DynamicsKernel::smesh; using DynamicsKernel::dgtransport; using typename DynamicsKernel::DataMap; - using DynamicsKernel::cosOceanAngle; - using DynamicsKernel::sinOceanAngle; public: - CGDynamicsKernel(double cosOceanAngleIn, double sinOceanAngleIn, - const DynamicsParameters& paramsIn, CGVector& uStressRef, - CGVector& vStressRef) - : DynamicsKernel(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 { void applyBoundaries() override; void prepareAdvection() override; - virtual std::pair, CGVector> getIceOceanStress() const; + virtual CGVector 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 { CGVector uAtmos; CGVector vAtmos; - // Velocities used to calculate ice-ocean stress. The references must be set in derived classes. - CGVector& uStress; - CGVector& vStress; - ParametricMomentumMap* pmap; }; diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 5340f1843..4ae72eabb 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file DynamicsKernel.hpp * - * @date 07 Oct 2024 + * @date 22 Oct 2024 * @author Tim Spain */ @@ -9,7 +9,6 @@ #define DYNAMICSKERNEL_HPP #include "DGTransport.hpp" -#include "DynamicsParameters.hpp" #include "Interpolations.hpp" #include "ParametricMesh.hpp" #include "ParametricTools.hpp" @@ -40,14 +39,7 @@ template class DynamicsKernel { typedef std::pair&> DataMapping; typedef std::map 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 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 class DynamicsKernel { */ virtual void prepareAdvection() = 0; - /*! - * Returns a const reference to the dynamics parameters - */ - const DynamicsParameters& getParams() const { return baseParams; } - private: std::unordered_map> advectedFields; // A map from field name to the type of std::unordered_map fieldType; - - const DynamicsParameters& baseParams; }; } diff --git a/dynamics/src/include/FreeDriftDynamicsKernel.hpp b/dynamics/src/include/FreeDriftDynamicsKernel.hpp index c0cdf1f29..8f8d9bd7b 100644 --- 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 07 Oct 2024 + * @date 22 Oct 2024 * @author Tim Spain * @author Einar Ólason */ @@ -22,8 +22,6 @@ template class FreeDriftDynamicsKernel : public CGDynamicsKern using DynamicsKernel::cice; using DynamicsKernel::advectionAndLimits; using DynamicsKernel::dgtransport; - using DynamicsKernel::cosOceanAngle; - using DynamicsKernel::sinOceanAngle; using CGDynamicsKernel::u; using CGDynamicsKernel::v; @@ -35,8 +33,7 @@ template class FreeDriftDynamicsKernel : public CGDynamicsKern public: FreeDriftDynamicsKernel(const DynamicsParameters& paramsIn) - : CGDynamicsKernel(cos(radians * paramsIn.ocean_turning_angle), - sin(radians * paramsIn.ocean_turning_angle), paramsIn, u, v) + : CGDynamicsKernel() , params(paramsIn) { } @@ -53,8 +50,10 @@ template 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 class FreeDriftDynamicsKernel : public CGDynamicsKern + NansenNumber * (-uAtmos(i) * sinOceanAngle + vAtmos(i) * cosOceanAngle); } } + + CGVector getIceOceanStress(const std::string& name) const override + { + CGVector 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 */ diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp index 5fc9c5d63..e07423377 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -1,7 +1,7 @@ /*! * @file VPCGDynamicsKernel.hpp * - * @date Aug 23, 2024 + * @date 22 Oct 2024 * @author Tim Spain * @author Einar Ólason */ @@ -54,13 +54,11 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel& stressStepIn, const DynamicsParameters& paramsIn) - // Note that the ocean turning angle is explicitly zero at present - : CGDynamicsKernel(1., 0., paramsIn, u, v) + : CGDynamicsKernel() , stressStep(stressStepIn) , params(reinterpret_cast(paramsIn)) { } - virtual ~VPCGDynamicsKernel() = default; void update(const TimestepTime& tst) override { @@ -95,9 +93,35 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::update(tst); } + CGVector getIceOceanStress(const std::string& name) const override + { + CGVector 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& stressStep; - const VPParameters params; + const VPParameters& params; const double alpha = 1500.; const double beta = 1500.;