diff --git a/core/src/include/gridNames.hpp b/core/src/include/gridNames.hpp index e3313f89a..76b1f423b 100644 --- a/core/src/include/gridNames.hpp +++ b/core/src/include/gridNames.hpp @@ -25,6 +25,9 @@ static const std::string shearName = "shear"; static const std::string divergenceName = "divergence"; static const std::string sigmaIName = "sigmaI"; static const std::string sigmaIIName = "sigmaII"; +static const std::string stress11Name = "stress11"; +static const std::string stress12Name = "stress12"; +static const std::string stress22Name = "stress22"; static const std::string uWindName = "uwind"; static const std::string vWindName = "vwind"; diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp index bffeb7f4a..3c38c319e 100644 --- a/core/src/modules/DynamicsModule/BBMDynamics.cpp +++ b/core/src/modules/DynamicsModule/BBMDynamics.cpp @@ -10,7 +10,11 @@ namespace Nextsim { static const std::vector namedFields = { uName, vName }; -static const std::map> defaultFields = {}; +static const std::map> defaultFields = { +// { stress11Name, { ModelArray::Type::DG, 0. } }, +// { stress12Name, { ModelArray::Type::DG, 0. } }, +// { stress22Name, { ModelArray::Type::DG, 0. } }, +}; // TODO: We should use getName() here, but it isn't static. static const std::string prefix = "BBMDynamics"; // MEVPDynamics::getName(); @@ -103,6 +107,16 @@ void BBMDynamics::setData(const ModelState::DataMap& ms) uice = ms.at(uName); vice = ms.at(vName); + if (ms.count(stress11Name)) { + stress11 = ms.at(stress11Name); + stress12 = ms.at(stress12Name); + stress22 = ms.at(stress22Name); + } else { + stress11 = 0.; + stress12 = 0.; + stress22 = 0.; + } + // Set the data in the kernel arrays. // Required data for (const auto& fieldName : namedFields) { @@ -130,6 +144,9 @@ void BBMDynamics::setData(const ModelState::DataMap& ms) kernel.setDGArray(ciceName, ciceDG.allComponents()); kernel.setDGArray(hsnowName, hsnowDG.allComponents()); kernel.setDGArray(damageName, damage.allComponents()); + kernel.setDGArray(stress11Name, stress11); + kernel.setDGArray(stress12Name, stress12); + kernel.setDGArray(stress22Name, stress22); } void BBMDynamics::update(const TimestepTime& tst) @@ -162,6 +179,19 @@ void BBMDynamics::update(const TimestepTime& tst) sigmaII = kernel.getDG0Data(sigmaIIName); } +ModelState BBMDynamics::getStatePrognostic() const +{ + ModelState state = { { + { stress11Name, stress11 }, + { stress12Name, stress12 }, + { stress22Name, stress22 }, + }, + getConfiguration() + }; + state.merge(IDynamics::getStatePrognostic()); + return state; +} + void BBMDynamics::prepareAdvection() { kernel.prepareAdvection(); } void BBMDynamics::advectField( diff --git a/core/src/modules/DynamicsModule/include/BBMDynamics.hpp b/core/src/modules/DynamicsModule/include/BBMDynamics.hpp index a72705ea4..05defe5de 100644 --- a/core/src/modules/DynamicsModule/include/BBMDynamics.hpp +++ b/core/src/modules/DynamicsModule/include/BBMDynamics.hpp @@ -21,6 +21,8 @@ class BBMDynamics : public IDynamics, public Configured { public: BBMDynamics(); + ModelState getStatePrognostic() const override; + std::string getName() const override { return "BBMDynamics"; } void prepareAdvection() override; void update(const TimestepTime& tst) override; diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp index 2be8c7ffd..7714db2a1 100644 --- a/core/src/modules/include/IDynamics.hpp +++ b/core/src/modules/include/IDynamics.hpp @@ -27,6 +27,9 @@ class IDynamics : public ModelComponent { , divergence(ModelArray::Type::H) , sigmaI(ModelArray::Type::H) , sigmaII(ModelArray::Type::H) + , stress11(ModelArray::Type::DGSTRESS) + , stress12(ModelArray::Type::DGSTRESS) + , stress22(ModelArray::Type::DGSTRESS) , damage(getStore()) , uwind(getStore()) , vwind(getStore()) @@ -76,6 +79,9 @@ class IDynamics : public ModelComponent { { divergenceName, divergence }, { sigmaIName, sigmaI }, { sigmaIIName, sigmaII }, + { stress11Name, stress11 }, + { stress11Name, stress12 }, + { stress11Name, stress22 }, }, {} }; return state.merge(getStatePrognostic()); @@ -90,6 +96,9 @@ class IDynamics : public ModelComponent { divergence.resize(); sigmaI.resize(); sigmaII.resize(); + stress11.resize(); + stress12.resize(); + stress22.resize(); } virtual void update(const TimestepTime& tst) = 0; @@ -119,6 +128,10 @@ class IDynamics : public ModelComponent { HField divergence; HField sigmaI; HField sigmaII; + // Stress components. Diagnostic for some dynamics, prognostic for brittle and others. + DGSField stress11; + DGSField stress12; + DGSField stress22; // References to the DG0 finite volume data arrays ModelArrayRef damage; diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 4aa2e8a6f..f07707599 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -12,6 +12,7 @@ #include "include/constants.hpp" #include "SlopeLimiter.hpp" +#include "include/DGModelArray.hpp" #include "include/Interpolations.hpp" #include "include/ParametricMap.hpp" #include "include/VectorManipulations.hpp" @@ -42,6 +43,10 @@ void CGDynamicsKernel::initialise( u.resize_by_mesh(*smesh); v.resize_by_mesh(*smesh); + e11.resize_by_mesh(*smesh); + e12.resize_by_mesh(*smesh); + e22.resize_by_mesh(*smesh); + cgH.resize_by_mesh(*smesh); cgA.resize_by_mesh(*smesh); @@ -82,6 +87,21 @@ void CGDynamicsKernel::setData(const std::string& name, const Model } } +template +void CGDynamicsKernel::setDGArray( + const std::string& name, ModelArray::DataType& dgData) +{ + if (name == stress11Name) { + s11 = DGVectorHolder(dgData); + } else if (name == stress12Name) { + s12 = DGVectorHolder(dgData); + } else if (name == stress22Name) { + s22 = DGVectorHolder(dgData); + } else { + DynamicsKernel::setDGArray(name, dgData); + } +} + template ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) const { @@ -105,6 +125,26 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co DGVector vtmp(*smesh); Nextsim::Interpolations::CG2DG(*smesh, vtmp, vIceOceanStress); return DGModelArray::dg2ma(vtmp, data); + } else if (name == shearName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::Shear(*smesh, e11, e12, e22), data); + } else if (name == divergenceName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, e11, e12, e22), data); + } else if (name == sigmaIName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma( + Tools::TensorInvI(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), + data); + } else if (name == sigmaIIName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma( + Tools::TensorInvII(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), + data); } else { return DynamicsKernel::getDG0Data(name); } diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 9dedaca46..44bb560ea 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -20,12 +20,6 @@ namespace Nextsim { // The brittle momentum solver for CG velocity fields template class BrittleCGDynamicsKernel : public CGDynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::hice; using DynamicsKernel::cice; using DynamicsKernel::smesh; @@ -37,6 +31,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using CGDynamicsKernel::u; using CGDynamicsKernel::v; + using CGDynamicsKernel::s11; + using CGDynamicsKernel::s12; + using CGDynamicsKernel::s22; + using CGDynamicsKernel::e11; + using CGDynamicsKernel::e12; + using CGDynamicsKernel::e22; using CGDynamicsKernel::xGradSeaSurfaceHeight; using CGDynamicsKernel::yGradSeaSurfaceHeight; using CGDynamicsKernel::uAtmos; @@ -157,6 +157,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern { if (name == damageName) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); + // } else if (name == stress11Name) { + // DGModelArray::ma2dg(data, s11); + // } else if (name == stress12Name) { + // DGModelArray::ma2dg(data, s12); + // } else if (name == stress22Name) { + // DGModelArray::ma2dg(data, s22); } else { CGDynamicsKernel::setData(name, data); } diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index 68e19dc6c..6449dce33 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -27,14 +27,9 @@ namespace Nextsim { template class CGDynamicsKernel : public DynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::smesh; using DynamicsKernel::dgtransport; + using DynamicsKernel::setDGArray; using typename DynamicsKernel::DataMap; public: @@ -43,6 +38,7 @@ class CGDynamicsKernel : public DynamicsKernel { void initialise(const ModelArray& coords, bool isSpherical, const ModelArray& mask) override; void setData(const std::string& name, const ModelArray& data) override; + void setDGArray(const std::string& name, ModelArray::DataType& dgData) override; ModelArray getDG0Data(const std::string& name) const override; void computeGradientOfSeaSurfaceHeight(const DGVector<1>& seaSurfaceHeight); void prepareIteration(const DataMap& data) override; @@ -68,6 +64,10 @@ class CGDynamicsKernel : public DynamicsKernel { CGVector cgA; CGVector cgH; + //! DGVectorHolders for strain and stress components + DGVector e11, e12, e22; + DGVectorHolder s11, s12, s22; + // CG gradient of the seaSurfaceHeight CGVector xGradSeaSurfaceHeight; CGVector yGradSeaSurfaceHeight; diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index fc228c4e9..8271b4aad 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -57,20 +57,20 @@ template class DynamicsKernel { seaSurfaceHeight.resize_by_mesh(*smesh); - e11.resize_by_mesh(*smesh); - e12.resize_by_mesh(*smesh); - e22.resize_by_mesh(*smesh); - s11.resize_by_mesh(*smesh); - s12.resize_by_mesh(*smesh); - s22.resize_by_mesh(*smesh); + // e11.resize_by_mesh(*smesh); + // e12.resize_by_mesh(*smesh); + // e22.resize_by_mesh(*smesh); + // s11.resize_by_mesh(*smesh); + // s12.resize_by_mesh(*smesh); + // s22.resize_by_mesh(*smesh); // Set initial values to zero. Prognostic fields will be filled from the restart file. - e11.zero(); - e12.zero(); - e22.zero(); - s11.zero(); - s12.zero(); - s22.zero(); + // e11.zero(); + // e12.zero(); + // e22.zero(); + // s11.zero(); + // s12.zero(); + // s22.zero(); } /*! @@ -91,9 +91,16 @@ template class DynamicsKernel { */ virtual void setData(const std::string& name, const ModelArray& data) { - + static const std::set vectorHolderNames = { + hiceName, + ciceName, + hsnowName, + stress11Name, + stress12Name, + stress22Name, + }; // Special cases: hice, cice, (damage, stress) <- not yet implemented - if (name == hiceName || name == ciceName || name == hsnowName) { + if (vectorHolderNames.count(name)) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); } else if (name == sshName) { DGModelArray::ma2dg(data, seaSurfaceHeight); @@ -125,14 +132,6 @@ template class DynamicsKernel { if (name == hiceName || name == ciceName || name == hsnowName) { throw std::runtime_error( std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); - } else if (name == shearName) { - return DGModelArray::dg2ma(Tools::Shear(*smesh, e11, e12, e22), data); - } else if (name == divergenceName) { - return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, e11, e12, e22), data); - } else if (name == sigmaIName) { - return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, s11, s12, s22), data); - } else if (name == sigmaIIName) { - return DGModelArray::dg2ma(Tools::TensorInvII(*smesh, s11, s12, s22), data); } else { // Any other named field must exist return ModelArray(ModelArray::component0Type(ModelArray::AdvectionType)); @@ -199,8 +198,9 @@ template class DynamicsKernel { DGVector<1> seaSurfaceHeight; //! Vectors storing strain and stress components - DGVector e11, e12, e22; - DGVector s11, s12, s22; + // DGVector e11, e12, e22; + // DGVector /*s11,*/ s12, s22; + // DGVectorHolder s11, s12, s22; size_t stepNumber = 0; diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp index 2ac19bcb3..fcdb0e107 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -19,12 +19,6 @@ namespace Nextsim { // The VP pseudo-timestepping momentum equation solver for CG velocities template class VPCGDynamicsKernel : public CGDynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::hice; using DynamicsKernel::cice; using DynamicsKernel::smesh; @@ -35,6 +29,12 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::u; using CGDynamicsKernel::v; + using CGDynamicsKernel::s11; + using CGDynamicsKernel::s12; + using CGDynamicsKernel::s22; + using CGDynamicsKernel::e11; + using CGDynamicsKernel::e12; + using CGDynamicsKernel::e22; using CGDynamicsKernel::xGradSeaSurfaceHeight; using CGDynamicsKernel::yGradSeaSurfaceHeight; using CGDynamicsKernel::uAtmos; diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index 8b87a3160..7b8285ac6 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -8,6 +8,8 @@ #include "dgVector.hpp" #include "include/ModelArray.hpp" +#include + namespace Nextsim { template class DGVectorHolder { @@ -20,6 +22,12 @@ template class DGVectorHolder { DGVectorHolder(ModelArray& ma) : DGVectorHolder(static_cast(ma)) { + if (ma.nComponents() != DG) { + std::string msg = "DGVectorHolder<" + std::to_string(DG) + + ">(ModelArray&): incorrect number of components = " + + std::to_string(ma.nComponents()); + throw std::length_error(msg); + } } DGVectorHolder(ModelArray::DataType& madt) : DGVectorHolder(reinterpret_cast(madt)) @@ -41,6 +49,9 @@ template class DGVectorHolder { void zero() { ref->setZero(); } + auto row(Eigen::Index i) { return static_cast&>(*ref).row(i); } + auto col(Eigen::Index i) { return static_cast&>(*ref).col(i); } + private: EigenDGVector* ref; };