Skip to content

Issue 779: stress restart #930

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
timspainNERSC wants to merge 6 commits into
base: develop
Choose a base branch
from
Open

Issue 779: stress restart #930

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5c83aa9
Read stress components from a restart file.
timspainNERSC Aug 11, 2025
809f590
Add exceptions to catch mismatches in the number of components.
timspainNERSC Sep 1, 2025
6a93434
Move strain & stress and change stress to a DGVectorHolder.
timspainNERSC Sep 24, 2025
c76fb3b
clang formatting
timspainNERSC Sep 24, 2025
0707c10
clang formatting
timspainNERSC Sep 24, 2025
365243c
Merge branch 'develop' into issue779_stressrestart
timspainNERSC Sep 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions core/src/include/gridNames.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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";
Expand Down
32 changes: 31 additions & 1 deletion core/src/modules/DynamicsModule/BBMDynamics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,11 @@
namespace Nextsim {

static const std::vector<std::string> namedFields = { uName, vName };
static const std::map<std::string, std::pair<ModelArray::Type, double>> defaultFields = {};
static const std::map<std::string, std::pair<ModelArray::Type, double>> defaultFields = {
// { stress11Name, { ModelArray::Type::DG, 0. } },
// { stress12Name, { ModelArray::Type::DG, 0. } },
// { stress22Name, { ModelArray::Type::DG, 0. } },
};
Comment on lines -13 to +17
Copy link
Member

@einola einola Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this should be reverted?


// TODO: We should use getName() here, but it isn't static.
static const std::string prefix = "BBMDynamics"; // MEVPDynamics::getName();
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand Down
2 changes: 2 additions & 0 deletions core/src/modules/DynamicsModule/include/BBMDynamics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,8 @@ class BBMDynamics : public IDynamics, public Configured<BBMDynamics> {
public:
BBMDynamics();

ModelState getStatePrognostic() const override;

std::string getName() const override { return "BBMDynamics"; }
void prepareAdvection() override;
void update(const TimestepTime& tst) override;
Expand Down
13 changes: 13 additions & 0 deletions core/src/modules/include/IDynamics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand Down Expand Up @@ -76,6 +79,9 @@ class IDynamics : public ModelComponent {
{ divergenceName, divergence },
{ sigmaIName, sigmaI },
{ sigmaIIName, sigmaII },
{ stress11Name, stress11 },
{ stress11Name, stress12 },
{ stress11Name, stress22 },
},
{} };
return state.merge(getStatePrognostic());
Expand All @@ -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;
Expand Down Expand Up @@ -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<Shared::DAMAGE, RW> damage;

Expand Down
40 changes: 40 additions & 0 deletions dynamics/src/CGDynamicsKernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -42,6 +43,10 @@ void CGDynamicsKernel<DGadvection>::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);

Expand Down Expand Up @@ -82,6 +87,21 @@ void CGDynamicsKernel<DGadvection>::setData(const std::string& name, const Model
}
}

template <int DGadvection>
void CGDynamicsKernel<DGadvection>::setDGArray(
const std::string& name, ModelArray::DataType& dgData)
{
if (name == stress11Name) {
s11 = DGVectorHolder<DGstressComp>(dgData);
} else if (name == stress12Name) {
s12 = DGVectorHolder<DGstressComp>(dgData);
} else if (name == stress22Name) {
s22 = DGVectorHolder<DGstressComp>(dgData);
} else {
DynamicsKernel<DGadvection, DGstressComp>::setDGArray(name, dgData);
}
}

template <int DGadvection>
ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) const
{
Expand All @@ -105,6 +125,26 @@ ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) co
DGVector<DGadvection> 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<const DGVector<DGstressComp>&>(s11),
static_cast<const DGVector<DGstressComp>&>(s12),
static_cast<const DGVector<DGstressComp>&>(s22)),
data);
} else if (name == sigmaIIName) {
ModelArray data(ModelArray::Type::H);
return DGModelArray::dg2ma(
Tools::TensorInvII(*smesh, static_cast<const DGVector<DGstressComp>&>(s11),
static_cast<const DGVector<DGstressComp>&>(s12),
static_cast<const DGVector<DGstressComp>&>(s22)),
data);
} else {
return DynamicsKernel<DGadvection, DGstressComp>::getDG0Data(name);
}
Expand Down
18 changes: 12 additions & 6 deletions dynamics/src/include/BrittleCGDynamicsKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,6 @@ namespace Nextsim {
// The brittle momentum solver for CG velocity fields
template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKernel<DGadvection> {
protected:
using DynamicsKernel<DGadvection, DGstressComp>::s11;
using DynamicsKernel<DGadvection, DGstressComp>::s12;
using DynamicsKernel<DGadvection, DGstressComp>::s22;
using DynamicsKernel<DGadvection, DGstressComp>::e11;
using DynamicsKernel<DGadvection, DGstressComp>::e12;
using DynamicsKernel<DGadvection, DGstressComp>::e22;
using DynamicsKernel<DGadvection, DGstressComp>::hice;
using DynamicsKernel<DGadvection, DGstressComp>::cice;
using DynamicsKernel<DGadvection, DGstressComp>::smesh;
Expand All @@ -37,6 +31,12 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern

using CGDynamicsKernel<DGadvection>::u;
using CGDynamicsKernel<DGadvection>::v;
using CGDynamicsKernel<DGadvection>::s11;
using CGDynamicsKernel<DGadvection>::s12;
using CGDynamicsKernel<DGadvection>::s22;
using CGDynamicsKernel<DGadvection>::e11;
using CGDynamicsKernel<DGadvection>::e12;
using CGDynamicsKernel<DGadvection>::e22;
using CGDynamicsKernel<DGadvection>::xGradSeaSurfaceHeight;
using CGDynamicsKernel<DGadvection>::yGradSeaSurfaceHeight;
using CGDynamicsKernel<DGadvection>::uAtmos;
Expand Down Expand Up @@ -157,6 +157,12 @@ template <int DGadvection> 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);
Comment on lines +160 to +165
Copy link
Member

@einola einola Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is also not needed

} else {
CGDynamicsKernel<DGadvection>::setData(name, data);
}
Expand Down
12 changes: 6 additions & 6 deletions dynamics/src/include/CGDynamicsKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,14 +27,9 @@ namespace Nextsim {
template <int DGadvection>
class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
protected:
using DynamicsKernel<DGadvection, DGstressComp>::s11;
using DynamicsKernel<DGadvection, DGstressComp>::s12;
using DynamicsKernel<DGadvection, DGstressComp>::s22;
using DynamicsKernel<DGadvection, DGstressComp>::e11;
using DynamicsKernel<DGadvection, DGstressComp>::e12;
using DynamicsKernel<DGadvection, DGstressComp>::e22;
using DynamicsKernel<DGadvection, DGstressComp>::smesh;
using DynamicsKernel<DGadvection, DGstressComp>::dgtransport;
using DynamicsKernel<DGadvection, DGstressComp>::setDGArray;
using typename DynamicsKernel<DGadvection, DGstressComp>::DataMap;

public:
Expand All @@ -43,6 +38,7 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
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;
Expand All @@ -68,6 +64,10 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
CGVector<CGdegree> cgA;
CGVector<CGdegree> cgH;

//! DGVectorHolders for strain and stress components
DGVector<DGstressComp> e11, e12, e22;
DGVectorHolder<DGstressComp> s11, s12, s22;

// CG gradient of the seaSurfaceHeight
CGVector<CGdegree> xGradSeaSurfaceHeight;
CGVector<CGdegree> yGradSeaSurfaceHeight;
Expand Down
48 changes: 24 additions & 24 deletions dynamics/src/include/DynamicsKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,20 +57,20 @@ template <int DGadvection, int DGstress> 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);
Comment on lines +60 to +65
Copy link
Member

