Skip to content

Commit 3f00564

Browse files
alicebarthelalicebarthel
committed
switched clamping to an option, and eos limits has 2 options
1 parent 30c5666 commit 3f00564

File tree

4 files changed

+90
-42
lines changed

4 files changed

+90
-42
lines changed

‎components/omega/configs/Default.yml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,8 @@ Omega:
5959
DRhoDT: -0.2
6060
DRhoDS: 0.8
6161
RhoT0S0: 1000.0
62+
ClampingEnable: false
63+
EosLimits: Funnel
6264
IOStreams:
6365
InitialVertCoord:
6466
UsePointerFile: false

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

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,12 @@ void Eos::init() {
8181
Err += EosConfig.get("EosType", EosTypeStr);
8282
CHECK_ERROR_ABORT(Err, "Eos::init: EosType subgroup not found in EosConfig");
8383

84+
std::string EosLimStr;
85+
Err += EosConfig.get("EosLimits", EosLimStr);
86+
CHECK_ERROR_ABORT(Err,
87+
"Eos::init: EosLimits subgroup not found in EosConfig");
88+
/// EosLimits EosLimChoice;
89+
8490
/// Set EosChoice and parameters based on EosTypeStr
8591
if (EosTypeStr == "Linear" or EosTypeStr == "linear") {
8692
Config EosLinConfig("Linear");
@@ -104,12 +110,23 @@ void Eos::init() {
104110
} else if ((EosTypeStr == "teos10") or (EosTypeStr == "teos-10") or
105111
(EosTypeStr == "TEOS-10")) {
106112
eos->EosChoice = EosType::Teos10Eos;
113+
if (EosLimStr == "Funnel") {
114+
eos->ComputeSpecVolTeos10.EosLimChoice = EosLimits::Funnel;
115+
} else if (EosLimStr == "Cube") {
116+
eos->ComputeSpecVolTeos10.EosLimChoice = EosLimits::Cube;
117+
} else {
118+
LOG_ERROR("Eos::init: Unknown EosLimits requested");
119+
Err += EosConfig.get("ClampingEnable",
120+
eos->ComputeSpecVolTeos10.ClampingEnable);
121+
CHECK_ERROR_ABORT(
122+
Err, "Eos::init: Parameter ClampingEnable not found in EosConfig");
123+
}
107124
} else {
108125
LOG_ERROR("Eos::init: Unknown EosType requested");
109126
}
110127
} // end init
111128

112-
/// Compute specific volume for all cells/layers (no displacement)
129+
/// Compute specific volume for all cells/levels (no displacement)
113130
void Eos::computeSpecVol(const Array2DReal &ConservTemp,
114131
const Array2DReal &AbsSalinity,
115132
const Array2DReal &Pressure) {

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

Lines changed: 63 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,11 @@ enum class EosType {
2626
Teos10Eos /// Roquet et al. 2015 75 term expansion
2727
};
2828

29+
enum class EosLimits {
30+
Funnel, ///
31+
Cube ///
32+
};
33+
2934
struct EosRange {
3035
Real Lo;
3136
Real Hi;
@@ -35,6 +40,8 @@ struct EosRange {
3540
class Teos10Eos {
3641
public:
3742
Array2DReal SpecVolPCoeffs;
43+
EosLimits EosLimChoice; ///< EOS clamping option in use
44+
bool ClampingEnable; ///< EOS clamping option in use
3845

3946
/// constructor declaration
4047
Teos10Eos(int NVertLayers);
@@ -90,10 +97,14 @@ class Teos10Eos {
9097
constexpr Real DeltaS = 24.0;
9198
EosRange SRange = calcSLimits(P);
9299
EosRange TRange = calcTLimits(Sa, P);
93-
Real SaInFunnel = Kokkos::clamp(
94-
Sa, SRange.Lo, SRange.Hi); // Salt limited to Poly75t valid range
95-
Real Ss = Kokkos::sqrt((SaInFunnel + DeltaS) / SAu);
96-
Real Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu;
100+
Real Ss = Kokkos::sqrt((Sa + DeltaS) / SAu);
101+
Real Tt = Ct / CTu;
102+
if (ClampingEnable) {
103+
Real SaInRange = Kokkos::clamp(
104+
Sa, SRange.Lo, SRange.Hi); // Salt limited to Poly75t valid range
105+
Ss = Kokkos::sqrt((SaInRange + DeltaS) / SAu);
106+
Tt = Kokkos::clamp(Ct, TRange.Lo, TRange.Hi) / CTu;
107+
}
97108

98109
/// Coefficients for the polynomial expansion
99110
constexpr Real V000 = 1.0769995862e-03;
@@ -215,10 +226,15 @@ class Teos10Eos {
215226
KOKKOS_FUNCTION Real calcDelta(const Array2DReal &SpecVolPCoeffs, const I4 K,
216227
const Real P) const {
217228

218-
constexpr Real Pu = 1e4;
219-
constexpr Real Pmax = 8000.0;
220-
Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
221-
229+
constexpr Real Pu = 1e4;
230+
const Real Pmax = (EosLimChoice == EosLimits::Funnel) ? 8000.0 : 10000.0;
231+
Real Pp = P / Pu;
232+
if (ClampingEnable) {
233+
Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
234+
LOG_INFO("P={} exceeds Pmax={}; Clamping the pressure", P, Pmax);
235+
} else if (P > Pmax) {
236+
LOG_WARN("P={} exceeds Pmax={}", P, Pmax);
237+
}
222238
Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp +
223239
SpecVolPCoeffs(3, K)) *
224240
Pp +
@@ -232,15 +248,20 @@ class Teos10Eos {
232248

233249
/// Calculate reference profile for TEOS-10
234250
KOKKOS_FUNCTION Real calcRefProfile(const Real P) const {
235-
constexpr Real Pu = 1e4;
236-
constexpr Real V00 = -4.4015007269e-05;
237-
constexpr Real V01 = 6.9232335784e-06;
238-
constexpr Real V02 = -7.5004675975e-07;
239-
constexpr Real V03 = 1.7009109288e-08;
240-
constexpr Real V04 = -1.6884162004e-08;
241-
constexpr Real V05 = 1.9613503930e-09;
242-
constexpr Real Pmax = 8000.0;
243-
Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
251+
constexpr Real Pu = 1e4;
252+
constexpr Real V00 = -4.4015007269e-05;
253+
constexpr Real V01 = 6.9232335784e-06;
254+
constexpr Real V02 = -7.5004675975e-07;
255+
constexpr Real V03 = 1.7009109288e-08;
256+
constexpr Real V04 = -1.6884162004e-08;
257+
constexpr Real V05 = 1.9613503930e-09;
258+
const Real Pmax = (EosLimChoice == EosLimits::Funnel) ? 8000.0 : 10000.0;
259+
Real Pp = P / Pu;
260+
if (ClampingEnable) {
261+
Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
262+
} else if (P > Pmax) {
263+
LOG_WARN("P={} exceeds Pmax={}", P, Pmax);
264+
}
244265

245266
Real V0 =
246267
(((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) *
@@ -250,15 +271,17 @@ class Teos10Eos {
250271

251272
/// Calculate S limits of validity given pressure p
252273
KOKKOS_FUNCTION EosRange calcSLimits(const Real P) const {
253-
Real Lo = 0.0;
254-
Real Hi = 42.0;
255-
Real Lo2 = P * 5e-3 - 2.5;
256-
Real Lo3 = 30.0;
257-
258-
if (P >= 500.0 && P < 6500.0) {
259-
Lo = Kokkos::max(Lo, Lo2);
260-
} else if (P >= 6500.0) {
261-
Lo = Kokkos::max(Lo, Lo3);
274+
Real Lo = 0.0;
275+
Real Hi = 42.0;
276+
if (EosLimChoice == EosLimits::Funnel) {
277+
Real Lo2 = P * 5e-3 - 2.5;
278+
Real Lo3 = 30.0;
279+
280+
if (P >= 500.0 && P < 6500.0) {
281+
Lo = Kokkos::max(Lo, Lo2);
282+
} else if (P >= 6500.0) {
283+
Lo = Kokkos::max(Lo, Lo3);
284+
}
262285
}
263286
return {Lo, Hi};
264287
}
@@ -268,14 +291,20 @@ class Teos10Eos {
268291
Real Lo = -15.0;
269292
Real Hi = 95.0;
270293

271-
if (P < 500.0) {
272-
Lo = calcCtFreezing(Sa, P, 0.0_Real);
273-
} else if (P < 6500.0) {
274-
Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real);
275-
Hi = 31.66666666666667 - P * 3.333333333333334e-3;
276-
} else {
277-
Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real);
278-
Hi = 10.0;
294+
if (EosLimChoice == EosLimits::Funnel) {
295+
if (P < 500.0) {
296+
Lo = calcCtFreezing(Sa, P, 0.0_Real);
297+
} else if (P < 6500.0) {
298+
Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real);
299+
Hi = 31.66666666666667 - P * 3.333333333333334e-3;
300+
} else {
301+
Lo = calcCtFreezing(Sa, 500.0_Real, 0.0_Real);
302+
Hi = 10.0;
303+
}
304+
}
305+
if (EosLimChoice == EosLimits::Cube) {
306+
Lo = -2.0;
307+
Hi = 40.0;
279308
}
280309
return {Lo, Hi};
281310
}

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

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -300,9 +300,9 @@ int testEosTeos10Clamping() {
300300
TestEos->EosChoice = EosType::Teos10Eos;
301301

302302
/// Create and fill ocean state arrays
303-
Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLevels);
304-
Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLevels);
305-
Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLevels);
303+
Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers);
304+
Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers);
305+
Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers);
306306
/// Use Kokkos::deep_copy to fill the entire view with the ref value
307307
/// -------------------------------------
308308
/// Test with valid poly75t values first (no clamping)
@@ -319,7 +319,7 @@ int testEosTeos10Clamping() {
319319
int numMismatches = 0;
320320
Array2DReal SpecVol = TestEos->SpecVol;
321321
parallelReduce(
322-
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels},
322+
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers},
323323
KOKKOS_LAMBDA(int i, int j, int &localCount) {
324324
if (!isApprox(SpecVol(i, j), TeosClampValue1, RTol)) {
325325
localCount++;
@@ -353,7 +353,7 @@ int testEosTeos10Clamping() {
353353
numMismatches = 0;
354354
SpecVol = TestEos->SpecVol;
355355
parallelReduce(
356-
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels},
356+
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers},
357357
KOKKOS_LAMBDA(int i, int j, int &localCount) {
358358
if (!isApprox(SpecVol(i, j), TeosClampValue2, RTol)) {
359359
localCount++;
@@ -387,7 +387,7 @@ int testEosTeos10Clamping() {
387387
numMismatches = 0;
388388
SpecVol = TestEos->SpecVol;
389389
parallelReduce(
390-
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels},
390+
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers},
391391
KOKKOS_LAMBDA(int i, int j, int &localCount) {
392392
if (!isApprox(SpecVol(i, j), TeosClampValue3, RTol)) {
393393
localCount++;
@@ -421,7 +421,7 @@ int testEosTeos10Clamping() {
421421
numMismatches = 0;
422422
SpecVol = TestEos->SpecVol;
423423
parallelReduce(
424-
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLevels},
424+
"CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers},
425425
KOKKOS_LAMBDA(int i, int j, int &localCount) {
426426
if (!isApprox(SpecVol(i, j), TeosClampValue4, RTol)) {
427427
localCount++;

0 commit comments

Comments
 (0)