Skip to content

Commit f6cb81b

Browse files
vanroekelvanroekel
committed
Fixes ssh calculation
Moves SSH calculation to the VertCoord class and computes it as a single layer instead of a multi level field.
1 parent 33decbc commit f6cb81b

File tree

8 files changed

+39
-64
lines changed

8 files changed

+39
-64
lines changed

‎components/omega/src/ocn/AuxiliaryState.cpp‎

Lines changed: 0 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -186,22 +186,6 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel,
186186
});
187187
Pacer::stop("AuxState:cellAuxState2", 2);
188188

189-
Pacer::start("AuxState:cellAuxState3", 2);
190-
parallelForOuter(
191-
"cellAuxState3", {Mesh->NCellsAll},
192-
KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
193-
const int KMin = MinLayerCell(ICell);
194-
const int KMax = MaxLayerCell(ICell);
195-
const int KRange = vertRangeChunked(KMin, KMax);
196-
197-
parallelForInner(
198-
Team, KRange, INNER_LAMBDA(int KChunk) {
199-
LocLayerThicknessAux.computeVarsOnCells(ICell, KChunk,
200-
LayerThickCell);
201-
});
202-
});
203-
Pacer::stop("AuxState:cellAuxState3", 2);
204-
205189
Pacer::stop("AuxState:computeMomAux", 1);
206190
}
207191

‎components/omega/src/ocn/Tendencies.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -332,6 +332,7 @@ void Tendencies::computeVelocityTendenciesOnly(
332332
OMEGA_SCOPE(LocBottomDrag, BottomDrag);
333333
OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot);
334334
OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop);
335+
OMEGA_SCOPE(LocSshCell, VCoord->SshCell);
335336

336337
Pacer::start("Tend:computeVelocityTendenciesOnly", 1);
337338

@@ -393,7 +394,7 @@ void Tendencies::computeVelocityTendenciesOnly(
393394
}
394395

395396
// Compute sea surface height gradient
396-
const Array2DReal &SSHCell = AuxState->LayerThicknessAux.SshCell;
397+
const Array1DReal &SSHCell = LocSshCell;
397398
if (LocSSHGrad.Enabled) {
398399
Pacer::start("Tend:SSHGrad", 2);
399400
parallelForOuter(

‎components/omega/src/ocn/TendencyTerms.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,7 @@ class SSHGradOnEdge {
174174
/// The functor takes edge index, vertical chunk index, and array of
175175
/// layer thickness/SSH, outputs tendency array
176176
KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk,
177-
const Array2DReal &SshCell) const {
177+
const Array1DReal &SshCell) const {
178178

179179
const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge));
180180
const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge));
@@ -185,7 +185,7 @@ class SSHGradOnEdge {
185185
for (int KVec = 0; KVec < KLen; ++KVec) {
186186
const I4 K = KStart + KVec;
187187
Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity *
188-
(SshCell(ICell1, K) - SshCell(ICell0, K)) *
188+
(SshCell(ICell1) - SshCell(ICell0)) *
189189
InvDcEdge;
190190
}
191191
}

