Skip to content

Issue 878 CG ice velocity components #879

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 14 commits into
base: develop
Choose a base branch
from
Open

Issue 878 CG ice velocity components #879

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
b230a73
Handle different ModelArray types in CGDynamicsKernel::ma2cg.
timspainNERSC Jun 23, 2025
3a88d27
Add the getCGData function.
timspainNERSC Jun 23, 2025
d309245
Use getCGData to retrieve full CG ice velocity components.
timspainNERSC Jun 23, 2025
8ad8702
Remove commented out code.
timspainNERSC Jun 23, 2025
21049d0
clang formatting
timspainNERSC Jun 23, 2025
4262ce6
Interpolate from DG0 to CG, if necessary.
timspainNERSC Jul 28, 2025
97a2a89
Read CG ice velocity data. Set the CG dimensions from the dynamics.
timspainNERSC Jul 28, 2025
90eb8de
Use the correct length ofr dimensions (to get CG dimensions to work).
timspainNERSC Jul 28, 2025
c168d80
Correct dimension names for CG fields.
timspainNERSC Jul 28, 2025
0cd2868
clang formatting
timspainNERSC Jul 28, 2025
6d846b8
Remove the U ModelArray type.
timspainNERSC Jul 28, 2025
ca3dddb
Merge remote-tracking branch 'origin/develop' into issue878_cguv
einola Sep 8, 2025
5862887
Merge branch 'develop' into issue878_cguv
TomMelt Oct 21, 2025
7fb04dd
Merge branch 'develop' into issue878_cguv
TomMelt Nov 3, 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
4 changes: 2 additions & 2 deletions core/src/ParaGridIO.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ ParaGridIO::ParaGridIO(ParametricGrid& grid)
{ "ydimxdimdg_comp", ModelArray::Type::DG },
{ "yxdgstress_comp", ModelArray::Type::DGSTRESS },
{ "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
{ "ycgxcg", ModelArray::Type::CG },
{ "y_cgx_cg", ModelArray::Type::CG },
{ "yvertexxvertexncoords", ModelArray::Type::VERTEX },
// clang-format on
})
Expand Down Expand Up @@ -318,7 +318,7 @@ void ParaGridIO::dumpModelState(const ModelState& state, const std::string& file
for (auto entry : ModelArray::definedDimensions) {
ModelArray::Dimension dim = entry.first;
size_t dimSz = (dimCompMap.count(dim)) ? ModelArray::nComponents(dimCompMap.at(dim))
: dimSz = entry.second.globalLength;
: dimSz = entry.second.localLength;
Copy link
Contributor

@TomMelt TomMelt Jul 29, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

found it 😓

localLength should be globalLength. This part sets the total dimension of the array ready to be written to file. Only the portion of ModelArray will actually be written, but netcdf parallel needs to know the total size so it can write to the correct portion.

Suggested change
: dimSz = entry.second.localLength;
: dimSz = entry.second.globalLength;

Copy link
Contributor

@TomMelt TomMelt Jul 29, 2025

Choose a reason for hiding this comment

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

I tested locally and this fixes the issue with the test

Copy link
Collaborator Author

@timspainNERSC timspainNERSC Aug 4, 2025

Choose a reason for hiding this comment

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

I made this change because a 1 day example run was failing with the following error message. Making the change to localLength let the file writing finish correctly.

2010-01-01T23:45:00Z
libc++abi: terminating due to uncaught exception of type netCDF::exceptions::NcEdge: NetCDF: Start+count exceeds dimension bound
file: /tmp/netcdf-cxx-20250303-65682-ahp603/netcdf-cxx4-4.3.1/cxx4/ncVar.cpp  line:1049
Abort trap: 6

ncFromMAMap[dim] = ncFile.addDim(entry.second.name, dimSz);
// TODO Do I need to add data, even if it is just integers 0...n-1?
}
Expand Down
18 changes: 2 additions & 16 deletions core/src/discontinuousgalerkin/ModelArrayDetails.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,16 +51,6 @@ ModelArray::TypeDimensions ModelArray::typeDimensions = {
ModelArray::Dimension::XVERTEX,
ModelArray::Dimension::YVERTEX,
} },
{ ModelArray::Type::U,
{
ModelArray::Dimension::X,
ModelArray::Dimension::Y,
} },
{ ModelArray::Type::V,
{
ModelArray::Dimension::X,
ModelArray::Dimension::Y,
} },
{ ModelArray::Type::DG,
{
ModelArray::Dimension::X,
Expand All @@ -81,8 +71,6 @@ ModelArray::TypeDimensions ModelArray::typeDimensions = {
const std::map<ModelArray::Type, std::string> ModelArray::typeNames = {
{ ModelArray::Type::H, "HField" },
{ ModelArray::Type::VERTEX, "VertexField" },
{ ModelArray::Type::U, "UField" },
{ ModelArray::Type::V, "VField" },
{ ModelArray::Type::DG, "DGField" },
{ ModelArray::Type::DGSTRESS, "DGStressField" },
{ ModelArray::Type::CG, "CGField" },
Expand All @@ -107,17 +95,15 @@ bool ModelArray::hasDoF(const Type type)
}

ModelArray::SizeMap::SizeMap()
: m_sizes({ { Type::H, 0 }, { Type::VERTEX, 1 }, { Type::U, 0 }, { Type::V, 0 },
{ Type::DG, 0 }, { Type::DGSTRESS, 0 }, { Type::CG, 1 } })
: m_sizes({ { Type::H, 0 }, { Type::VERTEX, 1 }, { Type::DG, 0 }, { Type::DGSTRESS, 0 },
{ Type::CG, 1 } })
{
}

ModelArray::DimensionMap::DimensionMap()
: m_dimensions({
{ Type::H, { 0, 0 } },
{ Type::VERTEX, { 1, 1 } },
{ Type::U, { 0, 0 } },
{ Type::V, { 0, 0 } },
{ Type::DG, { 0, 0 } },
{ Type::DGSTRESS, { 0, 0 } },
{ Type::CG, { 1, 1 } },
Expand Down
4 changes: 0 additions & 4 deletions core/src/discontinuousgalerkin/include/ModelArrayDetails.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,6 @@ enum class Dimension { X, Y, XVERTEX, YVERTEX, XCG, YCG, DG, DGSTRESS, NCOORDS,
enum class Type {
H,
VERTEX,
U,
V,
DG,
DGSTRESS,
CG,
Expand All @@ -29,8 +27,6 @@ static const Type AdvectionType = Type::DG;

static ModelArray HField() { return ModelArray(Type::H); }
static ModelArray VertexField() { return ModelArray(Type::VERTEX); }
static ModelArray UField() { return ModelArray(Type::U); }
static ModelArray VField() { return ModelArray(Type::V); }
static ModelArray DGField() { return ModelArray(Type::DG); }
static ModelArray DGSField() { return ModelArray(Type::DGSTRESS); }
static ModelArray CGField() { return ModelArray(Type::CG); }
Expand Down
2 changes: 0 additions & 2 deletions core/src/discontinuousgalerkin/include/ModelArrayTypedefs.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,6 @@

typedef ModelArray HField;
typedef ModelArray VertexField;
typedef ModelArray UField;
typedef ModelArray VField;
typedef ModelArray DGField;
typedef ModelArray DGSField;
typedef ModelArray CGField;
2 changes: 0 additions & 2 deletions core/src/finitevolume/include/ModelArrayTypedefs.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,4 @@
// ModelArrayDetails.cpp

typedef ModelArray HField;
typedef ModelArray UField;
typedef ModelArray VField;
typedef ModelArray VertexField;
30 changes: 30 additions & 0 deletions core/src/modules/DynamicsModule/BBMDynamics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,27 @@ void BBMDynamics::setData(const ModelState::DataMap& ms)
uice = ms.at(uName);
vice = ms.at(vName);

// Handle CG data field sizes
const ModelArray::MultiDim cgDims = kernel.getCGDimensions();
/* If the u is Type::CG and the dimensions read from the restart file do
* not match those of the CG dynamics kernel, throw an exception.
* TODO: Support interpolation between different CG number of Gauss points.
Copy link
Contributor

@joewallwork joewallwork Aug 26, 2025

Choose a reason for hiding this comment

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

Is this TODO for this PR or for a future one?

*/
if (uice.getType() == ModelArray::Type::CG && uice.dimensions()[0] != cgDims[0]) {
// CG degree of the read data
const unsigned int fileCGdegree = (ModelArray::dimensions(ModelArray::Type::CG)[0] - 1)
/ ModelArray::dimensions(ModelArray::Type::H)[0];
const unsigned int modelCGdegree = (cgDims[0] - 1)
/ ModelArray::dimensions(ModelArray::Type::H)[0];

throw std::runtime_error(
"Differing CG degrees between input data and model are not supported. File CG degree = "
+ std::to_string(fileCGdegree) + ", model CG degree = "
+ std::to_string(modelCGdegree) + ".");
}
// Set the dimensions of CG arrays
ModelArray::setDimensions(ModelArray::Type::CG, cgDims);

// Set the data in the kernel arrays.
// Required data
for (const auto& fieldName : namedFields) {
Expand Down Expand Up @@ -170,6 +191,15 @@ void BBMDynamics::advectField(
kernel.advectField(timestep, field, lowerLimit, upperLimit);
}

ModelState BBMDynamics::getStatePrognostic() const
{
return { {
{ uName, kernel.getCGData(uName) },
{ vName, kernel.getCGData(vName) },
{ damageName, damage },
}, { getConfiguration() } };
}

BBMDynamics::HelpMap& BBMDynamics::getHelpText(HelpMap& map, bool getAll)
{
map["BBMDynamics"] = {
Expand Down
29 changes: 29 additions & 0 deletions core/src/modules/DynamicsModule/MEVPDynamics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,27 @@ void MEVPDynamics::setData(const ModelState::DataMap& ms)
uice = ms.at(uName);
vice = ms.at(vName);

// Handle CG data field sizes
const ModelArray::MultiDim cgDims = kernel.getCGDimensions();
/* If the u is Type::CG and the dimensions read from the restart file do
* not match those of the CG dynamics kernel, throw an exception.
* TODO: Support interpolation between different CG number of Gauss points.
*/
if (uice.getType() == ModelArray::Type::CG && uice.dimensions()[0] != cgDims[0]) {
// CG degree of the read data
const unsigned int fileCGdegree = (ModelArray::dimensions(ModelArray::Type::CG)[0] - 1)
/ ModelArray::dimensions(ModelArray::Type::H)[0];
const unsigned int modelCGdegree = (cgDims[0] - 1)
/ ModelArray::dimensions(ModelArray::Type::H)[0];

throw std::runtime_error(
"Differing CG degrees between input data and model are not supported. File CG degree = "
+ std::to_string(fileCGdegree) + ", model CG degree = "
+ std::to_string(modelCGdegree) + ".");
}
// Set the dimensions of CG arrays
ModelArray::setDimensions(ModelArray::Type::CG, cgDims);

// Set the data in the kernel arrays.
for (const auto& fieldName : namedFields) {
kernel.setData(fieldName, ms.at(fieldName));
Expand Down Expand Up @@ -131,6 +152,14 @@ void MEVPDynamics::advectField(
kernel.advectField(timestep, field, lowerLimit, upperLimit);
}

ModelState MEVPDynamics::getStatePrognostic() const
{
return { {
{ uName, kernel.getCGData(uName) },
{ vName, kernel.getCGData(vName) },
}, { getConfiguration() } };
}

MEVPDynamics::HelpMap& MEVPDynamics::getHelpText(HelpMap& map, bool getAll)
{
map["MEVPDynamics"] = {
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 @@ -33,6 +33,8 @@ class BBMDynamics : public IDynamics, public Configured<BBMDynamics> {
void configure() override;
ConfigMap getConfiguration() const override;

ModelState getStatePrognostic() const override;

enum {
C_KEY,
NU_KEY,
Expand Down
8 changes: 8 additions & 0 deletions core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,14 @@ class FreeDriftDynamics : public IDynamics {

void prepareAdvection() override { kernel.prepareAdvection(); }

ModelState getStatePrognostic() const override
{
return { {
{ uName, kernel.getCGData(uName) },
{ vName, kernel.getCGData(vName) },
}, { getConfiguration() } };
}

private:
FreeDriftDynamicsKernel<DGCOMP> kernel;
DynamicsParameters params;
Expand Down
2 changes: 2 additions & 0 deletions core/src/modules/DynamicsModule/include/MEVPDynamics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ class MEVPDynamics : public IDynamics, public Configured<MEVPDynamics> {
void configure() override;
ConfigMap getConfiguration() const override;

ModelState getStatePrognostic() const override;

enum {
PSTAR_KEY,
DELTA_KEY,
Expand Down
18 changes: 0 additions & 18 deletions core/src/modules/include/IDynamics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,29 +49,11 @@ class IDynamics : public ModelComponent {
}
virtual ~IDynamics() = default;

ModelState getStatePrognostic() const override
{
ModelState state = { {
{ uName, uice },
{ vName, vice },
},
getConfiguration() };

if (m_usesDamage) {
ModelState::DataMap damageState = { { damageName, damage } };
state.merge(damageState);
}

return state;
}

ModelState getStateDiagnostic() const override
{
ModelState state = { {
{ uIOStressName, taux },
{ vIOStressName, tauy },
{ uName, uice },
{ vName, vice },
{ shearName, shear },
{ divergenceName, divergence },
{ sigmaIName, sigmaI },
Expand Down
27 changes: 23 additions & 4 deletions dynamics/src/CGDynamicsKernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,22 +100,22 @@ template <int DGadvection>
ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) const
{
if (name == uName) {
ModelArray data(ModelArray::Type::U);
ModelArray data(ModelArray::Type::H);
DGVector<DGadvection> utmp(*smesh);
Nextsim::Interpolations::CG2DG(*smesh, utmp, u);
return DGModelArray::dg2ma(utmp, data);
} else if (name == vName) {
ModelArray data(ModelArray::Type::V);
ModelArray data(ModelArray::Type::H);
DGVector<DGadvection> vtmp(*smesh);
Nextsim::Interpolations::CG2DG(*smesh, vtmp, v);
return DGModelArray::dg2ma(vtmp, data);
} else if (name == uIOStressName) {
ModelArray data(ModelArray::Type::U);
ModelArray data(ModelArray::Type::H);
DGVector<DGadvection> utmp(*smesh);
Nextsim::Interpolations::CG2DG(*smesh, utmp, uIceOceanStress);
return DGModelArray::dg2ma(utmp, data);
} else if (name == vIOStressName) {
ModelArray data(ModelArray::Type::V);
ModelArray data(ModelArray::Type::H);
Comment on lines +103 to +118
Copy link
Contributor

@joewallwork joewallwork Aug 26, 2025

Choose a reason for hiding this comment

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

You could now hoist the declaration out of the conditional block to avoid duplication.

DGVector<DGadvection> vtmp(*smesh);
Nextsim::Interpolations::CG2DG(*smesh, vtmp, vIceOceanStress);
return DGModelArray::dg2ma(vtmp, data);
Expand All @@ -124,6 +124,25 @@ ModelArray CGDynamicsKernel<DGadvection>::getDG0Data(const std::string& name) co
}
}

template <int DGadvection>
ModelArray CGDynamicsKernel<DGadvection>::getCGData(const std::string& name) const
{
CGField madata(ModelArray::Type::CG);
madata.resize();
if (name == uName) {
madata = u;
} else if (name == vName) {
madata = v;
}
return madata;
}

template <int DGadvection>
const ModelArray::MultiDim CGDynamicsKernel<DGadvection>::getCGDimensions() const
{
return { CGdegree * smesh->nx + 1, CGdegree * smesh->ny + 1 };
}

template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareAdvection()
{
dgtransport->prepareAdvection(u, v);
Expand Down
23 changes: 19 additions & 4 deletions dynamics/src/include/CGDynamicsKernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,9 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {

void setData(const std::string& name, const ModelArray& data) override;
ModelArray getDG0Data(const std::string& name) const override;
ModelArray getCGData(const std::string& name) const;

const ModelArray::MultiDim getCGDimensions() const;
void computeGradientOfSeaSurfaceHeight(const DGVector<1>& seaSurfaceHeight);
void prepareIteration(const DataMap& data) override;
void projectVelocityToStrain() override;
Expand Down Expand Up @@ -96,10 +99,22 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {

CGVector<CGdegree>& ma2cg(const ModelArray& maData, CGVector<CGdegree>& cgData)
{
DGVector<DGadvection> dgtmp(*smesh);
dgtmp.zero();
DGModelArray::ma2dg(maData, dgtmp);
Nextsim::Interpolations::DG2CG(*smesh, cgData, dgtmp);
if (maData.getType() == ModelArray::Type::H) {
DGVector<DGadvection> dgtmp(*smesh);
dgtmp.zero();
DGModelArray::ma2dg(maData, dgtmp);
Nextsim::Interpolations::DG2CG(*smesh, cgData, dgtmp);
} else if (maData.getType() == ModelArray::Type::DG) {
const DGVector<DGadvection>& asdgv
= reinterpret_cast<const DGVector<DGadvection>&>(maData.data());
Nextsim::Interpolations::DG2CG(*smesh, cgData, asdgv);
} else if (maData.getType() == ModelArray::Type::CG) {
cgData = maData.data().matrix();
} else {
throw std::runtime_error(
std::string("CGDynamicsKernel::ma2cg: unhandled ModelArray type ")
+ ModelArray::typeNames.at(maData.getType()));
}
return cgData;
}
};
Expand Down
20 changes: 17 additions & 3 deletions dynamics/src/include/CGModelArray.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,23 @@ class CGModelArray {
public:
template <int CG> static CGVector<CG>& ma2cg(const ModelArray& ma, CGVector<CG>& cg)
{
cg = ma.data().matrix();
//! Interpolation of DG0 to CGVector<CG>
// Nextsim::Interpolations::DG2CG(smesh, cg, ma);
if (ma.getType() != ModelArray::Type::CG) {
/*
* Create a ParametricMesh with the correct x and y dimensions, the
* only members used by the Interpolations functions. Constructed
* with Cartesian coordinates, but the coordinate system is not used.
*/
ParametricMesh smesh(Nextsim::CARTESIAN); // The coordinate system is unimportant here
smesh.nx = ma.dimensions()[0];
smesh.ny = ma.dimensions()[1];
// Assume the data is compatible with a DG0 array
DGVector<1> asDG(smesh);
asDG = ma.data().matrix();
Interpolations::DG2CG(smesh, cg, asDG);
} else {
// CG to CG. Assume the CG degrees are equal
cg = ma.data().matrix();
}
return cg;
}

Expand Down
1 change: 1 addition & 0 deletions dynamics/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ if(NOT ENABLE_MPI)
"${CoreDir}/ModelArray.cpp"
"${CoreDir}/ModelArraySlice.cpp"
"${CoreDir}/${ModelArrayStructure}/ModelArrayDetails.cpp"
"${SRC_DIR}/Interpolations.cpp"
)
target_include_directories(
testCGModelArray
Expand Down
8 changes: 4 additions & 4 deletions physics/src/modules/include/IAtmosphereBoundary.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ class IAtmosphereBoundary : public CheckingModelComponent {
, snow(ModelArray::Type::H, { 0, 1e-3 })
, rain(ModelArray::Type::H, { 0, 1e-3 })
, evap(ModelArray::Type::H, { -1e-3, 1e-3 })
, uwind(ModelArray::Type::U, { -100, 100 })
, vwind(ModelArray::Type::V, { -100, 100 })
, uwind(ModelArray::Type::H, { -100, 100 })
, vwind(ModelArray::Type::H, { -100, 100 })
, penSW(ModelArray::Type::H)
, tauXOW(ModelArray::Type::H, { -10, 10 })
, tauYOW(ModelArray::Type::H, { -10, 10 })
Expand Down Expand Up @@ -91,8 +91,8 @@ class IAtmosphereBoundary : public CheckingModelComponent {
HField snow;
HField rain;
HField evap;
UField uwind;
VField vwind;
HField uwind;
HField vwind;
HField penSW;
HField tauXOW; // x(east)-ward open ocean stress, Pa
HField tauYOW; // y(north)-ward open ocean stress, Pa
Expand Down
Loading
Loading