Skip to content

Change Brunt-Vaisala frequency vertical dimension and fill boundary values #404

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
mwarusz wants to merge 1 commit into
base: develop
Choose a base branch
from
+73 −59
Open

Change Brunt-Vaisala frequency vertical dimension and fill boundary values #404

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
40 changes: 36 additions & 4 deletions components/omega/src/ocn/Eos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,8 +46,8 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object
SpecVol = Array2DReal("SpecVol", Mesh->NCellsSize, VCoord->NVertLayers);
SpecVolDisplaced =
Array2DReal("SpecVolDisplaced", Mesh->NCellsSize, VCoord->NVertLayers);
BruntVaisalaFreqSq =
Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize, VCoord->NVertLayers);
BruntVaisalaFreqSq = Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize,
VCoord->NVertLayersP1);

deepCopy(SpecVol, 1.0_Real / RhoSw);
deepCopy(SpecVolDisplaced, 1.0_Real / RhoSw);
Expand Down Expand Up @@ -282,22 +282,37 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp,
parallelForOuter(
"bvf-linear", {Mesh->NCellsAll},
KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) {
const int KMin = MinLayerCell(ICell);
// Compute Brunt-Vaisala frequency at interior vertical interfaces
const int KMin = MinLayerCell(ICell) + 1;
const int KMax = MaxLayerCell(ICell);
const int KRange = vertRangeChunked(KMin, KMax);
parallelForInner(
Team, KRange, INNER_LAMBDA(int KChunk) {
LocComputeBruntVaisalaFreqSqLinear(LocBruntVaisalaFreqSq,
ICell, KChunk, SpecVol);
});

teamBarrier(Team);

// Fill Brunt-Vaisala frequency at vertical boundaries using the
// closest valid value. This is equivalent to doing one-sided
// differencing at the boundary.
Kokkos::single(
PerTeam(Team), INNER_LAMBDA() {
LocBruntVaisalaFreqSq(ICell, MinLayerCell(ICell)) =
LocBruntVaisalaFreqSq(ICell, KMin);
LocBruntVaisalaFreqSq(ICell, MaxLayerCell(ICell) + 1) =
LocBruntVaisalaFreqSq(ICell, KMax);
});
});
} else if (EosChoice == EosType::Teos10Eos) {
/// If TEOS-10 EOS, use TEOS-10 squared Brunt-Vaisala frequency
/// calculation
parallelForOuter(
"bvf-teos10", {Mesh->NCellsAll},
KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) {
const int KMin = MinLayerCell(ICell);
// Compute Brunt-Vaisala frequency at interior vertical interfaces
const int KMin = MinLayerCell(ICell) + 1;
const int KMax = MaxLayerCell(ICell);
const int KRange = vertRangeChunked(KMin, KMax);
parallelForInner(
Expand All @@ -306,6 +321,19 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp,
LocBruntVaisalaFreqSq, ICell, KChunk, ConservTemp,
AbsSalinity, Pressure, SpecVol);
});

teamBarrier(Team);

// Fill Brunt-Vaisala frequency at vertical boundaries using the
// closest valid value. This is equivalent to doing one-sided
// differencing at the boundary.
Kokkos::single(
PerTeam(Team), INNER_LAMBDA() {
LocBruntVaisalaFreqSq(ICell, MinLayerCell(ICell)) =
LocBruntVaisalaFreqSq(ICell, KMin);
LocBruntVaisalaFreqSq(ICell, MaxLayerCell(ICell) + 1) =
LocBruntVaisalaFreqSq(ICell, KMax);
});
});
}
}
Expand Down Expand Up @@ -355,6 +383,10 @@ void Eos::defineFields() {
NDims, // Number of dimensions
DimNames // Dimension names
);

// Brunt-Vaisala frequency is located at interfaces
DimNames[1] = "NVertLayersP1";

/// Create and register the BruntVaisalaFreqSq field
auto BruntVaisalaFreqSqField =
Field::create(BruntVaisalaFreqSqFldName, // Field name
Expand Down
62 changes: 25 additions & 37 deletions components/omega/src/ocn/Eos.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -324,33 +324,27 @@ class Teos10BruntVaisalaFreqSq {
const Array2DReal &Pressure,
const Array2DReal &SpecVol) const {

const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell));
const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1);
const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell));