‎components/omega/src/ocn/VertCoord.cpp‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -163,6 +163,8 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
163163
PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers);
164164
ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1);
165165
ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers);
166+
SshCell = Array1DReal("SshCell", NCellsSize);
167+
166168
GeopotentialMid = Array2DReal("GeopotentialMid", NCellsSize, NVertLayers);
167169
LayerThicknessTarget =
168170
Array2DReal("LayerThicknessTarget", NCellsSize, NVertLayers);
@@ -174,6 +176,8 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
174176
PressureMidH = createHostMirrorCopy(PressureMid);
175177
ZInterfaceH = createHostMirrorCopy(ZInterface);
176178
ZMidH = createHostMirrorCopy(ZMid);
179+
SshCellH = createHostMirrorCopy(SshCellH);
180+
177181
GeopotentialMidH = createHostMirrorCopy(GeopotentialMid);
178182
LayerThicknessTargetH = createHostMirrorCopy(LayerThicknessTarget);
179183
RefLayerThicknessH = createHostMirrorCopy(RefLayerThickness);
@@ -236,6 +240,8 @@ void VertCoord::defineFields() {
236240
PressMidFldName = "PressureMid";
237241
ZInterfFldName = "ZInterface";
238242
ZMidFldName = "ZMid";
243+
SshFldName = "SshCell";
244+
239245
GeopotFldName = "GeopotentialMid";
240246
LyrThickTargetFldName = "LayerThicknessTarget";
241247

@@ -249,6 +255,7 @@ void VertCoord::defineFields() {
249255
ZMidFldName.append(Name);
250256
GeopotFldName.append(Name);
251257
LyrThickTargetFldName.append(Name);
258+
SshFldName.append(Name);
252259
}
253260

254261
// Create fields for VertCoord variables
@@ -295,6 +302,19 @@ void VertCoord::defineFields() {
295302
DimNames // dimension names
296303
);
297304

305+
auto SshField = Field::create(
306+
SshFldName, // field name
307+
"sea surface height at cell center", // long Name or description
308+
"m", // units
309+
"sea_surface_height", // CF standard Name
310+
std::numeric_limits<Real>::min(), // min valid value
311+
std::numeric_limits<Real>::max(), // max valid value
312+
FillValueReal, // scalar for undefined entries
313+
NDims, // number of dimensions
314+
DimNames // dimension names
315+
);
316+
317+
298318
NDims = 2;
299319
DimNames.resize(NDims);
300320
DimNames[1] = "NVertLayersP1";
@@ -402,6 +422,7 @@ void VertCoord::defineFields() {
402422
VCoordGroup->addField(ZMidFldName);
403423
VCoordGroup->addField(GeopotFldName);
404424
VCoordGroup->addField(LyrThickTargetFldName);
425+
VCoordGroup->addField(SshFldName);
405426

406427
// Associate Field with data
407428
PressureInterfaceField->attachData<Array2DReal>(PressureInterface);
@@ -410,6 +431,7 @@ void VertCoord::defineFields() {
410431
ZMidField->attachData<Array2DReal>(ZMid);
411432
GeopotentialMidField->attachData<Array2DReal>(GeopotentialMid);
412433
LayerThicknessTargetField->attachData<Array2DReal>(LayerThicknessTarget);
434+
SshField->attachData<Array1DReal>(SshCell);
413435

414436
} // end defineFields
415437

@@ -431,6 +453,7 @@ VertCoord::~VertCoord() {
431453
Field::destroy(ZMidFldName);
432454
Field::destroy(GeopotFldName);
433455
Field::destroy(LyrThickTargetFldName);
456+
Field::destroy(SshFldName);
434457
FieldGroup::destroy(GroupName);
435458
}
436459

@@ -872,6 +895,7 @@ void VertCoord::computeZHeight(
872895
OMEGA_SCOPE(LocZInterf, ZInterface);
873896
OMEGA_SCOPE(LocZMid, ZMid);
874897
OMEGA_SCOPE(LocBotDepth, BottomDepth);
898+
OMEGA_SCOPE(LocSshCell, SshCell);
875899

876900
const auto Policy = TeamPolicy(NCellsAll, OMEGA_TEAMSIZE, 1);
877901
Kokkos::parallel_for(
@@ -893,6 +917,9 @@ void VertCoord::computeZHeight(
893917
LocZInterf(ICell, KLyr) = -LocBotDepth(ICell) + Accum;
894918
LocZMid(ICell, KLyr) =
895919
-LocBotDepth(ICell) + Accum - 0.5 * DZ;
920+
if (KLyr == 0) {
921+
LocSshCell(ICell) = LocZInterf(ICell,KLyr);
922+
}
896923
}
897924
});
898925
});
@@ -1005,6 +1032,8 @@ void VertCoord::copyToHost() {
10051032
deepCopy(PressureMidH, PressureMid);
10061033
deepCopy(ZInterfaceH, ZInterface);
10071034
deepCopy(ZMidH, ZMid);
1035+
deepCopy(SshCellH, SshCell);
1036+
10081037
deepCopy(GeopotentialMidH, GeopotentialMid);
10091038
deepCopy(LayerThicknessTargetH, LayerThicknessTarget);
10101039
deepCopy(RefLayerThicknessH, RefLayerThickness);
@@ -1018,6 +1047,8 @@ void VertCoord::copyToDevice() {
10181047
deepCopy(PressureMid, PressureMidH);
10191048
deepCopy(ZInterface, ZInterfaceH);
10201049
deepCopy(ZMid, ZMidH);
1050+
deepCopy(SshCell, SshCellH);
1051+
10211052
deepCopy(GeopotentialMid, GeopotentialMidH);
10221053
deepCopy(LayerThicknessTarget, LayerThicknessTargetH);
10231054
deepCopy(RefLayerThickness, RefLayerThicknessH);

‎components/omega/src/ocn/VertCoord.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -107,13 +107,15 @@ class VertCoord {
107107
Array2DReal ZMid;
108108
Array2DReal GeopotentialMid;
109109
Array2DReal LayerThicknessTarget;
110+
Array1DReal SshCell;
110111

111112
HostArray2DReal PressureInterfaceH;
112113
HostArray2DReal PressureMidH;
113114
HostArray2DReal ZInterfaceH;
114115
HostArray2DReal ZMidH;
115116
HostArray2DReal GeopotentialMidH;
116117
HostArray2DReal LayerThicknessTargetH;
118+
HostArray1DReal SshCellH;
117119

118120
// Vertical loop bounds
119121
Array1DI4 MinLayerCell;
@@ -181,6 +183,7 @@ class VertCoord {
181183
std::string ZMidFldName; ///< Field name for midpoint Z height
182184
std::string GeopotFldName; ///< Field name for geopotential
183185
std::string LyrThickTargetFldName; ///< Field name for target thickness
186+
std::string SshFldName; ///< Field name for sea surface height
184187

185188
// methods
186189

‎components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp‎

Lines changed: 0 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -12,8 +12,6 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix,
1212
Mesh->NEdgesSize, VCoord->NVertLayers),
1313
MeanLayerThickEdge("MeanLayerThickEdge" + AuxStateSuffix,
1414
Mesh->NEdgesSize, VCoord->NVertLayers),
15-
SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize,
16-
VCoord->NVertLayers),
1715
CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth),
1816
MinLayerEdgeBot(VCoord->MinLayerEdgeBot),
1917
MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop),
@@ -64,35 +62,18 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName,
6462
DimNames // dimension names
6563
);
6664

