Skip to content

Fix a bug in high-order tracer horizontal advection #403

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
hyungyukang wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Fix a bug in high-order tracer horizontal advection #403

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6c9e8f3
Add boundary masks
hyungyukang May 7, 2026
70c9857
Fix a bug in high-order tracer horizontal advection
hyungyukang May 9, 2026
61ceaac
Set the order option before initialization in TendencyTermsTest.cpp.
hyungyukang May 9, 2026
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
7 changes: 4 additions & 3 deletions components/omega/src/ocn/HorzOperators.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "HorzOperators.h"
#include "DataTypes.h"
#include "HorzMesh.h"
#include "VertCoord.h"

namespace OMEGA {

Expand Down Expand Up @@ -54,9 +55,9 @@ SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh)
}

MasksAndCoefficients::MasksAndCoefficients(
HorzMesh const *Mesh, const Array3DReal DerivTwo,
HorzMesh const *Mesh, VertCoord const *VCoord, const Array3DReal DerivTwo,
Array1DI4 NAdvCellsForEdge, Array2DI4 AdvCellsForEdge,
Array1DI4 AdvMaskHighOrder, Array2DReal AdvCoefs, Array2DReal AdvCoefs3rd)
Array2DI4 AdvMaskHighOrder, Array2DReal AdvCoefs, Array2DReal AdvCoefs3rd)
: NCellsGlobal(Mesh->NCellsGlobal), NCellsAll(Mesh->NCellsAll),
NAdvCellsMax(Mesh->MaxEdges2), // PatchCellLists("PatchCellLists",
// Mesh->NEdgesOwned, Mesh->NEdgesAll+1),
Expand All @@ -68,5 +69,5 @@ MasksAndCoefficients::MasksAndCoefficients(
EdgesOnEdge(Mesh->EdgesOnEdge), CellsOnCell(Mesh->CellsOnCell),
CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge),
DvEdge(Mesh->DvEdge), AdvCoefs(AdvCoefs), AdvCoefs3rd(AdvCoefs3rd),
DerivTwo(DerivTwo) {}
BoundaryCell(VCoord->BoundaryCell), DerivTwo(DerivTwo) {}
} // namespace OMEGA
46 changes: 33 additions & 13 deletions components/omega/src/ocn/HorzOperators.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "GlobalConstants.h"
#include "HorzMesh.h"
#include "HorzUtil.h"
#include "VertCoord.h"