for (int KVec = 0; KVec < KLen; ++KVec) {
const I4 K = KStart + KVec;
if (K == 0) {
// No Brunt-Vaisala frequency at surface
BruntVaisalaFreqSq(ICell, K) = 0.0_Real;
} else {
// Calculate squared Brunt-Vaisala frequency
Real CtInt =
0.5_Real * (ConservTemp(ICell, K) + ConservTemp(ICell, K - 1));
Real SaInt =
0.5_Real * (AbsSalinity(ICell, K) + AbsSalinity(ICell, K - 1));
Real PInt =
0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1));
Real SpInt = 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1));
Real AlphaInt = calcAlpha(SaInt, CtInt, PInt * Pa2Db, SpInt);
Real BetaInt = calcBeta(SaInt, CtInt, PInt * Pa2Db, SpInt);
Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1);
Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1);
Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1);

BruntVaisalaFreqSq(ICell, K) = Gravity * Gravity *
(BetaInt * DSa - AlphaInt * DCt) /
(SpInt * DP);
}
// Calculate squared Brunt-Vaisala frequency
Real CtInt =
0.5_Real * (ConservTemp(ICell, K) + ConservTemp(ICell, K - 1));
Real SaInt =
0.5_Real * (AbsSalinity(ICell, K) + AbsSalinity(ICell, K - 1));
Real PInt = 0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1));
Real SpInt = 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1));
Real AlphaInt = calcAlpha(SaInt, CtInt, PInt * Pa2Db, SpInt);
Real BetaInt = calcBeta(SaInt, CtInt, PInt * Pa2Db, SpInt);
Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1);
Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1);
Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1);

BruntVaisalaFreqSq(ICell, K) = Gravity * Gravity *
(BetaInt * DSa - AlphaInt * DCt) /
(SpInt * DP);
}
}

Expand Down Expand Up @@ -533,24 +527,18 @@ class LinearBruntVaisalaFreqSq {
I4 KChunk,
const Array2DReal &SpecVol) const {

const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell));
const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1);
const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell));

for (int KVec = 0; KVec < KLen; ++KVec) {
const I4 K = KStart + KVec;
if (K == 0) {
/// No Brunt-Vaisala frequency at the top level
BruntVaisalaFreqSq(ICell, K) = 0.0_Real;
} else {
/// Calculate squared Brunt-Vaisala frequency at mid-point between
/// K-1 and K Do not need to use displaced specific volume here
/// since only the linear EOS is used with this BVF formulation.
BruntVaisalaFreqSq(ICell, K) =
-(Gravity / RhoT0S0) *
((1.0_Real / SpecVol(ICell, K - 1)) -
(1.0_Real / SpecVol(ICell, K))) /
(ZMid(ICell, K - 1) - ZMid(ICell, K));
}
/// Calculate squared Brunt-Vaisala frequency at mid-point between
/// K-1 and K Do not need to use displaced specific volume here
/// since only the linear EOS is used with this BVF formulation.
BruntVaisalaFreqSq(ICell, K) = -(Gravity / RhoT0S0) *
((1.0_Real / SpecVol(ICell, K - 1)) -
(1.0_Real / SpecVol(ICell, K))) /
(ZMid(ICell, K - 1) - ZMid(ICell, K));
}
}