67-
// Sea surface height
68-
DimNames[0] = "NCells" + DimSuffix;
69-
auto SshCellField = Field::create(
70-
SshCell.label(), // field name
71-
"sea surface height at cell center", // long Name or description
72-
"m", // units
73-
"sea_surface_height", // CF standard Name
74-
0, // min valid value
75-
std::numeric_limits<Real>::max(), // max valid value
76-
FillValue, // scalar for undefined entries
77-
NDims, // number of dimensions
78-
DimNames // dimension names
79-
);
80-
8165
// Add fields to Aux field group
8266
FieldGroup::addFieldToGroup(FluxLayerThickEdge.label(), AuxGroupName);
8367
FieldGroup::addFieldToGroup(MeanLayerThickEdge.label(), AuxGroupName);
84-
FieldGroup::addFieldToGroup(SshCell.label(), AuxGroupName);
8568

8669
// Attach field data
8770
FluxLayerThickEdgeField->attachData<Array2DReal>(FluxLayerThickEdge);
8871
MeanLayerThickEdgeField->attachData<Array2DReal>(MeanLayerThickEdge);
89-
SshCellField->attachData<Array2DReal>(SshCell);
9072
}
9173

9274
void LayerThicknessAuxVars::unregisterFields() const {
9375
Field::destroy(FluxLayerThickEdge.label());
9476
Field::destroy(MeanLayerThickEdge.label());
95-
Field::destroy(SshCell.label());
9677
}
9778

9879
} // namespace OMEGA

‎components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h‎

Lines changed: 0 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,6 @@ class LayerThicknessAuxVars {
1616
public:
1717
Array2DReal FluxLayerThickEdge;
1818
Array2DReal MeanLayerThickEdge;
19-
Array2DReal SshCell;
2019

2120
FluxThickEdgeOption FluxThickEdgeChoice;
2221

@@ -63,29 +62,6 @@ class LayerThicknessAuxVars {
6362
}
6463
}
6564

66-
KOKKOS_FUNCTION void
67-
computeVarsOnCells(int ICell, int KChunk,
68-
const Array2DReal &LayerThickCell) const {
69-
70-
// Temporary for stacked shallow water
71-
const int KStart = chunkStart(KChunk, MinLayerCell(ICell));
72-
const int KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell));
73-
74-
for (int KVec = 0; KVec < KLen; ++KVec) {
75-
const int K = KStart + KVec;
76-
SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell);
77-
}
78-
79-
/*
80-
Real TotalThickness = 0.0;
81-
for (int K = 0; K < NVertLayers; K++) {
82-
TotalThickness += LayerThickCell(ICell, K);
83-
}
84-
85-
SshCell(ICell) = TotalThickness - BottomDepth(ICell);
86-
*/
87-
}
88-
8965
void registerFields(const std::string &AuxGroupName,
9066
const std::string &MeshName) const;
9167
void unregisterFields() const;

‎components/omega/test/ocn/TendencyTermsTest.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -529,8 +529,7 @@ int testSSHGrad(int NVertLayers, Real RTol) {
529529
},
530530
ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No);
531531

532-
// Set input array
533-
Array2DReal SSHCell("SSHCell", Mesh->NCellsSize, NVertLayers);
532+
Array1DReal SSHCell("SSHCell", Mesh->NCellsSize);
534533

535534
Err += setScalar(
536535
KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, SSHCell,

0 commit comments

Comments
 (0)