@einola einola Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove commented code


// 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();
Comment on lines +68 to +73
Copy link
Member

@einola einola Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove commented code

}

/*!
Expand All @@ -91,9 +91,16 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
*/
virtual void setData(const std::string& name, const ModelArray& data)
{

static const std::set<std::string> 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);
Expand Down Expand Up @@ -125,14 +132,6 @@ template <int DGadvection, int DGstress> 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));
Expand Down Expand Up @@ -199,8 +198,9 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
DGVector<1> seaSurfaceHeight;

//! Vectors storing strain and stress components
DGVector<DGstress> e11, e12, e22;
DGVector<DGstress> s11, s12, s22;
// DGVector<DGstress> e11, e12, e22;
// DGVector<DGstress> /*s11,*/ s12, s22;
// DGVectorHolder<DGstress> s11, s12, s22;
Comment on lines +201 to +203
Copy link
Member

@einola einola Oct 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove commented code


size_t stepNumber = 0;

Expand Down
12 changes: 6 additions & 6 deletions dynamics/src/include/VPCGDynamicsKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,6 @@ namespace Nextsim {
// The VP pseudo-timestepping momentum equation solver for CG velocities
template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DGadvection> {
protected:
using DynamicsKernel<DGadvection, DGstressComp>::s11;
using DynamicsKernel<DGadvection, DGstressComp>::s12;
using DynamicsKernel<DGadvection, DGstressComp>::s22;
using DynamicsKernel<DGadvection, DGstressComp>::e11;
using DynamicsKernel<DGadvection, DGstressComp>::e12;
using DynamicsKernel<DGadvection, DGstressComp>::e22;
using DynamicsKernel<DGadvection, DGstressComp>::hice;
using DynamicsKernel<DGadvection, DGstressComp>::cice;
using DynamicsKernel<DGadvection, DGstressComp>::smesh;
Expand All @@ -35,6 +29,12 @@ template <int DGadvection> class VPCGDynamicsKernel : public CGDynamicsKernel<DG

using CGDynamicsKernel<DGadvection>::u;
using CGDynamicsKernel<DGadvection>::v;
using CGDynamicsKernel<DGadvection>::s11;
using CGDynamicsKernel<DGadvection>::s12;
using CGDynamicsKernel<DGadvection>::s22;
using CGDynamicsKernel<DGadvection>::e11;
using CGDynamicsKernel<DGadvection>::e12;
using CGDynamicsKernel<DGadvection>::e22;
using CGDynamicsKernel<DGadvection>::xGradSeaSurfaceHeight;
using CGDynamicsKernel<DGadvection>::yGradSeaSurfaceHeight;
using CGDynamicsKernel<DGadvection>::uAtmos;
Expand Down
11 changes: 11 additions & 0 deletions dynamics/src/include/dgVectorHolder.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#include "dgVector.hpp"
#include "include/ModelArray.hpp"

#include <exception>

namespace Nextsim {

template <int DG> class DGVectorHolder {
Expand All @@ -20,6 +22,12 @@ template <int DG> class DGVectorHolder {
DGVectorHolder(ModelArray& ma)
: DGVectorHolder(static_cast<ModelArray::DataType&>(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<EigenDGVector&>(madt))
Expand All @@ -41,6 +49,9 @@ template <int DG> class DGVectorHolder {

void zero() { ref->setZero(); }

auto row(Eigen::Index i) { return static_cast<DGVector<DG>&>(*ref).row(i); }
auto col(Eigen::Index i) { return static_cast<DGVector<DG>&>(*ref).col(i); }

private:
EigenDGVector* ref;
};
Expand Down
Loading