diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index b2b1a7c37b16..620f790cdc65 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -28,6 +28,7 @@ Omega: WindStress: InterpType: Isotropic VertCoord: + Density0: 1026.0 MovementWeightType: Uniform Tendencies: ThicknessFluxTendencyEnable: true @@ -40,7 +41,6 @@ Omega: ViscDel4: 1.2e11 DivFactor: 1.0 WindForcingTendencyEnable: false - Density0: 1026.0 BottomDragTendencyEnable: false BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true diff --git a/components/omega/src/CMakeLists.txt b/components/omega/src/CMakeLists.txt index 4e635fc29545..4ea242099586 100644 --- a/components/omega/src/CMakeLists.txt +++ b/components/omega/src/CMakeLists.txt @@ -12,6 +12,7 @@ target_include_directories( ${OMEGA_SOURCE_DIR}/src/timeStepping ${E3SM_ROOT}/share/timing ${E3SM_ROOT}/share/pacer + ${E3SM_ROOT}/share/util_cxx ${PIOC_SOURCE_DIR} ${PIOC_BINARY_DIR} ${Parmetis_INCLUDE_DIRS} diff --git a/components/omega/src/ocn/CustomTendencyTerms.cpp b/components/omega/src/ocn/CustomTendencyTerms.cpp index 27f19a25ceb3..ff9c17da6435 100644 --- a/components/omega/src/ocn/CustomTendencyTerms.cpp +++ b/components/omega/src/ocn/CustomTendencyTerms.cpp @@ -8,6 +8,7 @@ #include "CustomTendencyTerms.h" #include "Config.h" +#include "GlobalConstants.h" #include "TimeStepper.h" namespace OMEGA { @@ -79,11 +80,9 @@ void ManufacturedSolution::init() { R8 H0 = DefHorzMesh->BottomDepthH(0); // Define and compute common constants - R8 Grav = 9.80665_Real; // Gravity acceleration - R8 Pii = 3.141592653589793_Real; // Pi - R8 Kx = 2.0_Real * Pii / WavelengthX; // Wave in X-dir - R8 Ky = 2.0_Real * Pii / WavelengthY; // Wave in Y-dir - R8 AngFreq = sqrt(H0 * Grav * (Kx * Kx + Ky * Ky)); // Angular frequency + R8 Kx = TwoPi / WavelengthX; // Wave in X-dir + R8 Ky = TwoPi / WavelengthY; // Wave in Y-dir + R8 AngFreq = sqrt(H0 * Gravity * (Kx * Kx + Ky * Ky)); // Angular frequency // Assign constants for thickness tendency function ManufacturedThickTend.H0 = H0; @@ -93,7 +92,7 @@ void ManufacturedSolution::init() { ManufacturedThickTend.AngFreq = AngFreq; // Assign constants for velocity tendency function - ManufacturedVelTend.Grav = Grav; + ManufacturedVelTend.Grav = Gravity; ManufacturedVelTend.Eta0 = Amplitude; ManufacturedVelTend.Kx = Kx; ManufacturedVelTend.Ky = Ky; @@ -160,7 +159,7 @@ void ManufacturedSolution::ManufacturedVelocityTendency::operator()( Array1DReal YEdge = Mesh->YEdge; Array1DReal AngleEdge = Mesh->AngleEdge; - OMEGA_SCOPE(LocGrav, Grav); + OMEGA_SCOPE(LocGrav, Gravity); OMEGA_SCOPE(LocEta0, Eta0); OMEGA_SCOPE(LocKx, Kx); OMEGA_SCOPE(LocKy, Ky); diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index b9eb0f949c31..a5588017ec6f 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -12,6 +12,7 @@ #include "AuxiliaryState.h" #include "Config.h" +#include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" #include "OmegaKokkos.h" @@ -205,8 +206,7 @@ class Teos10Eos { KOKKOS_FUNCTION Real calcDelta(const Array2DReal &SpecVolPCoeffs, const I4 K, const Real P) const { - constexpr Real Pu = 1e4; - Real Pp = P / Pu; + Real Pp = P * Pa2Db; Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp + SpecVolPCoeffs(3, K)) * @@ -221,14 +221,13 @@ class Teos10Eos { /// Calculate reference profile for TEOS-10 KOKKOS_FUNCTION Real calcRefProfile(Real P) const { - constexpr Real Pu = 1e4; constexpr Real V00 = -4.4015007269e-05; constexpr Real V01 = 6.9232335784e-06; constexpr Real V02 = -7.5004675975e-07; constexpr Real V03 = 1.7009109288e-08; constexpr Real V04 = -1.6884162004e-08; constexpr Real V05 = 1.9613503930e-09; - Real Pp = P / Pu; + Real Pp = P * Pa2Db; Real V0 = (((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) * @@ -244,9 +243,9 @@ class Teos10Eos { class LinearEos { public: /// Coefficients for LinearEos (overwritten by config file if set there) - Real DRhodT = -0.2; ///< Thermal expansion coefficient (kg m^-3 degC^-1) - Real DRhodS = 0.8; ///< Haline contraction coefficient (kg m^-3) - Real RhoT0S0 = 1000.0; ///< Reference density (kg m^-3) at (T,S)=(0,0) + Real DRhodT = -0.2; ///< Thermal expansion coefficient (kg m^-3 degC^-1) + Real DRhodS = 0.8; ///< Haline contraction coefficient (kg m^-3) + Real RhoT0S0 = RhoFw; ///< Reference density (kg m^-3) at (T,S)=(0,0) /// constructor declaration LinearEos(); diff --git a/components/omega/src/ocn/GlobalConstants.h b/components/omega/src/ocn/GlobalConstants.h new file mode 100644 index 000000000000..9d6fb25317bb --- /dev/null +++ b/components/omega/src/ocn/GlobalConstants.h @@ -0,0 +1,100 @@ +#ifndef OMEGA_GLOBALCONSTANTS_H +#define OMEGA_GLOBALCONSTANTS_H +//===-- ocn/GlobalConstants.h - Global Constants --------------------*- C++ +//-*-===// +// +/// \file +/// \brief Contains global constants for use across Omega +/// +/// This header defines constants to be used across Omega (to be replaced by a +/// shared E3SM constants file when that becomes available) +// +//===----------------------------------------------------------------------===// + +#include "DataTypes.h" +#include "pcd_const.h" +#include + +namespace OMEGA { + +// Mathematical constants +constexpr Real Pi = pcd::pi; // Pi (from Physical Constants Dictionary) +constexpr Real TwoPi = + 2.0 * pcd::pi; // 2*Pi (from Physical Constants Dictionary) +constexpr Real Radian = + pcd::radian; // Radians in one degree (from Physical Constants Dictionary) +constexpr Real Degree = + pcd::degree; // Degrees in one radian (from Physical Constants Dictionary) +constexpr Real SqrtTwo = + pcd::square_root_of_2; // Square root of 2 (from Physical Constants + // Dictionary) +constexpr Real SqrtThree = + pcd::square_root_of_3; // Square root of 3 (from Physical Constants + // Dictionary) + +// Earth constants +constexpr Real Gravity = + pcd::standard_acceleration_of_gravity; // Acceleration due to gravity ~ + // m/s^2 (from Physical Constants + // Dictionary) +constexpr Real CDay = 86400.0; // Seconds in a calendar day ~ sec +constexpr Real SDay = 86164.0; // Seconds in a sidereal day ~ sec +constexpr Real Omega = + 2.0 * Pi / SDay; // Angular velocity of the Earth ~ rad/sec +constexpr Real REarth = pcd::mean_radius; // Mean radius of the Earth ~ m (from + // Physical Constants Dictionary) + +/// Physical constants +constexpr Real TkTrip = + pcd::water_triple_point_temperature; // Triple point of fresh water ~ K + // (from Physical Constants Dictionary) +constexpr Real TkFrz = 273.15; // Freezing point of fresh water ~ K +constexpr Real TkFrzSw = TkFrz - 1.8; // Freezing point of seawater ~ K +constexpr Real RhoAir = 1.2; // Density of air ~ kg/m^3 +constexpr Real RhoFw = 1.000e3; // Density of fresh water ~ kg/m^3 +constexpr Real RhoSw = 1.026e3; // Density of seawater ~ kg/m^3 +constexpr Real RhoIce = 0.917e3; // Density of ice ~ kg/m^3 +constexpr Real CpAir = 1.005e3; // Specific heat capacity of air ~ J/(kg*K) +constexpr Real CpFw = + 4.188e3; // Specific heat capacity of fresh water ~ J/(kg*K) +constexpr Real CpSw = 3.996e3; // Specific heat capacity of seawater ~ J/(kg*K) +constexpr Real CpIce = 2.108e3; // Specific heat capacity of ice ~ J/(kg*K) +constexpr Real LatIce = 3.337e5; // Latent heat of fusion ~ J/kg +constexpr Real LatVap = 2.501e6; // Latent heat of vaporization ~ J/kg +constexpr Real LatSub = LatIce + LatVap; // Latent heat of sublimation ~ J/kg +constexpr Real CondIce = 2.1; // Universal gas constant ~ J/(mol*K) +constexpr Real OcnRefSal = 34.7; // Reference ocean salinity ~ psu +constexpr Real IceRefSal = 4.0; // Reference ice salinity ~ psu +constexpr Real Sound = 1.5e2; // Speed of sound ~ m/s +constexpr Real VonKar = 0.4; // Von Karman constant ~ dimensionless +constexpr Real Emiss = 1.0; // Emissivity ~ dimensionless +constexpr Real AtmRefP = + pcd::standard_atmosphere; // Reference atmospheric pressure ~ Pa (from + // Physical Constants Dictionary) + +// Conversion factors +constexpr Real Sec2Day = 1.0 / 86400.0; // Seconds to days +constexpr Real Day2Sec = 86400.0; // Days to seconds +constexpr Real Rad2Deg = + pcd::radian; // Radians to degrees (from Physical Constants Dictionary) +constexpr Real Deg2Rad = + pcd::degree; // Degrees to radians (from Physical Constants Dictionary) +constexpr Real Salt2PPt = 1000.0; // Salinity (kg/kg) to parts per thousand +constexpr Real PPt2Salt = 1.0e-3; // Parts per thousand to salinity (kg/kg) +constexpr Real Mass2Sv = 1.0e-12; // Mass flux (kg/s) to Sverdrup +constexpr Real Heat2Pw = 4.186e-15; // Heat flux (W) to Petawatt +constexpr Real Salt2SvPpt = 1.0e-9; // Salt flux (kg/s) to Sv*ppt +constexpr Real Salt2MmDay = 3.1536e+5; // Salt flux to water flux (mm/day) +constexpr Real Db2Pa = 1.0e4; // Decibar to Pascal +constexpr Real Pa2Db = 1.0e-4; // Pascal to Decibar +constexpr Real Cm2M = 1.0e-2; // Centimeters to meters +constexpr Real M2Cm = 1.0e2; // Meters to centimeters +constexpr Real HFluxFac = + 1.0 / (RhoSw * CpSw); // Heat flux (W/m^2) to temp flux (C*m/s) +constexpr Real FwFluxFac = 1.e-6; // Fw flux (kg/m^2/s) to salt((msu/psu)*m/s) +constexpr Real SaltFac = + -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s) +constexpr Real SFluxFac = 1.0; // Salt flux (kg/m^2/s) to salt flux (msu*m/s) + +} // namespace OMEGA +#endif diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 7d3efab243bd..9a0e40648402 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -179,9 +179,6 @@ void Tendencies::readTendConfig( CHECK_ERROR_ABORT( Err, "Tendencies: WindForcingTendencyEnable not found in TendConfig"); - Err += TendConfig->get("Density0", this->WindForcing.SaltWaterDensity); - CHECK_ERROR_ABORT(Err, "Tendencies: Density0 not found in TendConfig"); - Err += TendConfig->get("BottomDragTendencyEnable", this->BottomDrag.Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: BottomDragTendencyEnable not found in TendConfig"); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index aa4570d1f297..204cf7a65820 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -46,7 +46,7 @@ VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh) MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {} WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh) - : Enabled(false), EdgeMask(Mesh->EdgeMask) {} + : Enabled(false), LocRhoSw(RhoSw), EdgeMask(Mesh->EdgeMask) {} BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index a6a07574a931..bb718b12f05c 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -11,6 +11,7 @@ //===----------------------------------------------------------------------===// #include "AuxiliaryState.h" +#include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" #include "OceanState.h" @@ -166,14 +167,13 @@ class SSHGradOnEdge { for (int KVec = 0; KVec < VecLength; ++KVec) { const I4 K = KStart + KVec; - Tend(IEdge, K) -= EdgeMask(IEdge, K) * Grav * + Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity * (SshCell(ICell1, K) - SshCell(ICell0, K)) * InvDcEdge; } } private: - Real Grav = 9.80665_Real; Array2DI4 CellsOnEdge; Array1DReal DcEdge; Array2DReal EdgeMask; @@ -281,7 +281,7 @@ class VelocityHyperDiffOnEdge { class WindForcingOnEdge { public: bool Enabled; - Real SaltWaterDensity; + Real LocRhoSw; /// constructor declaration WindForcingOnEdge(const HorzMesh *Mesh); @@ -296,7 +296,7 @@ class WindForcingOnEdge { const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, K); Tend(IEdge, K) += EdgeMask(IEdge, K) * InvThickEdge * - NormalStressEdge(IEdge) / SaltWaterDensity; + NormalStressEdge(IEdge) / LocRhoSw; } } diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 7f3133aa22fd..771e8485ae91 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -10,6 +10,7 @@ #include "VertCoord.h" #include "Dimension.h" #include "Field.h" +#include "GlobalConstants.h" #include "IO.h" #include "IOStream.h" #include "OmegaKokkos.h" @@ -216,14 +217,14 @@ void VertCoord::completeSetup(Config *Options //< [in] configuration options MinLayerCellH = createHostMirrorCopy(MinLayerCell); BottomDepthH = createHostMirrorCopy(BottomDepth); - // Fetch reference desnity from Config - Config TendConfig("Tendencies"); + // Fetch reference density from Config + Config VCoordConfig("VertCoord"); Err.reset(); - Err += Options->get(TendConfig); - CHECK_ERROR_ABORT(Err, "VertCoord: Tendencies group not found in Config"); + Err += Options->get(VCoordConfig); + CHECK_ERROR_ABORT(Err, "VertCoord: VertCoord group not found in Config"); - Err += TendConfig.get("Density0", Rho0); - CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in TendConfig"); + Err += VCoordConfig.get("Density0", Rho0); + CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in VertCoord"); } // end completeSetup @@ -656,8 +657,6 @@ void VertCoord::computePressure( const Array1DReal &SurfacePressure // [in] surface pressure ) { - Real Gravity = 9.80616_Real; // gravitationl acceleration - OMEGA_SCOPE(LocRho0, Rho0); OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); @@ -745,8 +744,6 @@ void VertCoord::computeGeopotential( const Array1DReal &SelfAttractionLoading // [in] self attraction and loading ) { - Real Gravity = 9.80616_Real; // gravitationl acceleration - OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); OMEGA_SCOPE(LocGeopotMid, GeopotentialMid); @@ -784,8 +781,6 @@ void VertCoord::computeGeopotential( // reductions and a parallel_for over the active layers within a column. void VertCoord::computeTargetThickness() { - Real Gravity = 9.80616_Real; // gravitationl acceleration - OMEGA_SCOPE(LocRho0, Rho0); OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index 217c21de46df..119cb08fbc57 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -4,6 +4,7 @@ #include "Decomp.h" #include "Dimension.h" #include "Field.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -23,7 +24,7 @@ using namespace OMEGA; struct TestSetup { - Real Radius = 6371220; + Real Radius = REarth; KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 00505498439c..1b5e3f1368d7 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -3,6 +3,7 @@ #include "Decomp.h" #include "Dimension.h" #include "Field.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -26,10 +27,9 @@ using namespace OMEGA; struct TestSetupPlane { - Real Pi = M_PI; Real Lx = 1; - Real Ly = std::sqrt(3) / 2; + Real Ly = SqrtThree / 2; ErrorMeasures ExpectedKineticEnergyErrors = {0.00994439065100057897, 0.00703403756741667954}; @@ -69,54 +69,55 @@ struct TestSetupPlane { 0.0039954090464502795}; KOKKOS_FUNCTION Real layerThickness(Real X, Real Y) const { - return 2 + std::cos(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return 2 + std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real velocityX(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real velocityY(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real windStressX(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real windStressY(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real divergence(Real X, Real Y) const { - return 2 * Pi * (1. / Lx + 1. / Ly) * std::cos(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly); + return TwoPi * (1. / Lx + 1. / Ly) * std::cos(TwoPi * X / Lx) * + std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real relativeVorticity(Real X, Real Y) const { - return 2 * Pi * (-1. / Lx + 1. / Ly) * std::sin(2 * Pi * X / Lx) * - std::sin(2 * Pi * Y / Ly); + return TwoPi * (-1. / Lx + 1. / Ly) * std::sin(TwoPi * X / Lx) * + std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real velocityDel2X(Real X, Real Y) const { - return -4 * Pi * Pi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * velocityX(X, Y); + return -TwoPi * TwoPi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * velocityX(X, Y); } KOKKOS_FUNCTION Real velocityDel2Y(Real X, Real Y) const { - return -4 * Pi * Pi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * velocityY(X, Y); + return -TwoPi * TwoPi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * velocityY(X, Y); } KOKKOS_FUNCTION Real velocityDel2Div(Real X, Real Y) const { - return -4 * Pi * Pi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * divergence(X, Y); + return -TwoPi * TwoPi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * + divergence(X, Y); } KOKKOS_FUNCTION Real velocityDel2Curl(Real X, Real Y) const { - return -4 * Pi * Pi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * + return -TwoPi * TwoPi * (1 / (Lx * Lx) + 1 / (Ly * Ly)) * relativeVorticity(X, Y); } KOKKOS_FUNCTION Real planetaryVorticity(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real normalizedRelativeVorticity(Real X, Real Y) const { @@ -134,34 +135,34 @@ struct TestSetupPlane { } KOKKOS_FUNCTION Real tracer(Real X, Real Y) const { - return 2 - std::cos(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return 2 - std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real thickTracer(Real X, Real Y) const { - return 4 - std::pow(std::cos(2 * Pi * X / Lx), 2) * - std::pow(std::cos(2 * Pi * Y / Ly), 2); + return 4 - std::pow(std::cos(TwoPi * X / Lx), 2) * + std::pow(std::cos(TwoPi * Y / Ly), 2); } KOKKOS_FUNCTION Real del2Tracer(Real X, Real Y) const { - return 2 * Pi * Pi * - (4 * (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly) + - std::pow(std::cos(2 * Pi * X / Lx), 2) * + return TwoPi * Pi * + (4 * (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(TwoPi * X / Lx) * + std::cos(TwoPi * Y / Ly) + + std::pow(std::cos(TwoPi * X / Lx), 2) * (1 / Lx / Lx + - (2 / Ly / Ly + 1 / Lx / Lx) * std::cos(4 * Pi * Y / Ly)) - - (2 / Lx / Lx) * std::pow(std::sin(2 * Pi * X / Lx), 2) * - std::pow(std::cos(2 * Pi * Y / Ly), 2)); + (2 / Ly / Ly + 1 / Lx / Lx) * std::cos(2 * TwoPi * Y / Ly)) - + (2 / Lx / Lx) * std::pow(std::sin(TwoPi * X / Lx), 2) * + std::pow(std::cos(TwoPi * Y / Ly), 2)); } }; struct TestSetupSphere { - Real Radius = 6371220; + Real Radius = REarth; ErrorMeasures ExpectedKineticEnergyErrors = {0.0143579382532765844, 0.00681096618897046764}; - ErrorMeasures ExpectedVelocityDivErrors = {0.0136595773989793799, - 0.00367052484586382699}; + ErrorMeasures ExpectedVelocityDivErrors = {0.013652414501663887, + 0.00369043159835992}; ErrorMeasures ExpectedFluxThickErrors = {0.0159821090867812224, 0.010364511516135164}; @@ -175,22 +176,22 @@ struct TestSetupSphere { ErrorMeasures ExpectedNormPlanetVortVertexErrors = {0.00451268952953497778, 0.00101771171197261793}; - ErrorMeasures ExpectedNormRelVortEdgeErrors = {0.0125376497261775952, - 0.00307521304930552519}; + ErrorMeasures ExpectedNormRelVortEdgeErrors = {0.01255832726488633, + 0.003093726580373857}; ErrorMeasures ExpectedNormPlanetVortEdgeErrors = {0.00495174534686814403, 0.000855432390947949515}; - ErrorMeasures ExpectedDel2Errors = {0.00360406641962622652, - 0.00313406628499444213}; + ErrorMeasures ExpectedDel2Errors = {0.0036327578533323366, + 0.003161584491352198}; ErrorMeasures ExpectedDel2DivErrors = {0.0177782108439020134, - 0.00751922684420262138}; - ErrorMeasures ExpectedDel2RelVortErrors = {0.0915578492503972413, + 0.0075764886471247385}; + ErrorMeasures ExpectedDel2RelVortErrors = {0.09148975765556343, 0.0246736311927726465}; ErrorMeasures ExpectedHTracerErrors = {0.01603249913425972, 0.00546762028673672059}; ErrorMeasures ExpectedDel2TracerErrors = {0.0081206665417422382, - 0.004917863312407276}; + 0.004969575978774801}; ErrorMeasures ExpectedNormalStressErrors = {0.0038588958862868362, 0.003813760171030077}; diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index ca57e8e3ed8e..91010ce4094a 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -3,6 +3,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -22,13 +23,12 @@ using namespace OMEGA; // operator tests together with exact values of operators for computing errors // and expected error values struct TestSetupPlane { - Real Pi = M_PI; // lengths of periodic planar mesh // TODO: get this from the horizontal mesh once it supports periodic planar // meshes Real Lx = 1; - Real Ly = std::sqrt(3) / 2; + Real Ly = SqrtThree / 2; ErrorMeasures ExpectedDivErrors = {0.00124886886594427027, 0.00124886886590974385}; @@ -44,49 +44,47 @@ struct TestSetupPlane { 0.004200067675522098}; KOKKOS_FUNCTION Real exactScalar(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactGradScalarX(Real X, Real Y) const { - return 2 * Pi / Lx * std::cos(2 * Pi * X / Lx) * - std::sin(2 * Pi * Y / Ly); + return TwoPi / Lx * std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactGradScalarY(Real X, Real Y) const { - return 2 * Pi / Ly * std::sin(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly); + return TwoPi / Ly * std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactVecX(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactVecY(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactDivVec(Real X, Real Y) const { - return 2 * Pi * (1. / Lx + 1. / Ly) * std::cos(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly); + return TwoPi * (1. / Lx + 1. / Ly) * std::cos(TwoPi * X / Lx) * + std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real exactCurlVec(Real X, Real Y) const { - return 2 * Pi * (-1. / Lx + 1. / Ly) * std::sin(2 * Pi * X / Lx) * - std::sin(2 * Pi * Y / Ly); + return TwoPi * (-1. / Lx + 1. / Ly) * std::sin(TwoPi * X / Lx) * + std::sin(TwoPi * Y / Ly); } }; struct TestSetupSphere1 { // radius of spherical mesh // TODO: get this from the mesh - Real Radius = 6371220; - - ErrorMeasures ExpectedDivErrors = {0.013659577398978353, - 0.00367052484586382743}; - ErrorMeasures ExpectedGradErrors = {0.00187912292540628936, - 0.00149841802817334306}; - ErrorMeasures ExpectedCurlErrors = {0.0271404735181308317, - 0.025202316610921989}; + Real Radius = REarth; + + ErrorMeasures ExpectedDivErrors = {0.013652414501664885, + 0.0036904315983599676}; + ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, + 0.0015218320661105687}; + ErrorMeasures ExpectedCurlErrors = {0.02713957370128636, + 0.025202095212756463}; ErrorMeasures ExpectedReconErrors = {0.0206375134079833517, 0.00692590524910695858}; ErrorMeasures ExpectedAnisoInterpErrors = {0.0024015775047603197, @@ -129,7 +127,7 @@ struct TestSetupSphere1 { struct TestSetupSphere2 { // radius of spherical mesh // TODO: get this from the mesh - Real Radius = 6371220; + Real Radius = REarth; ErrorMeasures ExpectedDivErrors = {1.37734693033362766e-10, 0.000484370621558727582}; diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 2785876987ae..119af9221a26 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -6,6 +6,7 @@ #include "Dimension.h" #include "Error.h" #include "Field.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -24,7 +25,7 @@ using namespace OMEGA; struct TestSetup { - Real Radius = 6371220; + Real Radius = REarth; KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index e537e9e6341a..029a55bd5e0d 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -19,6 +19,7 @@ #include "Decomp.h" #include "Dimension.h" #include "Error.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -35,10 +36,9 @@ using namespace OMEGA; struct TestSetupPlane { - Real Pi = M_PI; Real Lx = 1; - Real Ly = std::sqrt(3) / 2; + Real Ly = SqrtThree / 2; ErrorMeasures ExpectedDivErrors = {0.00124886886594453264, 0.00124886886590977139}; @@ -59,52 +59,50 @@ struct TestSetupPlane { 0.01000133508329411}; KOKKOS_FUNCTION Real vectorX(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real vectorY(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real divergence(Real X, Real Y) const { - return 2 * Pi * (1. / Lx + 1. / Ly) * std::cos(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly); + return TwoPi * (1. / Lx + 1. / Ly) * std::cos(TwoPi * X / Lx) * + std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real scalar(Real X, Real Y) const { - return std::sin(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::sin(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real gradX(Real X, Real Y) const { - return 2 * Pi / Lx * std::cos(2 * Pi * X / Lx) * - std::sin(2 * Pi * Y / Ly); + return TwoPi / Lx * std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real gradY(Real X, Real Y) const { - return 2 * Pi / Ly * std::sin(2 * Pi * X / Lx) * - std::cos(2 * Pi * Y / Ly); + return TwoPi / Ly * std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real curl(Real X, Real Y) const { - return 2 * Pi * (-1. / Lx + 1. / Ly) * std::sin(2 * Pi * X / Lx) * - std::sin(2 * Pi * Y / Ly); + return TwoPi * (-1. / Lx + 1. / Ly) * std::sin(TwoPi * X / Lx) * + std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real laplaceVecX(Real X, Real Y) const { - return -4 * Pi * Pi * (1. / Lx / Lx + 1. / Ly / Ly) * - std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return -TwoPi * TwoPi * (1. / Lx / Lx + 1. / Ly / Ly) * + std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real laplaceVecY(Real X, Real Y) const { - return -4 * Pi * Pi * (1. / Lx / Lx + 1. / Ly / Ly) * - std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return -TwoPi * TwoPi * (1. / Lx / Lx + 1. / Ly / Ly) * + std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real layerThick(Real X, Real Y) const { - return 2. + std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return 2. + std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real planetaryVort(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real normRelVort(Real X, Real Y) const { @@ -116,40 +114,40 @@ struct TestSetupPlane { } KOKKOS_FUNCTION Real tracerFluxDiv(Real X, Real Y) const { - return (2 * Pi / (Lx * Ly)) * - (std::cos(2 * Pi * X / Lx) * - (2 * (Lx + Ly) * std::cos(2 * Pi * Y / Ly) + - (Lx + 2 * Ly) * std::sin(2 * Pi * X / Lx) * - std::pow(std::cos(2 * Pi * Y / Ly), 2) - - Lx * std::sin(2 * Pi * X / Lx) * - std::pow(std::sin(2 * Pi * Y / Ly), 2))); + return (TwoPi / (Lx * Ly)) * + (std::cos(TwoPi * X / Lx) * + (2 * (Lx + Ly) * std::cos(TwoPi * Y / Ly) + + (Lx + 2 * Ly) * std::sin(TwoPi * X / Lx) * + std::pow(std::cos(TwoPi * Y / Ly), 2) - + Lx * std::sin(TwoPi * X / Lx) * + std::pow(std::sin(TwoPi * Y / Ly), 2))); } KOKKOS_FUNCTION Real scalarA(Real X, Real Y) const { - return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + return std::cos(TwoPi * X / Lx) * std::sin(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real scalarB(Real X, Real Y) const { - return 2. + std::cos(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); + return 2. + std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } KOKKOS_FUNCTION Real tracerDiff(Real X, Real Y) const { - return -4 * Pi * Pi * std::sin(2 * Pi * Y / Ly) * - (2 * (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(2 * Pi * X / Lx) + + return -TwoPi * TwoPi * std::sin(TwoPi * Y / Ly) * + (2 * (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(TwoPi * X / Lx) + (1 / Ly / Ly + - (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(4 * Pi * X / Lx)) * - std::cos(2 * Pi * Y / Ly)); + (1 / Lx / Lx + 1 / Ly / Ly) * std::cos(2 * TwoPi * X / Lx)) * + std::cos(TwoPi * Y / Ly)); } KOKKOS_FUNCTION Real scalarC(Real X, Real Y) const { - return std::pow(std::cos(2 * Pi * X / Lx), 2) - - std::pow(std::sin(2 * Pi * Y / Ly), 2); + return std::pow(std::cos(TwoPi * X / Lx), 2) - + std::pow(std::sin(TwoPi * Y / Ly), 2); } KOKKOS_FUNCTION Real tracerHyperDiff(Real X, Real Y) const { - return -8 * Pi * Pi * - (std::cos(4 * Pi * X / Lx) / Lx / Lx + - std::cos(4 * Pi * Y / Ly) / Ly / Ly); + return -2 * TwoPi * TwoPi * + (std::cos(2 * TwoPi * X / Lx) / Lx / Lx + + std::cos(2 * TwoPi * Y / Ly) / Ly / Ly); } KOKKOS_FUNCTION Real windForcingX(Real X, Real Y, @@ -181,24 +179,22 @@ struct TestSetupPlane { struct TestSetupSphere { // radius of spherical mesh // TODO: get this from the mesh - Real Radius = 6371220; + Real Radius = REarth; - Real Pi = M_PI; - - ErrorMeasures ExpectedDivErrors = {0.0136595773989796766, - 0.00367052484586384131}; + ErrorMeasures ExpectedDivErrors = {0.013652414501664885, + 0.0036904315983599676}; ErrorMeasures ExpectedPVErrors = {0.0219217796608757037, 0.0122537418367830303}; - ErrorMeasures ExpectedGradErrors = {0.00187912292540623471, - 0.00149841802817334935}; - ErrorMeasures ExpectedLaplaceErrors = {0.281930203304510130, - 0.270530313560271740}; - ErrorMeasures ExpectedTrHAdvErrors = {0.0132310202299444034, - 0.0038523368564029538}; - ErrorMeasures ExpectedTrDel2Errors = {0.0486107109846934185, - 0.00507514214194892694}; - ErrorMeasures ExpectedTrDel4Errors = {0.000819552466009620408, - 0.00064700084412871962}; + ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, + 0.0015218320661105702}; + ErrorMeasures ExpectedLaplaceErrors = {0.28193638497826856, + 0.270546491554748}; + ErrorMeasures ExpectedTrHAdvErrors = {0.013227657020868148, + 0.0038723934782890863}; + ErrorMeasures ExpectedTrDel2Errors = {0.04865718541236144, + 0.005105510870642706}; + ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, + 0.0007118574326665881}; ErrorMeasures ExpectedWindForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, 0.0014897009917655022}; @@ -524,8 +520,8 @@ int testSSHGrad(int NVertLayers, Real RTol) { Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { - VecField[0] = -9.80665_Real * Setup.gradX(X, Y); - VecField[1] = -9.80665_Real * Setup.gradY(X, Y); + VecField[0] = -Gravity * Setup.gradX(X, Y); + VecField[1] = -Gravity * Setup.gradY(X, Y); }, ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); @@ -738,7 +734,7 @@ int testWindForcing(int NVertLayers) { Array2DReal NumWindForcing("NumWindForcing", Mesh->NEdgesOwned, NVertLayers); WindForcingOnEdge WindForcingOnE(Mesh); - WindForcingOnE.SaltWaterDensity = SaltWaterDensity; + WindForcingOnE.LocRhoSw = SaltWaterDensity; parallelFor( {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) { diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp index 827c288f0e4a..bee71c70481c 100644 --- a/components/omega/test/ocn/VertCoordTest.cpp +++ b/components/omega/test/ocn/VertCoordTest.cpp @@ -12,6 +12,7 @@ #include "Decomp.h" #include "Dimension.h" #include "Error.h" +#include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -107,8 +108,7 @@ int main(int argc, char *argv[]) { /// Initialize layer thickness and surface pressure so that resulting /// interface pressure is the number of layers above plus one - Real Gravity = 9.80616_Real; - Real Rho0 = DefVertCoord->Rho0; + Real Rho0 = DefVertCoord->Rho0; parallelFor( {NCellsAll}, KOKKOS_LAMBDA(int ICell) { SurfacePressure(ICell) = 1.0_Real;