Inactive layers introduce NaNs

[cbegeman]

For the test case E3SM-Project/polaris#353 where there are inactive layers (MaxLayerCell < NVertLayers - 1), NaNs are introduced from the boundary. Thanks to @brian-oneill for identifying this issue.

[brian-oneill]

These are caused by quantities (specifially tracers and vorticities) that are normalized by the layer thickness. The thickness is 0 in inactive layers, so these quantities become NaNs. Since 0 * NaN = NaN, the edge masks don't catch them and they propagate into the domain as the simulation advances. These should be resolved by moving to hierarchical parallelism and restricting the vertical loops to only active layers.

The NaNs are introduced here:

Omega/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h

Lines 49 to 59 in 43eb763

	for (int KVec = 0; KVec < VecLength; ++KVec) {
	const int K = KStart + KVec;
	const Real InvLayerThickVertex = 1._Real / LayerThickVertex[KVec];
	RelVortVertex(IVertex, K) = RelVortVertexTmp[KVec];
	NormRelVortVertex(IVertex, K) =
	RelVortVertexTmp[KVec] * InvLayerThickVertex;
	NormPlanetVortVertex(IVertex, K) =
	FVertex(IVertex) * InvLayerThickVertex;
	}
	}

Omega/components/omega/src/timeStepping/TimeStepper.cpp

Lines 461 to 469 in 43eb763

	parallelFor(
	"updateTracersByTend", {NTracers, Mesh->NCellsAll, NVertLayers},
	KOKKOS_LAMBDA(int L, int ICell, int K) {
	NextTracers(L, ICell, K) =
	(CurTracers(L, ICell, K) * LayerThick2(ICell, K) +
	CoeffSeconds * TracerTend(L, ICell, K)) /
	LayerThick1(ICell, K);
	});
	}

Omega/components/omega/src/timeStepping/TimeStepper.cpp

Lines 519 to 524 in 43eb763

	parallelFor(
	"finalizeTracersUpdate", {NTracers, Mesh->NCellsAll, NVertLayers},
	KOKKOS_LAMBDA(int L, int ICell, int K) {
	NextTracers(L, ICell, K) /= NextThick(ICell, K);
	});
	}

The vorticity normalizations should be resolved by #314, but not the tracer normalizations since this PR doesn't touch the TimeStepper loops yet. Not sure if we want to add some guards like:

if (LayerThickness == 0) InvLayerThickness = 0

since they won't be necessary once the calculations in inactive layers are skipped entirely.

[cbegeman]

Not sure if we want to add some guards like:
if (LayerThickness == 0) InvLayerThickness = 0
since they won't be necessary once the calculations in inactive layers are skipped entirely.

That's a good question. I'm wondering if we don't want to do this because it may obscure if/when we erroneously access inactive vertical layers.

[mwarusz]

That's a good question. I'm wondering if we don't want to do this because it may obscure if/when we erroneously access inactive vertical layers.

I agree, it would be great if we could avoid accessing inactive layers. I have updated #314 to limit vertical loops in the time stepping code. Hopefully, that PR can fix this issue once it gets merged.

[cbegeman]

@mwarusz and @brian-oneill I just wanted to put on your radar that #314 did not resolve this issue as we had hoped (at least based on the results of the barotropic channel case E3SM-Project/polaris#353, which still has NaNs)