Issue 878 CG ice velocity components

[timspainNERSC]

CG ice velocity components

Fixes #878

---

Change Description

When passing a ModelArray to the ma2cg() function, the operations that are performed depends on the ModelArray::Type. An HField is converted to a DGVector and then interpolated to a CGVector. A DGField will only be interpolated to a CGVector. A CGField will be directly copied into the CGVector.

Adds a getCGData() function to CGDynamicsKernel which will return the full CG array of the ice velocity components. These are then passed to the data maps collected as part of the getDataPrognostic() call tree.

Set the size of the ModelArray::Dimensions XCGand YCG to values derived from the CGdegree value and the size of a DG0 (Type::H) ModelArray.

Since the ice velocities are so completely removed from the thermodynamics/physics part of the model (that is where ModelArray is used), remove the special ModelArray types for U and V.

---

Test Description

1. Run the model from init_25km_NH.nc for 1 TS. Let the resulting restart file be 'file 1'.
2. Run the model from init_25km_NH.nc for 1 day. Let the resulting restart file be 'file 2'.
3. Run the model from 'file 2' for 1 TS. Let the resulting restart file be 'file 3'.
4. Compare the values of u and v in 'file 3' with the values in 'file 1' and 'file 2'. Note the velocity in 'file 3' much more closely matches that of 'file 2' than that of 'file 1' which is only one timestep from being initialized as all zero.

[TomMelt]

I had a look at the failing test and it looks like one of the MPI tests gets stuck in an infinite loop. I ran it locally and I get the same behaviour.

I will see if I can find out why

[TomMelt]

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;

[TomMelt]

I tested locally and this fixes the issue with the test

[timspainNERSC]

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

[joewallwork]

This looks good to me. I just have a query about TODO notes and a minor suggested refactoring.

[joewallwork]

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

[joewallwork]

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

[TomMelt]

@timspainNERSC I have had a more detailed look and I think the underlying problem is that the x_cg and y_cg dimensions are wrong in the initial model file. For example, when I run the config_june23.cfg it uses input file init_25km_NH.nc. In there the dimensions are listed as follows:

netcdf init_25km_NH {
dimensions:
	ydim = 121 ;
	xdim = 154 ;
	yvertex = 122 ;
	xvertex = 155 ;
	y_cg = 122 ;
	x_cg = 155 ;
	dg_comp = 1 ;
	dgstress_comp = 3 ;
	ncoords = 2 ;

They should be y_cg = 243 and x_cg = 309. If this the intention is not to hardcode these to a specific value because we dont want to make separate input files for each CGDegree, then we should just remove these completely and they should be set at run time using (CGDegree * Dim + 1).

These get read in here in ParaGridIO::getModelState and then set on line 137.

nextsimdg/core/src/ParaGridIO.cpp

Lines 81 to 137 in ca3dddb

	ModelState ParaGridIO::getModelState(const std::string& filePath)
	#endif
	{
	ModelState state;
	try {
	#ifdef USE_MPI
	netCDF::NcFilePar ncFile(filePath, netCDF::NcFile::read, metadata.mpiComm);
	#else
	netCDF::NcFile ncFile(filePath, netCDF::NcFile::read);
	#endif
	// Dimensions and DG components
	std::multimap<std::string, netCDF::NcDim> dimMap = ncFile.getDims();
	for (auto entry : ModelArray::definedDimensions) {
	auto dimType = entry.first;
	if (dimCompMap.count(dimType) > 0)
	// TODO Assertions that DG in the file equals the compile time DG in the model. See
	// #205
	continue;
	ModelArray::DimensionSpec& dimensionSpec = entry.second;
	// Find dimensions in the netCDF file by their name in the ModelArray details
	netCDF::NcDim dim = ncFile.getDim(dimensionSpec.name);
	// Also check the old name
	if (dim.isNull()) {
	dim = ncFile.getDim(dimensionSpec.altName);
	}
	// If we didn't find a dimension with the dimensions name or altName, throw.
	if (dim.isNull()) {
	throw std::out_of_range(
	std::string("No netCDF dimension found corresponding to the dimension named ")
	+ dimensionSpec.name + std::string(" or ") + dimensionSpec.altName);
	}
	#ifdef USE_MPI
	auto dimName = dim.getName();
	size_t localLength = 0;
	size_t start = 0;
	if (dimType == ModelArray::Dimension::X) {
	localLength = metadata.localExtentX;
	start = metadata.localCornerX;
	} else if (dimType == ModelArray::Dimension::Y) {
	localLength = metadata.localExtentY;
	start = metadata.localCornerY;
	} else if (dimType == ModelArray::Dimension::XVERTEX) {
	localLength = metadata.localExtentX + 1;
	start = metadata.localCornerX;
	} else if (dimType == ModelArray::Dimension::YVERTEX) {
	localLength = metadata.localExtentY + 1;
	start = metadata.localCornerY;
	} else {
	localLength = dim.getSize();
	start = 0;
	}
	ModelArray::setDimension(dimType, dim.getSize(), localLength, start);
	#else
	ModelArray::setDimension(dimType, dim.getSize());

The localLength is then overridden later in BBMDynamics.cpp on line 125:

nextsimdg/core/src/modules/DynamicsModule/BBMDynamics.cpp

Line 125 in ca3dddb

ModelArray::setDimensions(ModelArray::Type::CG, cgDims);

Note the function ModelArray::setDimensions only sets the localLength and not the globalLength and that's why you needed to make the changes here:

nextsimdg/core/src/ParaGridIO.cpp

Lines 266 to 267 in ca3dddb

	size_t dimSz = (dimCompMap.count(dim)) ? ModelArray::nComponents(dimCompMap.at(dim))
	: dimSz = entry.second.localLength;

I think this logic needs to be re-written to use the ModelMetadata singleton so that it will work for MPI and serial code. Then we need to also modify function ModelArray::setDimensions or create a new one that will also set the globalLength.

Let me know what you think. I can try and implement a fix if you want.