Expand Down
30 changes: 12 additions & 18 deletions components/omega/test/ocn/EosTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -379,16 +379,13 @@ void testBruntVaisalaFreqSqLinear() {
KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) {
int NumMismatchesCol;
const int KMin = MinLayerCell(ICell);
const int KMax = MaxLayerCell(ICell);
const int KMax = MaxLayerCell(ICell) + 1;
const int KRange = vertRange(KMin, KMax);
parallelReduceInner(
Team, KRange,
INNER_LAMBDA(int KOff, int &InnerCount) {
const int K = KMin + KOff;
if (K == 0) { // should be zero
if (BruntVaisalaFreqSq(ICell, K) != 0.0)
InnerCount++;
} else if (K == 1) { // should be ref value
if (K == 1 || K == 0) { // should be ref value
if (!isApprox(BruntVaisalaFreqSq(ICell, K),
LinearBVFExpValue, RTol))
InnerCount++;
Expand All @@ -410,18 +407,18 @@ void testBruntVaisalaFreqSqLinear() {
if (NumMismatches != 0) {
auto BruntVaisalaFreqSqH = createHostMirrorCopy(BruntVaisalaFreqSq);
for (int I = 0; I < Mesh->NCellsAll; ++I) {
// top layer should be zero
if (BruntVaisalaFreqSqH(I, 0) != 0.0)
// top layer should be ref value
if (!isApprox(BruntVaisalaFreqSqH(I, 0), LinearBVFExpValue, RTol))
LOG_ERROR("EosTest: Brunt-Vaisala Linear Bad Value: "
"BruntVaisala({},{}) = {}; Expected {}",
I, 0, BruntVaisalaFreqSqH(I, 0), 0.0);
I, 0, BruntVaisalaFreqSqH(I, 0), LinearBVFExpValue);
// K = 1 should be ref value
if (!isApprox(BruntVaisalaFreqSqH(I, 1), LinearBVFExpValue, RTol))
LOG_ERROR("EosTest: Brunt-Vaisala Linear Bad Value: "
"BruntVaisala({},{}) = {}; Expected {}",
I, 1, BruntVaisalaFreqSqH(I, 1), LinearBVFExpValue);
// remaining values just check for other conditions
for (int K = 2; K < NVertLayers; ++K) {
for (int K = 2; K < NVertLayers + 1; ++K) {
if (BruntVaisalaFreqSqH(I, K) == 0.0 or
Kokkos::isnan(BruntVaisalaFreqSqH(I, K)) or
Kokkos::isinf(BruntVaisalaFreqSqH(I, K)))
Expand Down Expand Up @@ -643,16 +640,13 @@ void testBruntVaisalaFreqSqTeos10() {
KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) {
int NumMismatchesCol;
const int KMin = MinLayerCell(ICell);
const int KMax = MaxLayerCell(ICell);
const int KMax = MaxLayerCell(ICell) + 1;
const int KRange = vertRange(KMin, KMax);
parallelReduceInner(
Team, KRange,
INNER_LAMBDA(int KOff, int &InnerCount) {
const int K = KMin + KOff;
if (K == 0) { // should be zero at top
if (BruntVaisalaFreqSq(ICell, K) != 0.0)
InnerCount++;
} else if (K == 1) { // should be ref value
if (K == 1 || K == 0) { // should be ref value
if (!isApprox(BruntVaisalaFreqSq(ICell, K), TeosBVFExpValue,
RTol))
InnerCount++;
Expand All @@ -674,18 +668,18 @@ void testBruntVaisalaFreqSqTeos10() {
if (NumMismatches != 0) {
auto BruntVaisalaFreqSqH = createHostMirrorCopy(BruntVaisalaFreqSq);
for (int ICell = 0; ICell < Mesh->NCellsAll; ++ICell) {
// top layer should be zero
if (BruntVaisalaFreqSqH(ICell, 0) != 0.0)
// top layer should be ref value
if (!isApprox(BruntVaisalaFreqSqH(ICell, 0), TeosBVFExpValue, RTol))
LOG_ERROR("EosTest: Brunt-Vaisala TEOS Bad Value: "
"BruntVaisala({},{}) = {}; Expected {}",
ICell, 0, BruntVaisalaFreqSqH(ICell, 0), 0.0);
ICell, 0, BruntVaisalaFreqSqH(ICell, 0), TeosBVFExpValue);
// K = 1 should be ref value
if (!isApprox(BruntVaisalaFreqSqH(ICell, 1), TeosBVFExpValue, RTol))
LOG_ERROR("EosTest: Brunt-Vaisala TEOS Bad Value: "
"BruntVaisala({},{}) = {}; Expected {}",
ICell, 1, BruntVaisalaFreqSqH(ICell, 1), TeosBVFExpValue);
// remaining values just check for other conditions
for (int K = 2; K < NVertLayers; ++K) {
for (int K = 2; K < NVertLayers + 1; ++K) {
if (BruntVaisalaFreqSqH(ICell, K) == 0.0 or
Kokkos::isnan(BruntVaisalaFreqSqH(ICell, K)) or
Kokkos::isinf(BruntVaisalaFreqSqH(ICell, K)))
Expand Down
Loading