namespace OMEGA {

Expand Down Expand Up @@ -467,10 +468,10 @@ class MasksAndCoefficients {
}

public:
MasksAndCoefficients(HorzMesh const *Mesh, const Array3DReal DerivTwo,
Array1DI4 NAdvCellsForEdge, Array2DI4 AdvCellsForEdge,
Array1DI4 AdvMaskHighOrder, Array2DReal AdvCoefs,
Array2DReal AdvCoefs3rd);
MasksAndCoefficients(HorzMesh const *Mesh, VertCoord const *VCoord,
const Array3DReal DerivTwo, Array1DI4 NAdvCellsForEdge,
Array2DI4 AdvCellsForEdge, Array2DI4 AdvMaskHighOrder,
Array2DReal AdvCoefs, Array2DReal AdvCoefs3rd);

KOKKOS_FUNCTION void operator()(const int IEdge) const {

Expand All @@ -484,17 +485,35 @@ class MasksAndCoefficients {
for (unsigned I = 0; I < CellIndexSorted.extent(0); ++I)
for (unsigned K = 0; K < CellIndexSorted.extent(1); ++K)
CellIndexSorted(I, K) = MaxI4;

const int Cell1 = CellsOnEdge(IEdge, 0);
const int Cell2 = CellsOnEdge(IEdge, 1);

// at boundaries, must stay at low order
AdvMaskHighOrder(IEdge) = 1;
for (int K = 0; K < NEdgesOnCell(Cell1); ++K)
if (CellsOnCell(Cell1, K) == NCellsGlobal)
AdvMaskHighOrder(IEdge) = 0;

for (int K = 0; K < NEdgesOnCell(Cell2); ++K)
if (CellsOnCell(Cell2, K) == NCellsGlobal)
AdvMaskHighOrder(IEdge) = 0;
const int NLyr = AdvMaskHighOrder.extent(1);
for (int K = 0; K < NLyr; ++K) {
if (BoundaryCell(Cell1, K) == 1 or BoundaryCell(Cell2, K) == 1) {
AdvMaskHighOrder(IEdge, K) = 0;
} else {
AdvMaskHighOrder(IEdge, K) = 1;
}
}

for (int J = 0; J < NEdgesOnCell(Cell1); ++J) {
if (CellsOnCell(Cell1, J) == NCellsGlobal) {
for (int K = 0; K < NLyr; ++K) {
AdvMaskHighOrder(IEdge, K) = 0;
}
}
}
for (int J = 0; J < NEdgesOnCell(Cell2); ++J) {
if (CellsOnCell(Cell2, J) == NCellsGlobal) {
for (int K = 0; K < NLyr; ++K) {
AdvMaskHighOrder(IEdge, K) = 0;
}
}
}

// do only if this edge flux is needed to update owned cells
if (Cell1 < NCellsAll && Cell2 < NCellsAll) {
// Insert cellsOnEdge to list of advection cells
Expand Down Expand Up @@ -622,14 +641,15 @@ class MasksAndCoefficients {
Array2DI4 CellIndx;
Array3DI4 CellIndxSorted;
Array1DI4 CellID;
Array1DI4 AdvMaskHighOrder;
Array2DI4 AdvMaskHighOrder;
Array2DI4 EdgesOnEdge;
Array2DI4 CellsOnCell;
Array2DI4 CellsOnEdge;
Array1DReal DcEdge;
Array1DReal DvEdge;
Array2DReal AdvCoefs;
Array2DReal AdvCoefs3rd;
Array2DReal BoundaryCell;
Array3DReal DerivTwo;
};
} // namespace OMEGA
Expand Down
27 changes: 18 additions & 9 deletions components/omega/src/ocn/TendencyTerms.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,12 +73,13 @@ BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh,

TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh,
const VertCoord *VCoord)
: HorzontalMesh(Mesh),
: HorzontalMesh(Mesh), VerticalCoord(VCoord),
NAdvCellsForEdge("NumberOfCellsContribToAdvectionAtEdge",
Mesh->NEdgesAll),
AdvCellsForEdge("IndexOfCellsContributingToAdvection", Mesh->NEdgesAll,
Mesh->MaxEdges2 + 2),
AdvMaskHighOrder("MaskForHighOrderAdvectionTerms", Mesh->NEdgesAll),
AdvMaskHighOrder("MaskForHighOrderAdvectionTerms", Mesh->NEdgesAll,
VCoord->NVertLayers),
AdvCoefs("CommonAdvectionCoefficients", Mesh->MaxEdges2 + 2,
Mesh->NEdgesAll),
AdvCoefs3rd("CommonAdvectionCoeffsForHighOrder", Mesh->MaxEdges2 + 2,
Expand Down Expand Up @@ -115,22 +116,30 @@ SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell(
const HorzMesh *Mesh) {}

void TracerHorzAdvOnCell::init() {
const HorzMesh *Mesh = this->HorzontalMesh;
const auto MaxEdges2 = Mesh->MaxEdges2;
const auto NEdgesAll = Mesh->NEdgesAll;
const auto NCellsAll = Mesh->NCellsAll;
const HorzMesh *Mesh = this->HorzontalMesh;
const VertCoord *VCoord = this->VerticalCoord;
const auto MaxEdges2 = Mesh->MaxEdges2;
const auto NEdgesAll = Mesh->NEdgesAll;
const auto NCellsAll = Mesh->NCellsAll;
// Allocate Kokkos arrays in member data

if (ForceLowOrder) {
// Return when the 2nd-order tracer horz adv
deepCopy(NAdvCellsForEdge, 0);
deepCopy(AdvMaskHighOrder, 0);
return;
}

SecondDerivativeOnCell secondDerivativeOnCell(Mesh);
Array3DReal DerivTwo("DerivTwo", MaxEdges2 + 2, 2, NEdgesAll);
parallelFor(
{NCellsAll},
KOKKOS_LAMBDA(int ICell) { secondDerivativeOnCell(DerivTwo, ICell); });
// Compute masks and coefficients
Kokkos::fence();
MasksAndCoefficients masksAndCoefficients(Mesh, DerivTwo, NAdvCellsForEdge,
AdvCellsForEdge, AdvMaskHighOrder,
AdvCoefs, AdvCoefs3rd);
MasksAndCoefficients masksAndCoefficients(
Mesh, VCoord, DerivTwo, NAdvCellsForEdge, AdvCellsForEdge,
AdvMaskHighOrder, AdvCoefs, AdvCoefs3rd);
Kokkos::fence();
parallelFor(
{NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); });
Expand Down
50 changes: 27 additions & 23 deletions components/omega/src/ocn/TendencyTerms.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -381,32 +381,35 @@ class TracerHorzAdvOnCell {
const I4 KEnd = KStart + VecLength;
for (int K = KStart; K < KEnd; ++K)
HighOrderFlxHorz(L, IEdge, K) = 0;
if (!ForceLowOrder && AdvMaskHighOrder(IEdge)) {
for (int I = 0; I < NAdvCellsForEdge(IEdge); ++I) {
const I4 ICell = AdvCellsForEdge(IEdge, I);
for (int K = KStart; K < KEnd; ++K) {
const Real NormalThicknessFlux =
FluxLayerThickEdge(IEdge, K) * NormVelEdge(IEdge, K);
const Real TracerWgt =
(AdvCoefs(I, IEdge) +
Coef3rdOrder * std::copysign(1._Real, NormalThicknessFlux) *
AdvCoefs3rd(I, IEdge)) *
NormalThicknessFlux;
HighOrderFlxHorz(L, IEdge, K) +=
TracerWgt * TracerCell(L, ICell, K);
}
}
} else {

// Stay at low order at boundaries
for (int K = KStart; K < KEnd; ++K) {
const I4 JCell0 = CellsOnEdge(IEdge, 0);
const I4 JCell1 = CellsOnEdge(IEdge, 1);
const Real NormalThicknessFlux =
FluxLayerThickEdge(IEdge, K) * NormVelEdge(IEdge, K);
const Real TracerWgt = DvEdge(IEdge) * 0.5_Real * NormalThicknessFlux;
HighOrderFlxHorz(L, IEdge, K) +=
TracerWgt * (1._Real - AdvMaskHighOrder(IEdge, K)) *
(TracerCell(L, JCell1, K) + TracerCell(L, JCell0, K));
}

// High order (3rd or 4th) fluxes elsewhere when requested
// - If HorzTracerFluxOrder = 2, NAdvCellsForEdge = 0 and
// this loop is skipped.
for (int I = 0; I < NAdvCellsForEdge(IEdge); ++I) {
const I4 ICell = AdvCellsForEdge(IEdge, I);
for (int K = KStart; K < KEnd; ++K) {
const I4 JCell0 = CellsOnEdge(IEdge, 0);
const I4 JCell1 = CellsOnEdge(IEdge, 1);
const Real NormalThicknessFlux =
FluxLayerThickEdge(IEdge, K) * NormVelEdge(IEdge, K);
const Real TracerWgt =
DvEdge(IEdge) * 0.5_Real * NormalThicknessFlux;
HighOrderFlxHorz(L, IEdge, K) +=
TracerWgt *
(TracerCell(L, JCell1, K) + TracerCell(L, JCell0, K));
(AdvCoefs(I, IEdge) +
Coef3rdOrder * std::copysign(1._Real, NormalThicknessFlux) *
AdvCoefs3rd(I, IEdge)) *
NormalThicknessFlux;
HighOrderFlxHorz(L, IEdge, K) += TracerWgt *
TracerCell(L, ICell, K) *
AdvMaskHighOrder(IEdge, K);
}
}
}
Expand All @@ -430,9 +433,10 @@ class TracerHorzAdvOnCell {

private:
const HorzMesh *HorzontalMesh;
const VertCoord *VerticalCoord;
Array1DI4 NAdvCellsForEdge;
Array2DI4 AdvCellsForEdge;
Array1DI4 AdvMaskHighOrder;
Array2DI4 AdvMaskHighOrder;
Array2DReal AdvCoefs;
Array2DReal AdvCoefs3rd;
Array3DReal HighOrderFlxHorz;
Expand Down
55 changes: 50 additions & 5 deletions components/omega/src/ocn/VertCoord.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,8 +148,9 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
VertexDegree = Decomp->VertexDegree;

// Retrieve connectivity arrays from Decomp
CellsOnEdge = Decomp->CellsOnEdge;
CellsOnVertex = Decomp->CellsOnVertex;
CellsOnEdge = Decomp->CellsOnEdge;
CellsOnVertex = Decomp->CellsOnVertex;
VerticesOnEdge = Decomp->VerticesOnEdge;

// Allocate device arrays
MaxLayerCell = Array1DI4("MaxLayerCell", NCellsSize);
Expand Down Expand Up @@ -758,14 +759,19 @@ void VertCoord::minMaxLayerVertex(Halo *MeshHalo) {
// set computational masks for mesh elements
void VertCoord::setMasks() {

EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLayers);
// EdgeMask & BoundaryEdge
EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLayers);
BoundaryEdge = Array2DReal("BoundaryEdge", NEdgesSize, NVertLayers);

OMEGA_SCOPE(LocEdgeMask, EdgeMask);
OMEGA_SCOPE(LocBoundaryEdge, BoundaryEdge);
OMEGA_SCOPE(LocMinLyrEdgeBot, MinLayerEdgeBot);
OMEGA_SCOPE(LocMaxLyrEdgeTop, MaxLayerEdgeTop);

// EdgeMask = 1 if active layers on both sides, 0 otherwise.
// BoundaryEdge = 0 if active layers on both sides, 1 otherwise.
deepCopy(EdgeMask, 0.);
deepCopy(BoundaryEdge, 1.);
parallelForOuter(
{NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
const I4 KMin = LocMinLyrEdgeBot(IEdge);
Expand All @@ -775,12 +781,50 @@ void VertCoord::setMasks() {
Team, KMax - KMin + 1, INNER_LAMBDA(int K) {
I4 KLyr = KMin + K;

LocEdgeMask(IEdge, KLyr) = 1._Real;
LocEdgeMask(IEdge, KLyr) = 1._Real;
LocBoundaryEdge(IEdge, KLyr) = 0._Real;
});
});

EdgeMaskH = createHostMirrorCopy(EdgeMask);
EdgeMaskH = createHostMirrorCopy(EdgeMask);
BoundaryEdgeH = createHostMirrorCopy(BoundaryEdge);

// BoundaryCell & BoundaryVertex
// = 1 in inactive layers, 0 otherwise.
BoundaryCell = Array2DReal("BoundaryCell", NCellsSize, NVertLayers);
BoundaryVertex = Array2DReal("BoundaryVertex", NVerticesSize, NVertLayers);

OMEGA_SCOPE(LocBoundaryCell, BoundaryCell);
OMEGA_SCOPE(LocBoundaryVertex, BoundaryVertex);
OMEGA_SCOPE(LocCellsOnEdge, CellsOnEdge);
OMEGA_SCOPE(LocVerticesOnEdge, CellsOnEdge);
OMEGA_SCOPE(LocNVertLayers, NVertLayers);

deepCopy(BoundaryCell, 0.);
deepCopy(BoundaryVertex, 0.);

parallelForOuter(
{NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
const I4 JCell0 = LocCellsOnEdge(IEdge, 0);
const I4 JCell1 = LocCellsOnEdge(IEdge, 1);
const I4 JVertex0 = LocVerticesOnEdge(IEdge, 0);
const I4 JVertex1 = LocVerticesOnEdge(IEdge, 1);

parallelForInner(
Team, LocNVertLayers, INNER_LAMBDA(int K) {
if (LocBoundaryEdge(IEdge, K) == 1._Real) {
LocBoundaryCell(JCell0, K) = 1._Real;
LocBoundaryCell(JCell1, K) = 1._Real;
LocBoundaryVertex(JVertex0, K) = 1._Real;
LocBoundaryVertex(JVertex1, K) = 1._Real;
}
});
});

BoundaryCellH = createHostMirrorCopy(BoundaryCell);
BoundaryVertexH = createHostMirrorCopy(BoundaryVertex);

// CellMask
CellMask = Array2DReal("CellMask", NCellsSize, NVertLayers);

OMEGA_SCOPE(LocCellMask, CellMask);
Expand All @@ -804,6 +848,7 @@ void VertCoord::setMasks() {

CellMaskH = createHostMirrorCopy(CellMask);

// VertexMask
VertexMask = Array2DReal("VertexMask", NVerticesSize, NVertLayers);

OMEGA_SCOPE(LocVrtxMask, VertexMask);
Expand Down
13 changes: 13 additions & 0 deletions components/omega/src/ocn/VertCoord.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,7 @@ class VertCoord {
I4 VertexDegree;
Array2DI4 CellsOnEdge;
Array2DI4 CellsOnVertex;
Array2DI4 VerticesOnEdge;

// Choice of VertCoorMovementWeight type
MovementWeightType MvmtWgtChoice;
Expand Down Expand Up @@ -161,6 +162,18 @@ class VertCoord {
HostArray2DReal VertexMaskH; ///< Mask to determine if computations should be
/// done on vertex

Array2DReal BoundaryEdge; ///< Mask to determine boundary edges
Copy link
Copy Markdown

@overfelt overfelt May 12, 2026

Choose a reason for hiding this comment

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

If these mask arrays are only 0s and 1s, would it be better to make them I4 instead of Real8.

Copy link
Copy Markdown
Author

@hyungyukang hyungyukang May 13, 2026

Choose a reason for hiding this comment

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

@overfelt , thanks a lot for reviewing!
It is true that the Boundary* arrays contain only 0 or 1 values. But I’m not sure whether I4 is the most appropriate type for these boundary mask arrays. In MPAS-Ocean, all mask arrays including AdvMaskHighOrder are real-valued. If we change the Boundary* arrays to I4, we may also need to consider changing other mask arrays, such as CellMask, EdgeMask, and VertexMask, to I4 for consistency, which would require a fair amount of work.
For now, I would suggest leaving this as is. We can revisit it later and update all related mask arrays together if needed.

Copy link
Copy Markdown

@xylar xylar May 13, 2026

Choose a reason for hiding this comment

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

@overfelt and @hyungyukang, please make an issue about this so it doesn't get lost. I do agree with @overfelt that these don't need to be Real, they should be I4.


HostArray2DReal BoundaryEdgeH; ///< Mask to determine boundary edges

Array2DReal BoundaryCell; ///< Mask to determine boundary cells

HostArray2DReal BoundaryCellH; ///< Mask to determine boundary cells

Array2DReal BoundaryVertex; ///< Mask to determine boundary vertex

HostArray2DReal BoundaryVertexH; ///< Mask to determine boundary vertex

// BottomDepth read from mesh file
Array1DReal BottomDepth;

Expand Down
2 changes: 1 addition & 1 deletion components/omega/test/ocn/TendencyTermsTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -879,8 +879,8 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) {
Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned,
NVertLayers);
TracerHorzAdvOnCell TrHorzAdvOnC(Mesh, VCoord);
TrHorzAdvOnC.init();
TrHorzAdvOnC.ForceLowOrder = true;
TrHorzAdvOnC.init();

parallelFor(
{NTracers, Mesh->NEdgesAll, NVertLayers},
Expand Down
Loading