diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 5e2a8fafc172..0ed03f37fda5 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -16,6 +16,8 @@ Omega: Advection: FluxThicknessType: Center FluxTracerType: Center + WindStress: + InterpType: Isotropic Tendencies: ThicknessFluxTendencyEnable: true PVTendencyEnable: true @@ -26,6 +28,10 @@ Omega: VelHyperDiffTendencyEnable: true ViscDel4: 1.2e11 DivFactor: 1.0 + WindForcingTendencyEnable: false + Density0: 1026.0 + BottomDragTendencyEnable: false + BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true TracerDiffTendencyEnable: true EddyDiff2: 10.0 @@ -91,7 +97,7 @@ Omega: Contents: - Tracers - State - - SshCellDefault + - SshCell Highfreq: UsePointerFile: false Filename: ocn.hifreq.$Y-$M diff --git a/components/omega/doc/devGuide/AuxiliaryVariables.md b/components/omega/doc/devGuide/AuxiliaryVariables.md index 25701c0ef275..f6802f9a9383 100644 --- a/components/omega/doc/devGuide/AuxiliaryVariables.md +++ b/components/omega/doc/devGuide/AuxiliaryVariables.md @@ -101,3 +101,6 @@ The following auxiliary variable groups are currently implemented: || Del2RelVortVertex || | TracerAuxVars | HTracersEdge | Center or Upwind| || Del2TracersCell || +| WindForcingAuxVars | ZonalStressCell || +|| MeridStressCell || +|| NormalStressEdge || diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index c47b6d2c0dfc..8280abc1a054 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -38,3 +38,5 @@ implemented: - `TracerHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` +- `WindForcingOnEdge` +- `BottomDragOnEdge` diff --git a/components/omega/doc/userGuide/AuxiliaryVariables.md b/components/omega/doc/userGuide/AuxiliaryVariables.md index 0a0fabc98fb8..4f6fbc18c616 100644 --- a/components/omega/doc/userGuide/AuxiliaryVariables.md +++ b/components/omega/doc/userGuide/AuxiliaryVariables.md @@ -29,3 +29,6 @@ The following auxiliary variables are currently available: | VelDel2RelVortVertex | curl of laplacian of horizontal velocity on cells | HTracersEdge | thickness-weighted tracers used for fluxes through edges. May be centered, upwinded or a combination of the two | Del2TracersCell | laplacian of tracers on cells +| ZonalStressCell | zonal component of wind stress on cells +| MeridStressCell | meridional component of wind stress on cells +| NormalStressEdge | normal component of wind stress on edge diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 37339be6a4d4..69ad8437d210 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -17,6 +17,8 @@ tendency terms are currently implemented: | TracerHorzAdvOnCell | horizontal advection of thickness-weighted tracers | TracerDiffOnCell | horizontal diffusion of thickness-weighted tracers | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers +| WindForcingOnEdge | forcing by wind stress, defined on edges +| BottomDragOnEdge | bottom drag, defined on edges Among the internal data stored by each functor is a `bool` which can enable or disable the contribution of that particular term to the tendency. These flags @@ -45,3 +47,6 @@ the currently available tendency terms: | | EddyDiff2 | horizontal diffusion coefficient | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term | | EddyDiff4 | biharmonic horizontal mixing coeffienct for tracers +| WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term +| BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term +| | BottomDragCoeff | bottom drag coefficient diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index e9427adf3260..e19275c03e74 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -10,15 +10,21 @@ AuxiliaryState *AuxiliaryState::DefaultAuxState = nullptr; std::map> AuxiliaryState::AllAuxStates; +static std::string stripDefault(const std::string &Name) { + return Name != "Default" ? Name : ""; +} + // Constructor. Constructs the member auxiliary variables and registers their // fields with IOStreams AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, - int NVertLevels, int NTracers) - : Mesh(Mesh), Name(Name), KineticAux(Name, Mesh, NVertLevels), - LayerThicknessAux(Name, Mesh, NVertLevels), - VorticityAux(Name, Mesh, NVertLevels), - VelocityDel2Aux(Name, Mesh, NVertLevels), - TracerAux(Name, Mesh, NVertLevels, NTracers) { + Halo *MeshHalo, int NVertLevels, int NTracers) + : Mesh(Mesh), MeshHalo(MeshHalo), Name(stripDefault(Name)), + KineticAux(stripDefault(Name), Mesh, NVertLevels), + LayerThicknessAux(stripDefault(Name), Mesh, NVertLevels), + VorticityAux(stripDefault(Name), Mesh, NVertLevels), + VelocityDel2Aux(stripDefault(Name), Mesh, NVertLevels), + WindForcingAux(stripDefault(Name), Mesh), + TracerAux(stripDefault(Name), Mesh, NVertLevels, NTracers) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -32,6 +38,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, LayerThicknessAux.registerFields(GroupName, AuxMeshName); VorticityAux.registerFields(GroupName, AuxMeshName); VelocityDel2Aux.registerFields(GroupName, AuxMeshName); + WindForcingAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); } @@ -42,6 +49,7 @@ AuxiliaryState::~AuxiliaryState() { LayerThicknessAux.unregisterFields(); VorticityAux.unregisterFields(); VelocityDel2Aux.unregisterFields(); + WindForcingAux.unregisterFields(); TracerAux.unregisterFields(); int Err = FieldGroup::destroy(GroupName); @@ -64,6 +72,7 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocVorticityAux, VorticityAux); OMEGA_SCOPE(LocVelocityDel2Aux, VelocityDel2Aux); + OMEGA_SCOPE(LocWindForcingAux, WindForcingAux); parallelFor( "vertexAuxState1", {Mesh->NVerticesAll, NChunks}, @@ -82,7 +91,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, const auto &RelVortVertex = VorticityAux.RelVortVertex; parallelFor( - "edgeAuxState1", {Mesh->NEdgesAll, NChunks}, + "edgeAuxState1", {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + LocWindForcingAux.computeVarsOnEdge(IEdge); + }); + + parallelFor( + "edgeAuxState2", {Mesh->NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { LocVorticityAux.computeVarsOnEdge(IEdge, KChunk); LocLayerThicknessAux.computeVarsOnEdge(IEdge, KChunk, LayerThickCell, @@ -153,8 +167,8 @@ void AuxiliaryState::computeAll(const OceanState *State, // Create a non-default auxiliary state AuxiliaryState *AuxiliaryState::create(const std::string &Name, - const HorzMesh *Mesh, int NVertLevels, - const int NTracers) { + const HorzMesh *Mesh, Halo *MeshHalo, + int NVertLevels, const int NTracers) { if (AllAuxStates.find(Name) != AllAuxStates.end()) { LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it " "already exists", @@ -162,7 +176,8 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, return nullptr; } - auto *NewAuxState = new AuxiliaryState(Name, Mesh, NVertLevels, NTracers); + auto *NewAuxState = + new AuxiliaryState(Name, Mesh, MeshHalo, NVertLevels, NTracers); AllAuxStates.emplace(Name, NewAuxState); return NewAuxState; @@ -172,12 +187,13 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, // initialized. void AuxiliaryState::init() { const HorzMesh *DefMesh = HorzMesh::getDefault(); + Halo *DefHalo = Halo::getDefault(); int NVertLevels = DefMesh->NVertLevels; int NTracers = Tracers::getNumTracers(); - AuxiliaryState::DefaultAuxState = - AuxiliaryState::create("Default", DefMesh, NVertLevels, NTracers); + AuxiliaryState::DefaultAuxState = AuxiliaryState::create( + "Default", DefMesh, DefHalo, NVertLevels, NTracers); Config *OmegaConfig = Config::getOmegaConfig(); DefaultAuxState->readConfigOptions(OmegaConfig); @@ -248,6 +264,37 @@ void AuxiliaryState::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("AuxiliaryState: Unknown FluxTracerType requested"); } + + Config WindStressConfig("WindStress"); + Err += OmegaConfig->get(WindStressConfig); + + std::string WindStressInterpTypeStr; + Err += WindStressConfig.get("InterpType", WindStressInterpTypeStr); + CHECK_ERROR_ABORT( + Err, "AuxiliaryState: InterpType not found in WindStressConfig"); + + if (WindStressInterpTypeStr == "Isotropic") { + this->WindForcingAux.InterpChoice = InterpCellToEdgeOption::Isotropic; + } else if (WindStressInterpTypeStr == "Anisotropic") { + this->WindForcingAux.InterpChoice = InterpCellToEdgeOption::Anisotropic; + } else { + ABORT_ERROR("AuxiliaryState: Unknown InterpType requested"); + } } +//------------------------------------------------------------------------------ +// Perform auxiliary state halo exchange +// Note that only non-computed auxiliary variables needs to be exchanged +I4 AuxiliaryState::exchangeHalo() { + I4 Err = 0; + + Err += + MeshHalo->exchangeFullArrayHalo(WindForcingAux.ZonalStressCell, OnCell); + Err += + MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell); + + return Err; + +} // end exchangeHalo + } // namespace OMEGA diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 13809387ee1b..20eb40caa694 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -3,6 +3,7 @@ #include "Config.h" #include "DataTypes.h" +#include "Halo.h" #include "HorzMesh.h" #include "OceanState.h" #include "Tracers.h" @@ -11,6 +12,7 @@ #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" +#include "auxiliaryVars/WindForcingAuxVars.h" #include #include @@ -36,6 +38,7 @@ class AuxiliaryState { TracerAuxVars TracerAux; VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; + WindForcingAuxVars WindForcingAux; ~AuxiliaryState(); @@ -46,7 +49,7 @@ class AuxiliaryState { // Create a non-default auxiliary state static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh, - int NVertLevels, int NTracers); + Halo *MeshHalo, int NVertLevels, int NTracers); /// Get the default auxiliary state static AuxiliaryState *getDefault(); @@ -63,6 +66,9 @@ class AuxiliaryState { /// Read and set config options void readConfigOptions(Config *OmegaConfig); + /// Exchange halo + I4 exchangeHalo(); + // Compute all auxiliary variables needed for momentum equation void computeMomAux(const OceanState *State, int ThickTimeLevel, int VelTimeLevel) const; @@ -75,13 +81,14 @@ class AuxiliaryState { int TimeLevel) const; private: - AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, + AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, int NVertLevels, int NTracers); AuxiliaryState(const AuxiliaryState &) = delete; AuxiliaryState(AuxiliaryState &&) = delete; const HorzMesh *Mesh; + Halo *MeshHalo; static AuxiliaryState *DefaultAuxState; static std::map> AllAuxStates; }; diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index f650f4377a30..b0f13dbe8d59 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -21,4 +21,10 @@ TangentialReconOnEdge::TangentialReconOnEdge(HorzMesh const *Mesh) : NEdgesOnEdge(Mesh->NEdgesOnEdge), EdgesOnEdge(Mesh->EdgesOnEdge), WeightsOnEdge(Mesh->WeightsOnEdge) {} +InterpCellToEdge::InterpCellToEdge(const HorzMesh *Mesh) + : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), + CellsOnVertex(Mesh->CellsOnVertex), + KiteAreasOnVertex(Mesh->KiteAreasOnVertex), + VertexDegree(Mesh->VertexDegree) {} + } // namespace OMEGA diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index 9d26f255d5cd..e891da0bd228 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -132,5 +132,57 @@ class TangentialReconOnEdge { Array2DReal WeightsOnEdge; }; +enum class InterpCellToEdgeOption { Anisotropic, Isotropic }; + +class InterpCellToEdge { + public: + InterpCellToEdge(const HorzMesh *Mesh); + + KOKKOS_FUNCTION Real operator()(int IEdge, const Array1DReal &ArrayCell, + InterpCellToEdgeOption Option) const { + switch (Option) { + case InterpCellToEdgeOption::Anisotropic: + return interpolateAnisotropic(IEdge, ArrayCell); + case InterpCellToEdgeOption::Isotropic: + return interpolateIsotropic(IEdge, ArrayCell); + } + }; + + KOKKOS_FUNCTION Real + interpolateAnisotropic(int IEdge, const Array1DReal &ArrayCell) const { + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + + return 0.5_Real * (ArrayCell(JCell0) + ArrayCell(JCell1)); + }; + + KOKKOS_FUNCTION Real + interpolateIsotropic(int IEdge, const Array1DReal &ArrayCell) const { + + Real Accum = 0; + Real AreaAccum = 0; + for (int J = 0; J < 2; ++J) { + const int JVertex = VerticesOnEdge(IEdge, J); + for (int L = 0; L < VertexDegree; ++L) { + const Real KiteArea = KiteAreasOnVertex(JVertex, L); + const int LCell = CellsOnVertex(JVertex, L); + + Accum += ArrayCell(LCell) * KiteArea; + AreaAccum += KiteArea; + } + } + + const Real InvAreaAccum = 1._Real / AreaAccum; + return Accum * InvAreaAccum; + }; + + private: + Array2DI4 CellsOnEdge; + Array2DI4 VerticesOnEdge; + Array2DI4 CellsOnVertex; + Array2DReal KiteAreasOnVertex; + I4 VertexDegree; +}; + } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 19491a7bf0b0..e92ce5d833fe 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -157,13 +157,17 @@ int initOmegaModules(MPI_Comm Comm) { } } - // Update Halos and Host arrays with new state and tracer fields + // Update Halos and Host arrays with new state, auxiliary state, and tracer + // fields OceanState *DefState = OceanState::getDefault(); I4 CurTimeLevel = 0; DefState->exchangeHalo(CurTimeLevel); DefState->copyToHost(CurTimeLevel); + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); + DefAuxState->exchangeHalo(); + // Now update tracers - assume using same time level index Err = Tracers::exchangeHalo(CurTimeLevel); if (Err != 0) { diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index d8d8b72ccfcf..bf4949f43949 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -172,6 +172,27 @@ void Tendencies::readTendConfig( CHECK_ERROR_ABORT( Err, "Tendencies: TracerDiffTendencyEnable not found in TendConfig"); + Err += + TendConfig->get("WindForcingTendencyEnable", this->WindForcing.Enabled); + 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"); + + Err += TendConfig->get("BottomDragCoeff", this->BottomDrag.Coeff); + CHECK_ERROR_ABORT(Err, + "Tendencies: BottomDragCoeff not found in TendConfig"); + + Err += TendConfig->get("TracerHorzAdvTendencyEnable", + this->TracerHorzAdv.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: TracerHorzAdvTendencyEnable not found in TendConfig"); + if (this->TracerDiffusion.Enabled) { Err += TendConfig->get("EddyDiff2", this->TracerDiffusion.EddyDiff2); CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff2 not found in TendConfig"); @@ -200,8 +221,8 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies CustomTendencyType InCustomVelocityTend) : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh), KEGrad(Mesh), SSHGrad(Mesh), VelocityDiffusion(Mesh), VelocityHyperDiff(Mesh), - TracerHorzAdv(Mesh), TracerDiffusion(Mesh), TracerHyperDiff(Mesh), - CustomThicknessTend(InCustomThicknessTend), + BottomDrag(Mesh), TracerHorzAdv(Mesh), TracerDiffusion(Mesh), + TracerHyperDiff(Mesh), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend) { // Tendency arrays @@ -280,9 +301,13 @@ void Tendencies::computeVelocityTendenciesOnly( OMEGA_SCOPE(LocSSHGrad, SSHGrad); OMEGA_SCOPE(LocVelocityDiffusion, VelocityDiffusion); OMEGA_SCOPE(LocVelocityHyperDiff, VelocityHyperDiff); + OMEGA_SCOPE(LocWindForcing, WindForcing); + OMEGA_SCOPE(LocBottomDrag, BottomDrag); deepCopy(LocNormalVelocityTend, 0); + const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + // Compute potential vorticity horizontal advection const Array2DReal &FluxLayerThickEdge = AuxState->LayerThicknessAux.FluxLayerThickEdge; @@ -340,6 +365,27 @@ void Tendencies::computeVelocityTendenciesOnly( }); } + // Compute wind forcing + const auto &NormalStressEdge = AuxState->WindForcingAux.NormalStressEdge; + const auto &MeanLayerThickEdge = + AuxState->LayerThicknessAux.MeanLayerThickEdge; + if (LocWindForcing.Enabled) { + parallelFor( + {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { + LocWindForcing(LocNormalVelocityTend, IEdge, KChunk, + NormalStressEdge, MeanLayerThickEdge); + }); + } + + // Compute bottom drag + if (LocBottomDrag.Enabled) { + parallelFor( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + LocBottomDrag(LocNormalVelocityTend, IEdge, NormalVelEdge, KECell, + MeanLayerThickEdge); + }); + } + if (CustomVelocityTend) { CustomVelocityTend(LocNormalVelocityTend, State, AuxState, ThickTimeLevel, VelTimeLevel, Time); diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index eae1744a514d..d27b7c08cf8c 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -62,6 +62,8 @@ class Tendencies { SSHGradOnEdge SSHGrad; VelocityDiffusionOnEdge VelocityDiffusion; VelocityHyperDiffOnEdge VelocityHyperDiff; + WindForcingOnEdge WindForcing; + BottomDragOnEdge BottomDrag; TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index cf9829e9c254..58a446560fa1 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -42,6 +42,10 @@ VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh) DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {} +BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh) + : Enabled(false), Coeff(0), CellsOnEdge(Mesh->CellsOnEdge), + NVertLevels(Mesh->NVertLevels) {} + TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 8cdd817a1847..3b645f497d02 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -271,6 +271,61 @@ class VelocityHyperDiffOnEdge { Array2DReal EdgeMask; }; +/// Wind forcing +class WindForcingOnEdge { + public: + bool Enabled; + Real SaltWaterDensity; + + /// The functor takes the edge index, vertical chunk index, and arrays for + /// normal wind stress and edge layer thickness, outputs tendency array + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, + const Array1DReal &NormalStressEdge, + const Array2DReal &LayerThickEdge) const { + if (KChunk == 0) { + const I4 K = 0; + + const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, K); + Tend(IEdge, K) += + InvThickEdge * NormalStressEdge(IEdge) / SaltWaterDensity; + } + } +}; + +/// Bottom drag +class BottomDragOnEdge { + public: + bool Enabled; + Real Coeff; + + /// constructor declaration + BottomDragOnEdge(const HorzMesh *Mesh); + + /// The functor takes the edge index and arrays for + /// horizontal velocity, kinetic energy, + /// and edge layer thickness, outputs tendency array + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, + const Array2DReal &NormalVelEdge, + const Array2DReal &KECell, + const Array2DReal &LayerThickEdge) const { + const I4 KBot = NVertLevels - 1; + + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + + const Real VelNormEdge = + Kokkos::sqrt(KECell(JCell0, KBot) + KECell(JCell1, KBot)); + + const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, KBot); + Tend(IEdge, KBot) -= + Coeff * VelNormEdge * InvThickEdge * NormalVelEdge(IEdge, KBot); + } + + private: + I4 NVertLevels; + Array2DI4 CellsOnEdge; +}; + // Tracer horizontal advection term class TracerHorzAdvOnCell { public: diff --git a/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.cpp new file mode 100644 index 000000000000..9f4fa52480a5 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.cpp @@ -0,0 +1,93 @@ +#include "WindForcingAuxVars.h" +#include "DataTypes.h" +#include "Field.h" + +#include + +namespace OMEGA { + +WindForcingAuxVars::WindForcingAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh) + : NormalStressEdge("NormalStressEdge" + AuxStateSuffix, Mesh->NEdgesSize), + ZonalStressCell("WindStressZonal" + AuxStateSuffix, Mesh->NCellsSize), + MeridStressCell("WindStressMeridional" + AuxStateSuffix, + Mesh->NCellsSize), + CellsOnEdge(Mesh->CellsOnEdge), AngleEdge(Mesh->AngleEdge), Interp(Mesh) { +} + +void WindForcingAuxVars::registerFields( + const std::string &AuxGroupName, // name of Auxiliary field group + const std::string &MeshName // name of horizontal mesh +) const { + + int Err = 0; // error flag for some calls + + // Create fields + const Real FillValue = -9.99e30; + int NDims = 1; + std::vector DimNames(NDims); + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + // Zonal wind stress + DimNames[0] = "NCells" + DimSuffix; + auto ZonalStressCellField = + Field::create(ZonalStressCell.label(), // field name + "zonal wind stress", // long name/describe + "N m^{-2}", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dim names + ); + + // Meridional wind stress + auto MeridStressCellField = + Field::create(MeridStressCell.label(), // field name + "meridional wind stress", // long Name or description + "N m^{-2}", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar used for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + // Add fields to FieldGroup + Err = FieldGroup::addFieldToGroup(ZonalStressCell.label(), AuxGroupName); + if (Err != 0) + LOG_ERROR("Error adding field {} to group {}", ZonalStressCell.label(), + AuxGroupName); + Err = FieldGroup::addFieldToGroup(MeridStressCell.label(), AuxGroupName); + if (Err != 0) + LOG_ERROR("Error adding field {} to group {}", MeridStressCell.label(), + AuxGroupName); + + // Attach data + Err = ZonalStressCellField->attachData(ZonalStressCell); + if (Err != 0) + LOG_ERROR("Error attaching data to field {}", ZonalStressCell.label()); + + Err = MeridStressCellField->attachData(MeridStressCell); + if (Err != 0) + LOG_ERROR("Error attaching data to field {}", MeridStressCell.label()); +} + +void WindForcingAuxVars::unregisterFields() const { + int Err = 0; + Err = Field::destroy(ZonalStressCell.label()); + if (Err != 0) + LOG_ERROR("Error destroying field {}", ZonalStressCell.label()); + Err = Field::destroy(MeridStressCell.label()); + if (Err != 0) + LOG_ERROR("Error destroying field {}", MeridStressCell.label()); +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h new file mode 100644 index 000000000000..b1546f441216 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h @@ -0,0 +1,42 @@ +#ifndef OMEGA_AUX_WIND_H +#define OMEGA_AUX_WIND_H + +#include "DataTypes.h" +#include "HorzMesh.h" +#include "HorzOperators.h" +#include "OmegaKokkos.h" + +#include + +namespace OMEGA { + +class WindForcingAuxVars { + public: + Array1DReal NormalStressEdge; + Array1DReal ZonalStressCell; + Array1DReal MeridStressCell; + InterpCellToEdgeOption InterpChoice; + + WindForcingAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh); + + KOKKOS_FUNCTION void computeVarsOnEdge(int IEdge) const { + const Real ZonalStressEdge = Interp(IEdge, ZonalStressCell, InterpChoice); + const Real MeridStressEdge = Interp(IEdge, MeridStressCell, InterpChoice); + + NormalStressEdge(IEdge) = + Kokkos::cos(AngleEdge(IEdge)) * ZonalStressEdge + + Kokkos::sin(AngleEdge(IEdge)) * MeridStressEdge; + } + + void registerFields(const std::string &AuxGroupName, + const std::string &MeshName) const; + void unregisterFields() const; + + private: + InterpCellToEdge Interp; + Array2DI4 CellsOnEdge; + Array1DReal AngleEdge; +}; + +} // namespace OMEGA +#endif diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index b7f1feaa9c7c..17051b3fd91d 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -317,6 +317,11 @@ add_omega_test( "-n 8;" single_precision ) +target_compile_definitions( + testTendencyTermsPlaneSinglePrecision.exe + PRIVATE + TENDENCYTERMS_TEST_PLANE +) add_omega_test( TEND_SPHERE_TEST diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index e54d9750f88b..3d42fbeaad28 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -63,19 +63,19 @@ int initState() { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickCell, Geom, Mesh, OnCell, NVertLevels); + LayerThickCell, Geom, Mesh, OnCell); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracerArray, Geom, Mesh, OnCell, NVertLevels, NTracers); + TracerArray, Geom, Mesh, OnCell); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { VecField[0] = Setup.velocityX(Lon, Lat); VecField[1] = Setup.velocityY(Lon, Lat); }, - NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::Yes, CartProjection::No); + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::Yes, + CartProjection::No); return Err; } @@ -145,8 +145,9 @@ int testAuxState() { } const auto *Mesh = HorzMesh::getDefault(); + auto *MeshHalo = Halo::getDefault(); // test creation of another auxiliary state - AuxiliaryState::create("AnotherAuxState", Mesh, 12, 3); + AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, 12, 3); // test retrievel of another if (AuxiliaryState::get("AnotherAuxState")) { diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 38915e14a118..f2abcc21df75 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -16,6 +16,7 @@ #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" +#include "auxiliaryVars/WindForcingAuxVars.h" #include "mpi.h" #include @@ -63,6 +64,9 @@ struct TestSetupPlane { ErrorMeasures ExpectedDel2TracerErrors = {0.0033346711042859123, 0.0029202923731303323}; + ErrorMeasures ExpectedNormalStressErrors = {0.0033910709836867704, + 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); } @@ -75,6 +79,14 @@ struct TestSetupPlane { return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); } + KOKKOS_FUNCTION Real windStressX(Real X, Real Y) const { + return std::cos(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); + } + + KOKKOS_FUNCTION Real windStressY(Real X, Real Y) const { + return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * 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); @@ -179,6 +191,9 @@ struct TestSetupSphere { ErrorMeasures ExpectedDel2TracerErrors = {0.0081206665417422382, 0.004917863312407276}; + ErrorMeasures ExpectedNormalStressErrors = {0.0038588958862868362, + 0.003813760171030077}; + KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); } @@ -191,6 +206,14 @@ struct TestSetupSphere { return -4 * std::sin(Lon) * std::cos(Lon) * std::pow(std::cos(Lat), 3) * std::sin(Lat); } + KOKKOS_FUNCTION Real windStressX(Real Lon, Real Lat) const { + return -4 * std::sin(Lon) * std::cos(Lon) * std::pow(std::cos(Lat), 3) * + std::sin(Lat); + } + + KOKKOS_FUNCTION Real windStressY(Real Lon, Real Lat) const { + return -std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); + } KOKKOS_FUNCTION Real relativeVorticity(Real Lon, Real Lat) const { return -4 * std::pow(std::cos(Lon), 2) * std::pow(std::cos(Lat), 2) * @@ -297,40 +320,21 @@ int initState(const Array2DReal &LayerThickCell, Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickCell, Geom, Mesh, OnCell, NVertLevels); + LayerThickCell, Geom, Mesh, OnCell); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { VecField[0] = Setup.velocityX(X, Y); VecField[1] = Setup.velocityY(X, Y); }, - NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh); // need to override FVertex with prescribed values - // cannot use setScalar because it doesn't support setting 1D arrays const auto &FVertex = Mesh->FVertex; - auto XVertex = createDeviceMirrorCopy(Mesh->XVertexH); - auto YVertex = createDeviceMirrorCopy(Mesh->YVertexH); - - auto LonVertex = createDeviceMirrorCopy(Mesh->LonVertexH); - auto LatVertex = createDeviceMirrorCopy(Mesh->LatVertexH); - - parallelFor( - {Mesh->NVerticesOwned}, KOKKOS_LAMBDA(int IVertex) { - if (Geom == Geometry::Planar) { - const Real XV = XVertex(IVertex); - const Real YV = YVertex(IVertex); - FVertex(IVertex) = Setup.planetaryVorticity(XV, YV); - } else { - const Real XV = LonVertex(IVertex); - const Real YV = LatVertex(IVertex); - FVertex(IVertex) = Setup.planetaryVorticity(XV, YV); - } - }); - - auto MyHalo = Halo::getDefault(); - Err += MyHalo->exchangeFullArrayHalo(FVertex, OnVertex); + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.planetaryVorticity(X, Y); }, + FVertex, Geom, Mesh, OnVertex); return Err; } @@ -348,15 +352,13 @@ int testKineticAuxVars(const Array2DReal &LayerThicknessCell, Mesh->NCellsOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.kineticEnergy(X, Y); }, - ExactKineticEnergyCell, Geom, Mesh, OnCell, NVertLevels, - ExchangeHalos::No); + ExactKineticEnergyCell, Geom, Mesh, OnCell, ExchangeHalos::No); Array2DReal ExactVelocityDivCell("ExactVelocityDivCell", Mesh->NCellsOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, - ExactVelocityDivCell, Geom, Mesh, OnCell, NVertLevels, - ExchangeHalos::No); + ExactVelocityDivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical result @@ -373,13 +375,13 @@ int testKineticAuxVars(const Array2DReal &LayerThicknessCell, ErrorMeasures KineticEnergyErrors; Err += computeErrors(KineticEnergyErrors, NumKineticEnergyCell, - ExactKineticEnergyCell, Mesh, OnCell, NVertLevels); + ExactKineticEnergyCell, Mesh, OnCell); Err += checkErrors("AuxVarsTest", "KineticEnergy", KineticEnergyErrors, Setup.ExpectedKineticEnergyErrors, RTol); ErrorMeasures VelocityDivErrors; Err += computeErrors(VelocityDivErrors, NumVelocityDivCell, - ExactVelocityDivCell, Mesh, OnCell, NVertLevels); + ExactVelocityDivCell, Mesh, OnCell); Err += checkErrors("AuxVarsTest", "VelocityDiv", VelocityDivErrors, Setup.ExpectedVelocityDivErrors, RTol); @@ -390,6 +392,57 @@ int testKineticAuxVars(const Array2DReal &LayerThicknessCell, return Err; } +int testWindForcingAuxVars(Real RTol) { + int Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + + // Compute exact result + + Array1DReal ExactNormalStressEdge("ExactNormalStressEdge", + Mesh->NEdgesOwned); + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.windStressX(X, Y); + VecField[1] = Setup.windStressY(X, Y); + }, + ExactNormalStressEdge, EdgeComponent::Normal, Geom, Mesh, + ExchangeHalos::No); + + WindForcingAuxVars WindForcingAux("", Mesh); + WindForcingAux.InterpChoice = InterpCellToEdgeOption::Anisotropic; + + // Set inputs + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.windStressX(X, Y); }, + WindForcingAux.ZonalStressCell, Geom, Mesh, OnCell); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.windStressY(X, Y); }, + WindForcingAux.MeridStressCell, Geom, Mesh, OnCell); + + // Compute numerical result + parallelFor( + {Mesh->NEdgesOwned}, + KOKKOS_LAMBDA(int IEdge) { WindForcingAux.computeVarsOnEdge(IEdge); }); + const auto &NumNormalStressEdge = WindForcingAux.NormalStressEdge; + + // Compute error measures and check error values + ErrorMeasures NormalStressErrors; + Err += computeErrors(NormalStressErrors, NumNormalStressEdge, + ExactNormalStressEdge, Mesh, OnEdge); + + Err += checkErrors("AuxVarsTest", "NormalStress", NormalStressErrors, + Setup.ExpectedNormalStressErrors, RTol); + + if (Err == 0) { + LOG_INFO("AuxVarsTest: WindForcingAuxVars PASS"); + } + + return Err; +} + int testLayerThicknessAuxVars(const Array2DReal &LayerThickCell, const Array2DReal &NormalVelEdge, Real RTol) { int Err = 0; @@ -402,7 +455,7 @@ int testLayerThicknessAuxVars(const Array2DReal &LayerThickCell, Array2DReal ExactThickEdge("ExactThickEdge", Mesh->NEdgesOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - ExactThickEdge, Geom, Mesh, OnEdge, NVertLevels, ExchangeHalos::No); + ExactThickEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); // Compute numerical result @@ -421,13 +474,13 @@ int testLayerThicknessAuxVars(const Array2DReal &LayerThickCell, ErrorMeasures FluxThickErrors; Err += computeErrors(FluxThickErrors, NumFluxLayerThickEdge, ExactThickEdge, - Mesh, OnEdge, NVertLevels); + Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "FluxThick", FluxThickErrors, Setup.ExpectedFluxThickErrors, RTol); ErrorMeasures MeanThickErrors; Err += computeErrors(MeanThickErrors, NumMeanLayerThickEdge, ExactThickEdge, - Mesh, OnEdge, NVertLevels); + Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "MeanThick", MeanThickErrors, Setup.ExpectedMeanThickErrors, RTol); @@ -453,8 +506,7 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.relativeVorticity(X, Y); }, - ExactRelVortVertex, Geom, Mesh, OnVertex, NVertLevels, - ExchangeHalos::No); + ExactRelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); Array2DReal ExactNormRelVortVertex("ExactNormRelVortVertex", Mesh->NVerticesOwned, NVertLevels); @@ -462,8 +514,7 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normalizedRelativeVorticity(X, Y); }, - ExactNormRelVortVertex, Geom, Mesh, OnVertex, NVertLevels, - ExchangeHalos::No); + ExactNormRelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); Array2DReal ExactNormPlanetVortVertex("ExactNormPlanetVortVertex", Mesh->NVerticesOwned, NVertLevels); @@ -471,8 +522,7 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normalizedPlanetaryVorticity(X, Y); }, - ExactNormPlanetVortVertex, Geom, Mesh, OnVertex, NVertLevels, - ExchangeHalos::No); + ExactNormPlanetVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); // Compute numerical results for vertex variables @@ -491,20 +541,20 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, ErrorMeasures RelVortVertexErrors; Err += computeErrors(RelVortVertexErrors, NumRelVortVertex, - ExactRelVortVertex, Mesh, OnVertex, NVertLevels); + ExactRelVortVertex, Mesh, OnVertex); Err += checkErrors("AuxVarsTest", "RelVortVertex", RelVortVertexErrors, Setup.ExpectedRelVortVertexErrors, RTol); ErrorMeasures NormRelVortVertexErrors; Err += computeErrors(NormRelVortVertexErrors, NumNormRelVortVertex, - ExactNormRelVortVertex, Mesh, OnVertex, NVertLevels); + ExactNormRelVortVertex, Mesh, OnVertex); Err += checkErrors("AuxVarsTest", "NormRelVortVertex", NormRelVortVertexErrors, Setup.ExpectedNormRelVortVertexErrors, RTol); ErrorMeasures NormPlanetVortVertexErrors; Err += computeErrors(NormPlanetVortVertexErrors, NumNormPlanetVortVertex, - ExactNormPlanetVortVertex, Mesh, OnVertex, NVertLevels); + ExactNormPlanetVortVertex, Mesh, OnVertex); Err += checkErrors("AuxVarsTest", "NormPlanetVortVertex", NormPlanetVortVertexErrors, Setup.ExpectedNormPlanetVortVertexErrors, RTol); @@ -517,8 +567,7 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normalizedRelativeVorticity(X, Y); }, - ExactNormRelVortEdge, Geom, Mesh, OnEdge, NVertLevels, - ExchangeHalos::No); + ExactNormRelVortEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); Array2DReal ExactNormPlanetVortEdge("ExactNormPlanetVortEdge", Mesh->NEdgesOwned, NVertLevels); @@ -526,8 +575,7 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normalizedPlanetaryVorticity(X, Y); }, - ExactNormPlanetVortEdge, Geom, Mesh, OnEdge, NVertLevels, - ExchangeHalos::No); + ExactNormPlanetVortEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); // Compute numerical results for vertex variables @@ -542,13 +590,13 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, ErrorMeasures NormRelVortEdgeErrors; Err += computeErrors(NormRelVortEdgeErrors, NumNormRelVortEdge, - ExactNormRelVortEdge, Mesh, OnEdge, NVertLevels); + ExactNormRelVortEdge, Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "NormRelVortEdge", NormRelVortEdgeErrors, Setup.ExpectedNormRelVortEdgeErrors, RTol); ErrorMeasures NormPlanetVortEdgeErrors; Err += computeErrors(NormPlanetVortEdgeErrors, NumNormPlanetVortEdge, - ExactNormPlanetVortEdge, Mesh, OnEdge, NVertLevels); + ExactNormPlanetVortEdge, Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "NormPlanetVortEdge", NormPlanetVortEdgeErrors, Setup.ExpectedNormPlanetVortEdgeErrors, RTol); @@ -574,13 +622,13 @@ int testVelocityDel2AuxVars(Real RTol) { NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, - ExactVelocityDivCell, Geom, Mesh, OnCell, NVertLevels); + ExactVelocityDivCell, Geom, Mesh, OnCell); Array2DReal ExactRelVortVertex("ExactRelVortVertex", Mesh->NVerticesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.relativeVorticity(X, Y); }, - ExactRelVortVertex, Geom, Mesh, OnVertex, NVertLevels); + ExactRelVortVertex, Geom, Mesh, OnVertex); // Compute exact Del2 @@ -590,8 +638,7 @@ int testVelocityDel2AuxVars(Real RTol) { VecField[0] = Setup.velocityDel2X(X, Y); VecField[1] = Setup.velocityDel2Y(X, Y); }, - ExactDel2Edge, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactDel2Edge, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Compute numerical Del2 @@ -606,8 +653,7 @@ int testVelocityDel2AuxVars(Real RTol) { // Compute error measures and check errors for Del2 ErrorMeasures Del2Errors; - Err += computeErrors(Del2Errors, NumDel2Edge, ExactDel2Edge, Mesh, OnEdge, - NVertLevels); + Err += computeErrors(Del2Errors, NumDel2Edge, ExactDel2Edge, Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "Del2", Del2Errors, Setup.ExpectedDel2Errors, RTol); @@ -618,7 +664,7 @@ int testVelocityDel2AuxVars(Real RTol) { NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.velocityDel2Div(X, Y); }, - ExactDel2DivCell, Geom, Mesh, OnCell, NVertLevels, ExchangeHalos::No); + ExactDel2DivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical Del2Div @@ -632,7 +678,7 @@ int testVelocityDel2AuxVars(Real RTol) { ErrorMeasures Del2DivErrors; Err += computeErrors(Del2DivErrors, NumDel2DivCell, ExactDel2DivCell, Mesh, - OnCell, NVertLevels); + OnCell); Err += checkErrors("AuxVarsTest", "Del2Div", Del2DivErrors, Setup.ExpectedDel2DivErrors, RTol); @@ -642,8 +688,7 @@ int testVelocityDel2AuxVars(Real RTol) { Mesh->NVerticesOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.velocityDel2Curl(X, Y); }, - ExactDel2RelVortVertex, Geom, Mesh, OnVertex, NVertLevels, - ExchangeHalos::No); + ExactDel2RelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); // Compute numerical Del2RelVort @@ -658,7 +703,7 @@ int testVelocityDel2AuxVars(Real RTol) { ErrorMeasures Del2RelVortErrors; Err += computeErrors(Del2RelVortErrors, NumDel2RelVortVertex, - ExactDel2RelVortVertex, Mesh, OnVertex, NVertLevels); + ExactDel2RelVortVertex, Mesh, OnVertex); Err += checkErrors("AuxVarsTest", "Del2RelVort", Del2RelVortErrors, Setup.ExpectedDel2RelVortErrors, RTol); @@ -686,12 +731,12 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracersOnCell, Geom, Mesh, OnCell, NVertLevels, NTracers); + TracersOnCell, Geom, Mesh, OnCell); Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickEdge, Geom, Mesh, OnEdge, NVertLevels); + LayerThickEdge, Geom, Mesh, OnEdge); // Compute exact HTracerEdge @@ -699,8 +744,7 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.thickTracer(X, Y); }, - ExactHTrEdge, Geom, Mesh, OnEdge, NVertLevels, NTracers, - ExchangeHalos::No); + ExactHTrEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); // Compute numerical HTracersEdge @@ -716,8 +760,7 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, const auto &NumHTrEdge = TracerAux.HTracersEdge; ErrorMeasures HTracerErrors; - Err += computeErrors(HTracerErrors, NumHTrEdge, ExactHTrEdge, Mesh, OnEdge, - NVertLevels, NTracers); + Err += computeErrors(HTracerErrors, NumHTrEdge, ExactHTrEdge, Mesh, OnEdge); Err += checkErrors("AuxVarsTest", "HTracers", HTracerErrors, Setup.ExpectedHTracerErrors, RTol); @@ -727,8 +770,7 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.del2Tracer(X, Y); }, - ExactDel2TrCell, Geom, Mesh, OnCell, NVertLevels, NTracers, - ExchangeHalos::No); + ExactDel2TrCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical Del2TracerCell @@ -745,7 +787,7 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, ErrorMeasures Del2TracerErrors; Err += computeErrors(Del2TracerErrors, NumDel2TrCell, ExactDel2TrCell, Mesh, - OnCell, NVertLevels, NTracers); + OnCell); Err += checkErrors("AuxVarsTest", "Del2Tracers", Del2TracerErrors, Setup.ExpectedDel2TracerErrors, RTol); @@ -823,6 +865,8 @@ int auxVarsTest(const std::string &mesh = DefaultMeshFile) { Err += testTracerAuxVars(LayerThickCell, NormalVelEdge, RTol); + Err += testWindForcingAuxVars(RTol); + if (Err == 0) { LOG_INFO("AuxVarsTest: Successful completion"); } diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index f70cbc859a92..bf559e43c515 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -29,14 +29,18 @@ struct TestSetupPlane { Real Lx = 1; Real Ly = std::sqrt(3) / 2; - ErrorMeasures ExpectedDivErrors = {0.00124886886594427027, - 0.00124886886590974385}; - ErrorMeasures ExpectedGradErrors = {0.00125026071878537952, - 0.00134354611117262204}; - ErrorMeasures ExpectedCurlErrors = {0.161365663569699946, - 0.161348016897141039}; - ErrorMeasures ExpectedReconErrors = {0.00450897496974901352, - 0.00417367308684470691}; + ErrorMeasures ExpectedDivErrors = {0.00124886886594427027, + 0.00124886886590974385}; + ErrorMeasures ExpectedGradErrors = {0.00125026071878537952, + 0.00134354611117262204}; + ErrorMeasures ExpectedCurlErrors = {0.161365663569699946, + 0.161348016897141039}; + ErrorMeasures ExpectedReconErrors = {0.00450897496974901352, + 0.00417367308684470691}; + ErrorMeasures ExpectedAnisoInterpErrors = {0.0026762081503380526, + 0.003058198461518835}; + ErrorMeasures ExpectedIsoInterpErrors = {0.004279097382993937, + 0.004200067675522098}; KOKKOS_FUNCTION Real exactScalar(Real X, Real Y) const { return std::sin(2 * Pi * X / Lx) * std::sin(2 * Pi * Y / Ly); @@ -76,14 +80,18 @@ struct TestSetupSphere1 { // 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}; - ErrorMeasures ExpectedReconErrors = {0.0206375134079833517, - 0.00692590524910695858}; + ErrorMeasures ExpectedDivErrors = {0.013659577398978353, + 0.00367052484586382743}; + ErrorMeasures ExpectedGradErrors = {0.00187912292540628936, + 0.00149841802817334306}; + ErrorMeasures ExpectedCurlErrors = {0.0271404735181308317, + 0.025202316610921989}; + ErrorMeasures ExpectedReconErrors = {0.0206375134079833517, + 0.00692590524910695858}; + ErrorMeasures ExpectedAnisoInterpErrors = {0.0024015775047603197, + 0.0018490649516209202}; + ErrorMeasures ExpectedIsoInterpErrors = {0.007438367234983312, + 0.0029921955942401697}; KOKKOS_FUNCTION Real exactScalar(Real Lon, Real Lat) const { return Radius * std::cos(Lon) * std::pow(std::cos(Lat), 4); @@ -122,14 +130,18 @@ struct TestSetupSphere2 { // TODO: get this from the mesh Real Radius = 6371220; - ErrorMeasures ExpectedDivErrors = {1.37734693033362766e-10, - 0.000484370621558727582}; - ErrorMeasures ExpectedGradErrors = {0.000906351303388669991, - 0.000949206041390823676}; - ErrorMeasures ExpectedCurlErrors = {0.00433205620592059647, - 0.00204725417666192042}; - ErrorMeasures ExpectedReconErrors = {0.0254271921029878764, - 0.00419630561428921064}; + ErrorMeasures ExpectedDivErrors = {1.37734693033362766e-10, + 0.000484370621558727582}; + ErrorMeasures ExpectedGradErrors = {0.000906351303388669991, + 0.000949206041390823676}; + ErrorMeasures ExpectedCurlErrors = {0.00433205620592059647, + 0.00204725417666192042}; + ErrorMeasures ExpectedReconErrors = {0.0254271921029878764, + 0.00419630561428921064}; + ErrorMeasures ExpectedAnisoInterpErrors = {0.0014465229922953644, + 0.001643777653612931}; + ErrorMeasures ExpectedIsoInterpErrors = {0.004755875091568591, + 0.0025556382734782538}; KOKKOS_FUNCTION Real exactScalar(Real Lon, Real Lat) const { return -Radius * std::pow(std::sin(Lat), 2); @@ -184,13 +196,13 @@ int testDivergence(Real RTol) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, - VecEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + VecEdge, EdgeComponent::Normal, Geom, Mesh); // Compute exact result Array2DReal ExactDivCell("ExactDivCell", Mesh->NCellsOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.exactDivVec(X, Y); }, - ExactDivCell, Geom, Mesh, OnCell, NVertLevels, ExchangeHalos::No); + ExactDivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical result Array2DReal NumDivCell("NumDivCell", Mesh->NCellsOwned, NVertLevels); @@ -202,8 +214,7 @@ int testDivergence(Real RTol) { // Compute error measures ErrorMeasures DivErrors; - Err += computeErrors(DivErrors, NumDivCell, ExactDivCell, Mesh, OnCell, - NVertLevels); + Err += computeErrors(DivErrors, NumDivCell, ExactDivCell, Mesh, OnCell); // Check error values Err += checkErrors("OperatorsTest", "Divergence", DivErrors, Setup.ExpectedDivErrors, RTol); @@ -228,7 +239,7 @@ int testGradient(Real RTol) { KOKKOS_LAMBDA(Real Coord1, Real Coord2) { return Setup.exactScalar(Coord1, Coord2); }, - ScalarCell, Geom, Mesh, OnCell, NVertLevels); + ScalarCell, Geom, Mesh, OnCell); // Compute exact result Array2DReal ExactGradEdge("ExactGradEdge", Mesh->NEdgesOwned, NVertLevels); @@ -237,8 +248,7 @@ int testGradient(Real RTol) { VecField[0] = Setup.exactGradScalarX(X, Y); VecField[1] = Setup.exactGradScalarY(X, Y); }, - ExactGradEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactGradEdge, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Compute numerical result GradientOnEdge GradientEdge(Mesh); @@ -250,8 +260,7 @@ int testGradient(Real RTol) { // Compute error measures ErrorMeasures GradErrors; - Err += computeErrors(GradErrors, NumGradEdge, ExactGradEdge, Mesh, OnEdge, - NVertLevels); + Err += computeErrors(GradErrors, NumGradEdge, ExactGradEdge, Mesh, OnEdge); // Check error values Err += checkErrors("OperatorsTest", "Gradient", GradErrors, Setup.ExpectedGradErrors, RTol); @@ -276,14 +285,14 @@ int testCurl(Real RTol) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, - VecEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + VecEdge, EdgeComponent::Normal, Geom, Mesh); // Compute exact result Array2DReal ExactCurlVertex("ExactCurlVertex", Mesh->NVerticesOwned, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.exactCurlVec(X, Y); }, - ExactCurlVertex, Geom, Mesh, OnVertex, NVertLevels, ExchangeHalos::No); + ExactCurlVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); // Compute numerical result Array2DReal NumCurlVertex("NumCurlVertex", Mesh->NVerticesOwned, @@ -297,7 +306,7 @@ int testCurl(Real RTol) { // Compute error measures ErrorMeasures CurlErrors; Err += computeErrors(CurlErrors, NumCurlVertex, ExactCurlVertex, Mesh, - OnVertex, NVertLevels); + OnVertex); // Check error values Err += checkErrors("OperatorsTest", "Curl", CurlErrors, Setup.ExpectedCurlErrors, RTol); @@ -323,7 +332,7 @@ int testRecon(Real RTol) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, - VecEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + VecEdge, EdgeComponent::Normal, Geom, Mesh); // Compute exact result Array2DReal ExactReconEdge("ExactReconEdge", Mesh->NEdgesOwned, NVertLevels); @@ -333,7 +342,7 @@ int testRecon(Real RTol) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, - ExactReconEdge, EdgeComponent::Tangential, Geom, Mesh, NVertLevels, + ExactReconEdge, EdgeComponent::Tangential, Geom, Mesh, ExchangeHalos::No); // Compute numerical result @@ -346,8 +355,8 @@ int testRecon(Real RTol) { // Compute error measures ErrorMeasures ReconErrors; - Err += computeErrors(ReconErrors, NumReconEdge, ExactReconEdge, Mesh, OnEdge, - NVertLevels); + Err += + computeErrors(ReconErrors, NumReconEdge, ExactReconEdge, Mesh, OnEdge); // Check error values Err += checkErrors("OperatorsTest", "Recon", ReconErrors, Setup.ExpectedReconErrors, RTol); @@ -359,6 +368,62 @@ int testRecon(Real RTol) { return Err; } +int testInterpCellToEdge(Real RTol) { + int Err = 0; + TestSetup Setup; + const auto &Mesh = HorzMesh::getDefault(); + + // Prepare operator input + Array1DReal ScalarCell("ScalarCell", Mesh->NCellsSize); + Err += setScalar( + KOKKOS_LAMBDA(Real Coord1, Real Coord2) { + return Setup.exactScalar(Coord1, Coord2); + }, + ScalarCell, Geom, Mesh, OnCell); + + // Compute exact result + Array1DReal ExactScalarEdge("ExactScalarEdge", Mesh->NEdgesOwned); + Err += setScalar( + KOKKOS_LAMBDA(Real Coord1, Real Coord2) { + return Setup.exactScalar(Coord1, Coord2); + }, + ExactScalarEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); + + // Compute numerical result + Array1DReal IsoNumScalarEdge("IsoNumScalarEdge", Mesh->NEdgesOwned); + Array1DReal AnisoNumScalarEdge("AnisoNumScalarEdge", Mesh->NEdgesOwned); + InterpCellToEdge Interp(Mesh); + parallelFor( + {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { + AnisoNumScalarEdge(IEdge) = + Interp(IEdge, ScalarCell, InterpCellToEdgeOption::Anisotropic); + IsoNumScalarEdge(IEdge) = + Interp(IEdge, ScalarCell, InterpCellToEdgeOption::Isotropic); + }); + + // Compute error measures + ErrorMeasures AnisoInterpErrors; + Err += computeErrors(AnisoInterpErrors, AnisoNumScalarEdge, ExactScalarEdge, + Mesh, OnEdge); + + ErrorMeasures IsoInterpErrors; + Err += computeErrors(IsoInterpErrors, IsoNumScalarEdge, ExactScalarEdge, + Mesh, OnEdge); + + // Check error values + Err += checkErrors("OperatorsTest", "AnisoInterpCellToEdge", + AnisoInterpErrors, Setup.ExpectedAnisoInterpErrors, RTol); + + Err += checkErrors("OperatorsTest", "IsoInterpCellToEdge", IsoInterpErrors, + Setup.ExpectedIsoInterpErrors, RTol); + + if (Err == 0) { + LOG_INFO("OperatorsTest: InterpCellToEdge PASS"); + } + + return Err; +} + //------------------------------------------------------------------------------ // The initialization routine for Operators testing int initOperatorsTest(const std::string &MeshFile) { @@ -414,6 +479,7 @@ int operatorsTest(const std::string &MeshFile = DefaultMeshFile) { Err += testGradient(RTol); Err += testCurl(RTol); Err += testRecon(RTol); + Err += testInterpCellToEdge(RTol); if (Err == 0) { LOG_INFO("OperatorsTest: Successful completion"); diff --git a/components/omega/test/ocn/OceanTestCommon.h b/components/omega/test/ocn/OceanTestCommon.h index 6e0fc05a927f..0e095dc93283 100644 --- a/components/omega/test/ocn/OceanTestCommon.h +++ b/components/omega/test/ocn/OceanTestCommon.h @@ -11,8 +11,9 @@ namespace OMEGA { // check if two real numbers are equal with a given relative tolerance -inline bool isApprox(Real X, Real Y, Real RTol) { - return std::abs(X - Y) <= RTol * std::max(std::abs(X), std::abs(Y)); +inline bool isApprox(Real X, Real Y, Real RTol, Real ATol = 0) { + return std::abs(X - Y) <= + std::max(ATol, RTol * std::max(std::abs(X), std::abs(Y))); } // convert spherical components of a vector to Cartesian @@ -62,10 +63,9 @@ enum class ExchangeHalos { Yes, No }; // set scalar field on chosen elements (cells/vertices/edges) based on // analytical formula and optionally exchange halos -template -int setScalar(const Functor &Fun, const Array2DReal &ScalarElement, - Geometry Geom, const HorzMesh *Mesh, MeshElement Element, - int NVertLevels, +template +int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, + const HorzMesh *Mesh, MeshElement Element, ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes) { int Err = 0; @@ -102,83 +102,55 @@ int setScalar(const Functor &Fun, const Array2DReal &ScalarElement, return 1; } - parallelFor( - {NElementsOwned, NVertLevels}, KOKKOS_LAMBDA(int IElement, int K) { - if (Geom == Geometry::Planar) { - const Real X = XElement(IElement); - const Real Y = YElement(IElement); - ScalarElement(IElement, K) = Fun(X, Y); - } else { - const Real Lon = LonElement(IElement); - const Real Lat = LatElement(IElement); - ScalarElement(IElement, K) = Fun(Lon, Lat); - } - }); - - if (ExchangeHalosOpt == ExchangeHalos::Yes) { - auto MyHalo = Halo::getDefault(); - Err = MyHalo->exchangeFullArrayHalo(ScalarElement, Element); - if (Err != 0) - LOG_ERROR("setScalarElement: error in halo exchange"); + if constexpr (Array::rank == 1) { + parallelFor( + {NElementsOwned}, KOKKOS_LAMBDA(int IElement) { + if (Geom == Geometry::Planar) { + const Real X = XElement(IElement); + const Real Y = YElement(IElement); + ScalarElement(IElement) = Fun(X, Y); + } else { + const Real Lon = LonElement(IElement); + const Real Lat = LatElement(IElement); + ScalarElement(IElement) = Fun(Lon, Lat); + } + }); } - return Err; -} -// set scalar field on chosen elements (cells/vertices/edges) based on -// analytical formula and optionally exchange halos -template -int setScalar(const Functor &Fun, const Array3DReal &ScalarElement, - Geometry Geom, const HorzMesh *Mesh, MeshElement Element, - int NVertLevels, int NTracers, - ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes) { - - int Err = 0; - - int NElementsOwned; - Array1DReal XElement, YElement; - Array1DReal LonElement, LatElement; + if constexpr (Array::rank == 2) { + const int NVertLevels = ScalarElement.extent_int(1); - switch (Element) { - case OnCell: - NElementsOwned = Mesh->NCellsOwned; - XElement = createDeviceMirrorCopy(Mesh->XCellH); - YElement = createDeviceMirrorCopy(Mesh->YCellH); - LonElement = createDeviceMirrorCopy(Mesh->LonCellH); - LatElement = createDeviceMirrorCopy(Mesh->LatCellH); - break; - case OnVertex: - NElementsOwned = Mesh->NVerticesOwned; - XElement = createDeviceMirrorCopy(Mesh->XVertexH); - YElement = createDeviceMirrorCopy(Mesh->YVertexH); - LonElement = createDeviceMirrorCopy(Mesh->LonVertexH); - LatElement = createDeviceMirrorCopy(Mesh->LatVertexH); - break; - case OnEdge: - NElementsOwned = Mesh->NEdgesOwned; - XElement = createDeviceMirrorCopy(Mesh->XEdgeH); - YElement = createDeviceMirrorCopy(Mesh->YEdgeH); - LonElement = createDeviceMirrorCopy(Mesh->LonEdgeH); - LatElement = createDeviceMirrorCopy(Mesh->LatEdgeH); - break; - default: - LOG_ERROR("setScalar: element needs to be one of (OnCell, OnVertex, " - "OnEdge)"); - return 1; + parallelFor( + {NElementsOwned, NVertLevels}, KOKKOS_LAMBDA(int IElement, int K) { + if (Geom == Geometry::Planar) { + const Real X = XElement(IElement); + const Real Y = YElement(IElement); + ScalarElement(IElement, K) = Fun(X, Y); + } else { + const Real Lon = LonElement(IElement); + const Real Lat = LatElement(IElement); + ScalarElement(IElement, K) = Fun(Lon, Lat); + } + }); } - parallelFor( - {NTracers, NElementsOwned, NVertLevels}, - KOKKOS_LAMBDA(int L, int IElement, int K) { - if (Geom == Geometry::Planar) { - const Real X = XElement(IElement); - const Real Y = YElement(IElement); - ScalarElement(L, IElement, K) = Fun(X, Y); - } else { - const Real Lon = LonElement(IElement); - const Real Lat = LatElement(IElement); - ScalarElement(L, IElement, K) = Fun(Lon, Lat); - } - }); + if constexpr (Array::rank == 3) { + const int NTracers = ScalarElement.extent_int(0); + const int NVertLevels = ScalarElement.extent_int(2); + parallelFor( + {NTracers, NElementsOwned, NVertLevels}, + KOKKOS_LAMBDA(int L, int IElement, int K) { + if (Geom == Geometry::Planar) { + const Real X = XElement(IElement); + const Real Y = YElement(IElement); + ScalarElement(L, IElement, K) = Fun(X, Y); + } else { + const Real Lon = LonElement(IElement); + const Real Lat = LatElement(IElement); + ScalarElement(L, IElement, K) = Fun(Lon, Lat); + } + }); + } if (ExchangeHalosOpt == ExchangeHalos::Yes) { auto MyHalo = Halo::getDefault(); @@ -193,10 +165,9 @@ enum class CartProjection { Yes, No }; // set vector field on edges based on analytical formula and optionally // exchange halos -template -int setVectorEdge(const Functor &Fun, const Array2DReal &VectorFieldEdge, +template +int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, EdgeComponent EdgeComp, Geometry Geom, const HorzMesh *Mesh, - int NVertLevels, ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, CartProjection CartProjectionOpt = CartProjection::Yes) { @@ -221,84 +192,98 @@ int setVectorEdge(const Functor &Fun, const Array2DReal &VectorFieldEdge, auto &CellsOnEdge = Mesh->CellsOnEdge; auto &VerticesOnEdge = Mesh->VerticesOnEdge; - parallelFor( - {Mesh->NEdgesOwned, NVertLevels}, KOKKOS_LAMBDA(int IEdge, int K) { - Real VecFieldEdge; - if (Geom == Geometry::Planar) { - const Real XE = XEdge(IEdge); - const Real YE = YEdge(IEdge); - - Real VecField[2]; - Fun(VecField, XE, YE); - - if (EdgeComp == EdgeComponent::Normal) { - const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); - const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); - VecFieldEdge = - EdgeNormalX * VecField[0] + EdgeNormalY * VecField[1]; - } + auto ProjectVector = KOKKOS_LAMBDA(int IEdge) { + Real VecFieldEdge; + if (Geom == Geometry::Planar) { + const Real XE = XEdge(IEdge); + const Real YE = YEdge(IEdge); + + Real VecField[2]; + Fun(VecField, XE, YE); + + if (EdgeComp == EdgeComponent::Normal) { + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + VecFieldEdge = + EdgeNormalX * VecField[0] + EdgeNormalY * VecField[1]; + } + + if (EdgeComp == EdgeComponent::Tangential) { + const Real EdgeTangentX = -std::sin(AngleEdge(IEdge)); + const Real EdgeTangentY = std::cos(AngleEdge(IEdge)); + VecFieldEdge = + EdgeTangentX * VecField[0] + EdgeTangentY * VecField[1]; + } + } else { + const Real LonE = LonEdge(IEdge); + const Real LatE = LatEdge(IEdge); + + Real VecField[2]; + Fun(VecField, LonE, LatE); + + if (CartProjectionOpt == CartProjection::Yes) { + Real VecFieldCart[3]; + sphereToCartVec(VecFieldCart, VecField, LonE, LatE); + + const Real EdgeCoords[3] = {XEdge[IEdge], YEdge[IEdge], + ZEdge[IEdge]}; + + if (EdgeComp == EdgeComponent::Normal) { + const int JCell1 = CellsOnEdge(IEdge, 1); + const Real CellCoords[3] = {XCell(JCell1), YCell(JCell1), + ZCell(JCell1)}; + + Real EdgeNormal[3]; + tangentVector(EdgeNormal, EdgeCoords, CellCoords); + VecFieldEdge = EdgeNormal[0] * VecFieldCart[0] + + EdgeNormal[1] * VecFieldCart[1] + + EdgeNormal[2] * VecFieldCart[2]; + } + + if (EdgeComp == EdgeComponent::Tangential) { + const int JVertex1 = VerticesOnEdge(IEdge, 1); + const Real VertexCoords[3] = { + XVertex(JVertex1), YVertex(JVertex1), ZVertex(JVertex1)}; + + Real EdgeTangent[3]; + tangentVector(EdgeTangent, EdgeCoords, VertexCoords); + VecFieldEdge = EdgeTangent[0] * VecFieldCart[0] + + EdgeTangent[1] * VecFieldCart[1] + + EdgeTangent[2] * VecFieldCart[2]; + } + } else { + if (EdgeComp == EdgeComponent::Normal) { + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + VecFieldEdge = + EdgeNormalX * VecField[0] + EdgeNormalY * VecField[1]; + } + + if (EdgeComp == EdgeComponent::Tangential) { + const Real EdgeTangentX = -std::sin(AngleEdge(IEdge)); + const Real EdgeTangentY = std::cos(AngleEdge(IEdge)); + VecFieldEdge = + EdgeTangentX * VecField[0] + EdgeTangentY * VecField[1]; + } + } + } + return VecFieldEdge; + }; + + if constexpr (Array::rank == 1) { + parallelFor( + {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { + VectorFieldEdge(IEdge) = ProjectVector(IEdge); + }); + } - if (EdgeComp == EdgeComponent::Tangential) { - const Real EdgeTangentX = -std::sin(AngleEdge(IEdge)); - const Real EdgeTangentY = std::cos(AngleEdge(IEdge)); - VecFieldEdge = - EdgeTangentX * VecField[0] + EdgeTangentY * VecField[1]; - } - } else { - const Real LonE = LonEdge(IEdge); - const Real LatE = LatEdge(IEdge); - - Real VecField[2]; - Fun(VecField, LonE, LatE); - - if (CartProjectionOpt == CartProjection::Yes) { - Real VecFieldCart[3]; - sphereToCartVec(VecFieldCart, VecField, LonE, LatE); - - const Real EdgeCoords[3] = {XEdge[IEdge], YEdge[IEdge], - ZEdge[IEdge]}; - - if (EdgeComp == EdgeComponent::Normal) { - const int JCell1 = CellsOnEdge(IEdge, 1); - const Real CellCoords[3] = {XCell(JCell1), YCell(JCell1), - ZCell(JCell1)}; - - Real EdgeNormal[3]; - tangentVector(EdgeNormal, EdgeCoords, CellCoords); - VecFieldEdge = EdgeNormal[0] * VecFieldCart[0] + - EdgeNormal[1] * VecFieldCart[1] + - EdgeNormal[2] * VecFieldCart[2]; - } - - if (EdgeComp == EdgeComponent::Tangential) { - const int JVertex1 = VerticesOnEdge(IEdge, 1); - const Real VertexCoords[3] = { - XVertex(JVertex1), YVertex(JVertex1), ZVertex(JVertex1)}; - - Real EdgeTangent[3]; - tangentVector(EdgeTangent, EdgeCoords, VertexCoords); - VecFieldEdge = EdgeTangent[0] * VecFieldCart[0] + - EdgeTangent[1] * VecFieldCart[1] + - EdgeTangent[2] * VecFieldCart[2]; - } - } else { - if (EdgeComp == EdgeComponent::Normal) { - const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); - const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); - VecFieldEdge = - EdgeNormalX * VecField[0] + EdgeNormalY * VecField[1]; - } - - if (EdgeComp == EdgeComponent::Tangential) { - const Real EdgeTangentX = -std::sin(AngleEdge(IEdge)); - const Real EdgeTangentY = std::cos(AngleEdge(IEdge)); - VecFieldEdge = - EdgeTangentX * VecField[0] + EdgeTangentY * VecField[1]; - } - } - } - VectorFieldEdge(IEdge, K) = VecFieldEdge; - }); + if constexpr (Array::rank == 2) { + const int NVertLevels = VectorFieldEdge.extent_int(1); + parallelFor( + {Mesh->NEdgesOwned, NVertLevels}, KOKKOS_LAMBDA(int IEdge, int K) { + VectorFieldEdge(IEdge, K) = ProjectVector(IEdge); + }); + } if (ExchangeHalosOpt == ExchangeHalos::Yes) { auto MyHalo = Halo::getDefault(); @@ -309,6 +294,19 @@ int setVectorEdge(const Functor &Fun, const Array2DReal &VectorFieldEdge, return Err; } +inline Real maxVal(const Array1DReal &Arr) { + Real MaxVal; + + parallelReduce( + {Arr.extent_int(0)}, + KOKKOS_LAMBDA(int I, Real &Accum) { + Accum = Kokkos::max(Arr(I), Accum); + }, + Kokkos::Max(MaxVal)); + + return MaxVal; +} + inline Real maxVal(const Array2DReal &Arr) { Real MaxVal; @@ -335,6 +333,17 @@ inline Real maxVal(const Array3DReal &Arr) { return MaxVal; } +inline Real sum(const Array1DReal &Arr, int Extent0) { + Real Sum; + + parallelReduce( + {Extent0}, KOKKOS_LAMBDA(int I, Real &Accum) { Accum += Arr(I); }, Sum); + + return Sum; +} + +inline Real sum(const Array1DReal &Arr) { return sum(Arr, Arr.extent_int(0)); } + inline Real sum(const Array2DReal &Arr, int Extent0, int Extent1) { Real Sum; @@ -385,11 +394,10 @@ struct ErrorMeasures { // compute global normalized error measures based on the difference // between two fields -inline int computeErrors(ErrorMeasures &ErrorMeasures, - const Array2DReal &NumFieldElement, - const Array2DReal &ExactFieldElement, - const HorzMesh *Mesh, MeshElement Element, - int NVertLevels) { +template +int computeErrors(ErrorMeasures &ErrorMeasures, const Array &NumFieldElement, + const Array &ExactFieldElement, const HorzMesh *Mesh, + MeshElement Element) { int Err = 0; @@ -425,135 +433,93 @@ inline int computeErrors(ErrorMeasures &ErrorMeasures, return 1; } - // Compute element-wise errors - Array2DReal LInfElement("LInfElement", NElementsOwned, NVertLevels); - Array2DReal L2Element("L2Element", NElementsOwned, NVertLevels); - - Array2DReal LInfScaleElement("LInfScaleElement", NElementsOwned, - NVertLevels); - Array2DReal L2ScaleElement("L2ScaleElement", NElementsOwned, NVertLevels); - - parallelFor( - {NElementsOwned, NVertLevels}, KOKKOS_LAMBDA(int IElement, int K) { - const Real NumValElement = NumFieldElement(IElement, K); - const Real ExactValElement = ExactFieldElement(IElement, K); - - // Errors - LInfElement(IElement, K) = std::abs(NumValElement - ExactValElement); - LInfScaleElement(IElement, K) = std::abs(ExactValElement); - L2Element(IElement, K) = AreaElement(IElement) * - LInfElement(IElement, K) * - LInfElement(IElement, K); - L2ScaleElement(IElement, K) = AreaElement(IElement) * - LInfScaleElement(IElement, K) * - LInfScaleElement(IElement, K); - }); + Array LInfElement; + Array L2Element; + Array LInfScaleElement; + Array L2ScaleElement; - // Compute global normalized error norms - const Real LInfErrorLoc = maxVal(LInfElement); - const Real L2ErrorLoc = sum(L2Element); - const Real LInfScaleLoc = maxVal(LInfScaleElement); - const Real L2ScaleLoc = sum(L2ScaleElement); - - MPI_Comm Comm = MachEnv::getDefault()->getComm(); - Real LInfError, LInfScale; - Err += - MPI_Allreduce(&LInfErrorLoc, &LInfError, 1, MPI_RealKind, MPI_MAX, Comm); - Err += - MPI_Allreduce(&LInfScaleLoc, &LInfScale, 1, MPI_RealKind, MPI_MAX, Comm); - if (LInfScale > 0) { - LInfError /= LInfScale; + // Compute element-wise errors + if constexpr (Array::rank == 1) { + LInfElement = Array("LInfElement", NElementsOwned); + L2Element = Array("L2Element", NElementsOwned); + + LInfScaleElement = Array("LInfScaleElement", NElementsOwned); + L2ScaleElement = Array("L2ScaleElement", NElementsOwned); + + parallelFor( + {NElementsOwned}, KOKKOS_LAMBDA(int IElement) { + const Real NumValElement = NumFieldElement(IElement); + const Real ExactValElement = ExactFieldElement(IElement); + + // Errors + LInfElement(IElement) = std::abs(NumValElement - ExactValElement); + LInfScaleElement(IElement) = std::abs(ExactValElement); + L2Element(IElement) = AreaElement(IElement) * + LInfElement(IElement) * + LInfElement(IElement); + L2ScaleElement(IElement) = AreaElement(IElement) * + LInfScaleElement(IElement) * + LInfScaleElement(IElement); + }); } - - Real L2Error, L2Scale; - Err += MPI_Allreduce(&L2ErrorLoc, &L2Error, 1, MPI_RealKind, MPI_SUM, Comm); - Err += MPI_Allreduce(&L2ScaleLoc, &L2Scale, 1, MPI_RealKind, MPI_SUM, Comm); - - if (Err != 0) - LOG_ERROR("computeErrors: MPI Allreduce error"); - - if (L2Scale > 0) { - L2Error = std::sqrt(L2Error / L2Scale); - } else { - L2Error = std::sqrt(L2Error); + if constexpr (Array::rank == 2) { + const int NVertLevels = NumFieldElement.extent_int(1); + + LInfElement = Array("LInfElement", NElementsOwned, NVertLevels); + L2Element = Array("L2Element", NElementsOwned, NVertLevels); + + LInfScaleElement = Array("LInfScaleElement", NElementsOwned, NVertLevels); + L2ScaleElement = Array("L2ScaleElement", NElementsOwned, NVertLevels); + + parallelFor( + {NElementsOwned, NVertLevels}, KOKKOS_LAMBDA(int IElement, int K) { + const Real NumValElement = NumFieldElement(IElement, K); + const Real ExactValElement = ExactFieldElement(IElement, K); + + // Errors + LInfElement(IElement, K) = + std::abs(NumValElement - ExactValElement); + LInfScaleElement(IElement, K) = std::abs(ExactValElement); + L2Element(IElement, K) = AreaElement(IElement) * + LInfElement(IElement, K) * + LInfElement(IElement, K); + L2ScaleElement(IElement, K) = AreaElement(IElement) * + LInfScaleElement(IElement, K) * + LInfScaleElement(IElement, K); + }); } - ErrorMeasures.L2 = L2Error; - ErrorMeasures.LInf = LInfError; - - return Err; -} - -// compute global normalized error measures based on the difference -// between two fields -inline int computeErrors(ErrorMeasures &ErrorMeasures, - const Array3DReal &NumFieldElement, - const Array3DReal &ExactFieldElement, - const HorzMesh *Mesh, MeshElement Element, - int NVertLevels, int NTracers) { - - int Err = 0; - - int NElementsOwned; - Array1DReal AreaElement; - - switch (Element) { - case OnCell: - NElementsOwned = Mesh->NCellsOwned; - AreaElement = Mesh->AreaCell; - break; - case OnVertex: - NElementsOwned = Mesh->NVerticesOwned; - AreaElement = Mesh->AreaTriangle; - break; - case OnEdge: - NElementsOwned = Mesh->NEdgesOwned; - // need to compute areas associated with edges since we don't store those - // in the mesh class - { - auto &DcEdge = Mesh->DcEdge; - auto &DvEdge = Mesh->DvEdge; - AreaElement = Array1DReal("AreaEdge", Mesh->NEdgesOwned); - parallelFor( - {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { - AreaElement(IEdge) = DcEdge(IEdge) * DvEdge(IEdge) / 2; - }); - } - break; - default: - LOG_ERROR("computeErrors: element needs to be one of (OnCell, OnVertex, " - "OnEdge)"); - return 1; + if constexpr (Array::rank == 3) { + const int NTracers = NumFieldElement.extent_int(0); + const int NVertLevels = NumFieldElement.extent_int(2); + + LInfElement = Array("LInfElement", NTracers, NElementsOwned, NVertLevels); + L2Element = Array("L2Element", NTracers, NElementsOwned, NVertLevels); + + LInfScaleElement = + Array("LInfScaleElement", NTracers, NElementsOwned, NVertLevels); + L2ScaleElement = + Array("L2ScaleElement", NTracers, NElementsOwned, NVertLevels); + + parallelFor( + {NTracers, NElementsOwned, NVertLevels}, + KOKKOS_LAMBDA(int L, int IElement, int K) { + const Real NumValElement = NumFieldElement(L, IElement, K); + const Real ExactValElement = ExactFieldElement(L, IElement, K); + + // Errors + LInfElement(L, IElement, K) = + std::abs(NumValElement - ExactValElement); + LInfScaleElement(L, IElement, K) = std::abs(ExactValElement); + L2Element(L, IElement, K) = AreaElement(IElement) * + LInfElement(L, IElement, K) * + LInfElement(L, IElement, K); + L2ScaleElement(L, IElement, K) = AreaElement(IElement) * + LInfScaleElement(L, IElement, K) * + LInfScaleElement(L, IElement, K); + }); } - // Compute element-wise errors - Array3DReal LInfElement("LInfElement", NTracers, NElementsOwned, - NVertLevels); - Array3DReal L2Element("L2Element", NTracers, NElementsOwned, NVertLevels); - - Array3DReal LInfScaleElement("LInfScaleElement", NTracers, NElementsOwned, - NVertLevels); - Array3DReal L2ScaleElement("L2ScaleElement", NTracers, NElementsOwned, - NVertLevels); - - parallelFor( - {NTracers, NElementsOwned, NVertLevels}, - KOKKOS_LAMBDA(int L, int IElement, int K) { - const Real NumValElement = NumFieldElement(L, IElement, K); - const Real ExactValElement = ExactFieldElement(L, IElement, K); - - // Errors - LInfElement(L, IElement, K) = - std::abs(NumValElement - ExactValElement); - LInfScaleElement(L, IElement, K) = std::abs(ExactValElement); - L2Element(L, IElement, K) = AreaElement(IElement) * - LInfElement(L, IElement, K) * - LInfElement(L, IElement, K); - L2ScaleElement(L, IElement, K) = AreaElement(IElement) * - LInfScaleElement(L, IElement, K) * - LInfScaleElement(L, IElement, K); - }); - // Compute global normalized error norms const Real LInfErrorLoc = maxVal(LInfElement); const Real L2ErrorLoc = sum(L2Element); @@ -591,14 +557,15 @@ inline int computeErrors(ErrorMeasures &ErrorMeasures, inline int checkErrors(const std::string &TestSuite, const std::string &Variable, const ErrorMeasures &Errors, - const ErrorMeasures &ExpectedErrors, Real RTol) { + const ErrorMeasures &ExpectedErrors, Real RTol, + Real ATol = 0) { int Err = 0; - if (!isApprox(Errors.LInf, ExpectedErrors.LInf, RTol)) { + if (!isApprox(Errors.LInf, ExpectedErrors.LInf, RTol, ATol)) { Err++; LOG_ERROR("{}: {} LInf FAIL, expected {}, got {}", TestSuite, Variable, ExpectedErrors.LInf, Errors.LInf); } - if (!isApprox(Errors.L2, ExpectedErrors.L2, RTol)) { + if (!isApprox(Errors.L2, ExpectedErrors.L2, RTol, ATol)) { Err++; LOG_ERROR("{}: {} L2 FAIL, expected {}, got {}", TestSuite, Variable, ExpectedErrors.L2, Errors.L2); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index fab00781c36d..5c60741c50a0 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -66,19 +66,19 @@ int initState() { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickCell, Geom, Mesh, OnCell, NVertLevels); + LayerThickCell, Geom, Mesh, OnCell); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracersCell, Geom, Mesh, OnCell, NVertLevels, NTracers); + TracersCell, Geom, Mesh, OnCell); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { VecField[0] = Setup.velocityX(Lon, Lat); VecField[1] = Setup.velocityY(Lon, Lat); }, - NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::Yes, CartProjection::No); + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::Yes, + CartProjection::No); return Err; } diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index be30b4610bed..29699e5c73af 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -29,6 +29,7 @@ #include "mpi.h" #include +#include using namespace OMEGA; @@ -38,20 +39,23 @@ struct TestSetupPlane { Real Lx = 1; Real Ly = std::sqrt(3) / 2; - ErrorMeasures ExpectedDivErrors = {0.00124886886594453264, - 0.00124886886590977139}; - ErrorMeasures ExpectedPVErrors = {0.00807347170900282914, - 0.00794755105765788429}; - ErrorMeasures ExpectedGradErrors = {0.00125026071878537952, - 0.00134354611117262161}; - ErrorMeasures ExpectedLaplaceErrors = {0.00113090174765822192, - 0.00134324628763667899}; - ErrorMeasures ExpectedTrHAdvErrors = {0.00205864372747571571, - 0.00172418025417940784}; - ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, - 0.00290978146207349032}; - ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, - 0.00523080740758275625}; + ErrorMeasures ExpectedDivErrors = {0.00124886886594453264, + 0.00124886886590977139}; + ErrorMeasures ExpectedPVErrors = {0.00807347170900282914, + 0.00794755105765788429}; + ErrorMeasures ExpectedGradErrors = {0.00125026071878537952, + 0.00134354611117262161}; + ErrorMeasures ExpectedLaplaceErrors = {0.00113090174765822192, + 0.00134324628763667899}; + ErrorMeasures ExpectedTrHAdvErrors = {0.00205864372747571571, + 0.00172418025417940784}; + ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, + 0.00290978146207349032}; + ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, + 0.00523080740758275625}; + ErrorMeasures ExpectedWindForcingErrors = {0, 0}; + ErrorMeasures ExpectedBottomDragErrors = {0.033848740052302935, + 0.01000133508329411}; KOKKOS_FUNCTION Real vectorX(Real X, Real Y) const { return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); @@ -147,6 +151,30 @@ struct TestSetupPlane { std::cos(4 * Pi * Y / Ly) / Ly / Ly); } + KOKKOS_FUNCTION Real windForcingX(Real X, Real Y, + Real SaltWaterDensity) const { + const Real StressU = vectorX(X, Y); + const Real Thick = scalarB(X, Y); + return StressU / (Thick * SaltWaterDensity); + } + + KOKKOS_FUNCTION Real windForcingY(Real X, Real Y, + Real SaltWaterDensity) const { + const Real StressV = vectorY(X, Y); + const Real Thick = scalarB(X, Y); + return StressV / (Thick * SaltWaterDensity); + } + + KOKKOS_FUNCTION Real bottomDragX(Real X, Real Y, Real Coeff) const { + const Real UVel = vectorX(X, Y); + return -Coeff * std::abs(scalarA(X, Y)) / scalarB(X, Y) * UVel; + } + + KOKKOS_FUNCTION Real bottomDragY(Real X, Real Y, Real Coeff) const { + const Real VVel = vectorY(X, Y); + return -Coeff * std::abs(scalarA(X, Y)) / scalarB(X, Y) * VVel; + } + }; // end TestSetupPlane struct TestSetupSphere { @@ -156,20 +184,23 @@ struct TestSetupSphere { Real Pi = M_PI; - ErrorMeasures ExpectedDivErrors = {0.0136595773989796766, - 0.00367052484586384131}; - 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 ExpectedDivErrors = {0.0136595773989796766, + 0.00367052484586384131}; + 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 ExpectedWindForcingErrors = {0, 0}; + ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, + 0.0014897009917655022}; KOKKOS_FUNCTION Real vectorX(Real Lon, Real Lat) const { return -Radius * std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); @@ -265,6 +296,30 @@ struct TestSetupSphere { return std::sqrt(3 / 2 / Pi) * std::cos(Lat) * std::cos(Lon) / Radius; } + KOKKOS_FUNCTION Real windForcingX(Real Lon, Real Lat, + Real SaltWaterDensity) const { + const Real StressU = vectorX(Lon, Lat); + const Real Thick = scalarB(Lon, Lat); + return StressU / (Thick * SaltWaterDensity); + } + + KOKKOS_FUNCTION Real windForcingY(Real Lon, Real Lat, + Real SaltWaterDensity) const { + const Real StressV = vectorY(Lon, Lat); + const Real Thick = scalarB(Lon, Lat); + return StressV / (Thick * SaltWaterDensity); + } + + KOKKOS_FUNCTION Real bottomDragX(Real Lon, Real Lat, Real Coeff) const { + const Real UVel = vectorX(Lon, Lat); + return -Coeff * std::abs(scalarA(Lon, Lat)) / scalarB(Lon, Lat) * UVel; + } + + KOKKOS_FUNCTION Real bottomDragY(Real Lon, Real Lat, Real Coeff) const { + const Real VVel = vectorY(Lon, Lat); + return -Coeff * std::abs(scalarA(Lon, Lat)) / scalarB(Lon, Lat) * VVel; + } + }; // end TestSetupSphere #ifdef TENDENCYTERMS_TEST_PLANE @@ -290,7 +345,7 @@ int testThickFluxDiv(int NVertLevels, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.divergence(X, Y); }, - ExactThickFluxDiv, Geom, Mesh, OnCell, NVertLevels, ExchangeHalos::No); + ExactThickFluxDiv, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input array Array2DReal ThickFluxEdge("ThickFluxEdge", Mesh->NEdgesSize, NVertLevels); @@ -304,7 +359,7 @@ int testThickFluxDiv(int NVertLevels, Real RTol) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, - ThickFluxEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + ThickFluxEdge, EdgeComponent::Normal, Geom, Mesh); // Compute numerical result Array2DReal NumThickFluxDiv("NumThickFluxDiv", Mesh->NCellsOwned, @@ -319,7 +374,7 @@ int testThickFluxDiv(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures TFDivErrors; Err += computeErrors(TFDivErrors, NumThickFluxDiv, ExactThickFluxDiv, Mesh, - OnCell, NVertLevels); + OnCell); // Check error values Err += checkErrors("TendencyTermsTest", "ThickFluxDiv", TFDivErrors, @@ -350,7 +405,7 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { VecField[1] = (Setup.normRelVort(X, Y) + Setup.normPlanetVort(X, Y)) * Setup.layerThick(X, Y) * Setup.vectorY(X, Y); }, - ExactPotVortHAdv, EdgeComponent::Tangential, Geom, Mesh, NVertLevels, + ExactPotVortHAdv, EdgeComponent::Tangential, Geom, Mesh, ExchangeHalos::No); // Set input arrays @@ -359,19 +414,19 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normRelVort(X, Y); }, - NormRelVortEdge, Geom, Mesh, OnEdge, NVertLevels); + NormRelVortEdge, Geom, Mesh, OnEdge); Array2DReal NormPlanetVortEdge("NormPlanetVortEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normPlanetVort(X, Y); }, - NormPlanetVortEdge, Geom, Mesh, OnEdge, NVertLevels); + NormPlanetVortEdge, Geom, Mesh, OnEdge); Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThick(X, Y); }, - LayerThickEdge, Geom, Mesh, OnEdge, NVertLevels); + LayerThickEdge, Geom, Mesh, OnEdge); Array2DReal NormVelEdge("NormVelEdge", Mesh->NEdgesSize, NVertLevels); @@ -380,7 +435,7 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, - NormVelEdge, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + NormVelEdge, EdgeComponent::Normal, Geom, Mesh); // Compute numerical result Array2DReal NumPotVortHAdv("NumPotVortHAdv", Mesh->NEdgesOwned, NVertLevels); @@ -395,7 +450,7 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures PotVortHAdvErrors; Err += computeErrors(PotVortHAdvErrors, NumPotVortHAdv, ExactPotVortHAdv, - Mesh, OnEdge, NVertLevels); + Mesh, OnEdge); // Check error values Err += checkErrors("TendencyTermsTest", "PotVortHAdv", PotVortHAdvErrors, @@ -423,15 +478,14 @@ int testKEGrad(int NVertLevels, Real RTol) { VecField[0] = -Setup.gradX(X, Y); VecField[1] = -Setup.gradY(X, Y); }, - ExactKEGrad, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactKEGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Set input array Array2DReal KECell("KECell", Mesh->NCellsSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, KECell, - Geom, Mesh, OnCell, NVertLevels); + Geom, Mesh, OnCell); // Compute numerical result Array2DReal NumKEGrad("NumKEGrad", Mesh->NEdgesOwned, NVertLevels); @@ -444,8 +498,7 @@ int testKEGrad(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures KEGradErrors; - Err += computeErrors(KEGradErrors, NumKEGrad, ExactKEGrad, Mesh, OnEdge, - NVertLevels); + Err += computeErrors(KEGradErrors, NumKEGrad, ExactKEGrad, Mesh, OnEdge); // Check error values Err += checkErrors("TendencyTermsTest", "KEGrad", KEGradErrors, @@ -473,15 +526,14 @@ int testSSHGrad(int NVertLevels, Real RTol) { VecField[0] = -9.80665_Real * Setup.gradX(X, Y); VecField[1] = -9.80665_Real * Setup.gradY(X, Y); }, - ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Set input array Array2DReal SSHCell("SSHCell", Mesh->NCellsSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, SSHCell, - Geom, Mesh, OnCell, NVertLevels); + Geom, Mesh, OnCell); // Compute numerical result Array2DReal NumSSHGrad("NumSSHGrad", Mesh->NEdgesOwned, NVertLevels); @@ -494,8 +546,7 @@ int testSSHGrad(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures SSHGradErrors; - Err += computeErrors(SSHGradErrors, NumSSHGrad, ExactSSHGrad, Mesh, OnEdge, - NVertLevels); + Err += computeErrors(SSHGradErrors, NumSSHGrad, ExactSSHGrad, Mesh, OnEdge); // Check error values Err += checkErrors("TendencyTermsTest", "SSHGrad", SSHGradErrors, @@ -531,21 +582,20 @@ int testVelDiff(int NVertLevels, Real RTol) { VecField[0] = ViscDel2 * Setup.laplaceVecX(X, Y); VecField[1] = ViscDel2 * Setup.laplaceVecY(X, Y); }, - ExactVelDiff, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactVelDiff, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Set input arrays Array2DReal DivCell("DivCell", Mesh->NCellsSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, - DivCell, Geom, Mesh, OnCell, NVertLevels); + DivCell, Geom, Mesh, OnCell); Array2DReal RVortVertex("RVortVertex", Mesh->NVerticesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.curl(X, Y); }, RVortVertex, - Geom, Mesh, OnVertex, NVertLevels); + Geom, Mesh, OnVertex); // Compute numerical result Array2DReal NumVelDiff("NumVelDiff", Mesh->NEdgesOwned, NVertLevels); @@ -557,8 +607,7 @@ int testVelDiff(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures VelDiffErrors; - Err += computeErrors(VelDiffErrors, NumVelDiff, ExactVelDiff, Mesh, OnEdge, - NVertLevels); + Err += computeErrors(VelDiffErrors, NumVelDiff, ExactVelDiff, Mesh, OnEdge); // Check error values Err += checkErrors("TendencyTermsTest", "VelDiff", VelDiffErrors, @@ -602,21 +651,20 @@ int testVelHyperDiff(int NVertLevels, Real RTol) { VecField[0] = -ViscDel4 * Setup.laplaceVecX(X, Y); VecField[1] = -ViscDel4 * Setup.laplaceVecY(X, Y); }, - ExactVelHyperDiff, EdgeComponent::Normal, Geom, Mesh, NVertLevels, - ExchangeHalos::No); + ExactVelHyperDiff, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); // Set input arrays Array2DReal DivCell("DivCell", Mesh->NCellsSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, - DivCell, Geom, Mesh, OnCell, NVertLevels); + DivCell, Geom, Mesh, OnCell); Array2DReal RVortVertex("RVortVertex", Mesh->NVerticesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.curl(X, Y); }, RVortVertex, - Geom, Mesh, OnVertex, NVertLevels); + Geom, Mesh, OnVertex); // Compute numerical result Array2DReal NumVelHyperDiff("NumVelHyperDiff", Mesh->NEdgesOwned, @@ -630,7 +678,7 @@ int testVelHyperDiff(int NVertLevels, Real RTol) { // Compute errors ErrorMeasures VelHyperDiffErrors; Err += computeErrors(VelHyperDiffErrors, NumVelHyperDiff, ExactVelHyperDiff, - Mesh, OnEdge, NVertLevels); + Mesh, OnEdge); // Check error values Err += checkErrors("TendencyTermsTest", "VelHyperDiff", VelHyperDiffErrors, @@ -643,6 +691,148 @@ int testVelHyperDiff(int NVertLevels, Real RTol) { return Err; } // end testVelHyperDiff +int testWindForcing(int NVertLevels) { + + int Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + + const Real SaltWaterDensity = 0.987654321; + + // Compute exact result + Array2DReal ExactWindForcing("ExactWindForcing", Mesh->NEdgesOwned, + NVertLevels); + + // Note: this computes wind forcing at every level + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.windForcingX(X, Y, SaltWaterDensity); + VecField[1] = Setup.windForcingY(X, Y, SaltWaterDensity); + }, + ExactWindForcing, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); + + // Reset wind forcing to zero below the surface + deepCopy(Kokkos::subview(ExactWindForcing, Kokkos::ALL, + Kokkos::make_pair(1, NVertLevels)), + 0); + + // Set input arrays + Array1DReal NormalStressEdge("NormalStressEdge", Mesh->NEdgesSize); + + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.vectorX(X, Y); + VecField[1] = Setup.vectorY(X, Y); + }, + NormalStressEdge, EdgeComponent::Normal, Geom, Mesh); + + Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + LayerThickEdge, Geom, Mesh, OnEdge); + + // Compute numerical result + Array2DReal NumWindForcing("NumWindForcing", Mesh->NEdgesOwned, NVertLevels); + + WindForcingOnEdge WindForcingOnE; + WindForcingOnE.SaltWaterDensity = SaltWaterDensity; + + parallelFor( + {Mesh->NEdgesOwned, NVertLevels}, KOKKOS_LAMBDA(int IEdge, int KLevel) { + WindForcingOnE(NumWindForcing, IEdge, KLevel, NormalStressEdge, + LayerThickEdge); + }); + + // Compute errors + ErrorMeasures WindForcingErrors; + Err += computeErrors(WindForcingErrors, NumWindForcing, ExactWindForcing, + Mesh, OnEdge); + + // Check error values + const Real RTol = 0; + const Real ATol = 100 * std::numeric_limits::epsilon(); + Err += checkErrors("TendencyTermsTest", "WindForcing", WindForcingErrors, + Setup.ExpectedWindForcingErrors, RTol, ATol); + + return Err; +} // end testWindForcing + +int testBottomDrag(int NVertLevels, Real RTol) { + + int Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + + const Real Coeff = 1.123456789; + + // Compute exact result + Array2DReal ExactBottomDrag("ExactBottomDrag", Mesh->NEdgesOwned, + NVertLevels); + + // Note: this computes bottom drag at every level + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.bottomDragX(X, Y, Coeff); + VecField[1] = Setup.bottomDragY(X, Y, Coeff); + }, + ExactBottomDrag, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); + + // Reset bottom drag to zero above the lowest level + deepCopy(Kokkos::subview(ExactBottomDrag, Kokkos::ALL, + Kokkos::make_pair(0, NVertLevels - 1)), + 0); + + // Set input arrays + Array2DReal NormalVelEdge("NormalVelEdge", Mesh->NEdgesSize, NVertLevels); + + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.vectorX(X, Y); + VecField[1] = Setup.vectorY(X, Y); + }, + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh); + + Array2DReal KECell("KECell", Mesh->NCellsSize, NVertLevels); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { + return Setup.scalarA(X, Y) * Setup.scalarA(X, Y) / 2; + }, + KECell, Geom, Mesh, OnCell); + + Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + LayerThickEdge, Geom, Mesh, OnEdge); + + // Compute numerical result + Array2DReal NumBottomDrag("NumBottomDrag", Mesh->NEdgesOwned, NVertLevels); + + BottomDragOnEdge BottomDragOnE(Mesh); + BottomDragOnE.Coeff = Coeff; + + parallelFor( + {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { + BottomDragOnE(NumBottomDrag, IEdge, NormalVelEdge, KECell, + LayerThickEdge); + }); + + // Compute errors + ErrorMeasures BottomDragErrors; + Err += computeErrors(BottomDragErrors, NumBottomDrag, ExactBottomDrag, Mesh, + OnEdge); + + // Check error values + Err += checkErrors("TendencyTermsTest", "BottomDrag", BottomDragErrors, + Setup.ExpectedBottomDragErrors, RTol); + + return Err; +} // end testBottomDrag + int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { I4 Err = 0; @@ -656,8 +846,7 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracerFluxDiv(X, Y); }, - ExactTrFluxDiv, Geom, Mesh, OnCell, NVertLevels, NTracers, - ExchangeHalos::No); + ExactTrFluxDiv, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays Array2DReal NormalVelocity("NormalVelocity", Mesh->NEdgesSize, NVertLevels); @@ -667,13 +856,13 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, - NormalVelocity, EdgeComponent::Normal, Geom, Mesh, NVertLevels); + NormalVelocity, EdgeComponent::Normal, Geom, Mesh); Array3DReal HTrOnEdge("HTrOnEdge", NTracers, Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.layerThick(X, Y); }, - HTrOnEdge, Geom, Mesh, OnEdge, NVertLevels, NTracers); + HTrOnEdge, Geom, Mesh, OnEdge); // Compute numerical result Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, @@ -687,8 +876,8 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { }); ErrorMeasures TrHAdvErrors; - Err += computeErrors(TrHAdvErrors, NumTrFluxDiv, ExactTrFluxDiv, Mesh, - OnCell, NVertLevels, NTracers); + Err += + computeErrors(TrHAdvErrors, NumTrFluxDiv, ExactTrFluxDiv, Mesh, OnCell); Err += checkErrors("TendencyTermsTest", "TracerHorzAdv", TrHAdvErrors, Setup.ExpectedTrHAdvErrors, RTol); @@ -713,8 +902,7 @@ int testTracerDiffOnCell(int NVertLevels, int NTracers, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracerDiff(X, Y); }, - ExactTracerDiff, Geom, Mesh, OnCell, NVertLevels, NTracers, - ExchangeHalos::No); + ExactTracerDiff, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays Array3DReal TracerCell("TracerCell", NTracers, Mesh->NCellsSize, @@ -722,13 +910,13 @@ int testTracerDiffOnCell(int NVertLevels, int NTracers, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, - TracerCell, Geom, Mesh, OnCell, NVertLevels, NTracers); + TracerCell, Geom, Mesh, OnCell); Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, - LayerThickEdge, Geom, Mesh, OnEdge, NVertLevels); + LayerThickEdge, Geom, Mesh, OnEdge); // Compute numerical result Array3DReal NumTracerDiff("NumTracerDiff", NTracers, Mesh->NCellsOwned, @@ -745,7 +933,7 @@ int testTracerDiffOnCell(int NVertLevels, int NTracers, Real RTol) { ErrorMeasures TrDiffErrors; Err += computeErrors(TrDiffErrors, NumTracerDiff, ExactTracerDiff, Mesh, - OnCell, NVertLevels, NTracers); + OnCell); Err += checkErrors("TendencyTermsTest", "TracerDiff", TrDiffErrors, Setup.ExpectedTrDel2Errors, RTol); @@ -770,8 +958,7 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.tracerHyperDiff(X, Y); }, - ExactTracerHyperDiff, Geom, Mesh, OnCell, NVertLevels, NTracers, - ExchangeHalos::No); + ExactTracerHyperDiff, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays Array3DReal TrDel2Cell("TracerCell", NTracers, Mesh->NCellsSize, @@ -779,7 +966,7 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarC(X, Y); }, - TrDel2Cell, Geom, Mesh, OnCell, NVertLevels, NTracers); + TrDel2Cell, Geom, Mesh, OnCell); // Compute numerical result Array3DReal NumTracerHyperDiff("NumTracerHyperDiff", NTracers, @@ -793,9 +980,8 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { }); ErrorMeasures TrHyperDiffErrors; - Err += - computeErrors(TrHyperDiffErrors, NumTracerHyperDiff, - ExactTracerHyperDiff, Mesh, OnCell, NVertLevels, NTracers); + Err += computeErrors(TrHyperDiffErrors, NumTracerHyperDiff, + ExactTracerHyperDiff, Mesh, OnCell); Err += checkErrors("TendencyTermsTest", "TracerHyperDiff", TrHyperDiffErrors, Setup.ExpectedTrDel4Errors, RTol); @@ -807,9 +993,9 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { return Err; } // end testTracerHyperDiffOnCell -int initTendTest(const std::string &mesh) { +void initTendTest(const std::string &MeshFile, int NVertLevels) { - I4 Err = 0; + Error Err; MachEnv::init(MPI_COMM_WORLD); MachEnv *DefEnv = MachEnv::getDefault(); @@ -822,23 +1008,29 @@ int initTendTest(const std::string &mesh) { Config("Omega"); Config::readAll("omega.yml"); + // Reset NVertLevels to the test value + Config *OmegaConfig = Config::getOmegaConfig(); + Config DimConfig("Dimension"); + Err += OmegaConfig->get(DimConfig); + CHECK_ERROR_ABORT(Err, + "TendencyTermsTest: Dimension group not found in Config"); + + DimConfig.set("NVertLevels", NVertLevels); + I4 IOErr = IO::init(DefComm); if (IOErr != 0) { - Err++; - LOG_ERROR("TendencyTermsTest: error initializing parallel IO"); + ABORT_ERROR("TendencyTermsTest: error initializing parallel IO"); } - Decomp::init(mesh); + Decomp::init(MeshFile); int HaloErr = Halo::init(); if (HaloErr != 0) { - Err++; - LOG_ERROR("TendencyTermsTest: error initializing default halo"); + ABORT_ERROR("TendencyTermsTest: error initializing default halo"); } HorzMesh::init(); - return Err; } // end initTendTest void finalizeTendTest() { @@ -851,15 +1043,13 @@ void finalizeTendTest() { } // end finalizeTendTest -int tendencyTermsTest(const std::string &mesh = DefaultMeshFile) { +int tendencyTermsTest(const std::string &MeshFile = DefaultMeshFile) { + int Err = 0; + int NVertLevels = 16; - int Err = initTendTest(mesh); - if (Err != 0) { - LOG_CRITICAL("TendencyTermsTest: Error initializing"); - } + initTendTest(MeshFile, NVertLevels); const auto &Mesh = HorzMesh::getDefault(); - int NVertLevels = 16; int NTracers = 3; const Real RTol = sizeof(Real) == 4 ? 2e-2 : 1e-5; @@ -876,6 +1066,10 @@ int tendencyTermsTest(const std::string &mesh = DefaultMeshFile) { Err += testVelHyperDiff(NVertLevels, RTol); + Err += testWindForcing(NVertLevels); + + Err += testBottomDrag(NVertLevels, RTol); + Err += testTracerHorzAdvOnCell(NVertLevels, NTracers, RTol); Err += testTracerDiffOnCell(NVertLevels, NTracers, RTol); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 3284e46e1214..2f5e8203e5b1 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -28,6 +28,7 @@ #include "MachEnv.h" #include "OceanState.h" #include "OmegaKokkos.h" +#include "Pacer.h" #include "TendencyTerms.h" #include "TimeMgr.h" #include "Tracers.h" @@ -132,8 +133,7 @@ ErrorMeasures computeErrors() { // Only velocity errors matters, because thickness remains constant ErrorMeasures VelErrors; - computeErrors(VelErrors, NormalVelEdge, ExactNormalVelEdge, DefMesh, OnEdge, - NVertLevels); + computeErrors(VelErrors, NormalVelEdge, ExactNormalVelEdge, DefMesh, OnEdge); return VelErrors; } @@ -220,8 +220,11 @@ int initTimeStepperTest(const std::string &mesh) { LOG_ERROR("TimeStepperTest: error creating test state"); } - auto *TestAuxState = - AuxiliaryState::create("TestAuxState", DefMesh, NVertLevels, NTracers); + auto *TestAuxState = AuxiliaryState::create("TestAuxState", DefMesh, DefHalo, + NVertLevels, NTracers); + + TestAuxState->readConfigOptions(OmegaConfig); + if (!TestAuxState) { Err++; LOG_ERROR("TimeStepperTest: error creating test auxiliary state"); @@ -248,6 +251,8 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerHorzAdv.Enabled = false; TestTendencies->TracerDiffusion.Enabled = false; TestTendencies->TracerHyperDiff.Enabled = false; + TestTendencies->WindForcing.Enabled = false; + TestTendencies->BottomDrag.Enabled = false; return Err; } @@ -408,6 +413,8 @@ int main(int argc, char *argv[]) { MPI_Init(&argc, &argv); Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); LOG_INFO("----- Time Stepper Unit Test -----");