diff --git a/components/omega/src/infra/OmegaKokkos.h b/components/omega/src/infra/OmegaKokkos.h index 537e2312204f..5a3eb0b82499 100644 --- a/components/omega/src/infra/OmegaKokkos.h +++ b/components/omega/src/infra/OmegaKokkos.h @@ -359,8 +359,7 @@ KOKKOS_FUNCTION void parallelForInner(const TeamMember &Team, int UpperBound, template struct AccumTypeHelper; template -struct AccumTypeHelper< - T, std::enable_if_t>>> { +struct AccumTypeHelper>> { using Type = T; }; @@ -386,7 +385,8 @@ inline void parallelReduceOuter(const std::string &Label, auto Policy = TeamPolicy(LinBound, OMEGA_TEAMSIZE); Kokkos::parallel_reduce( Label, Policy, - KOKKOS_LAMBDA(const TeamMember &Team, AccumType &...Accums) { + KOKKOS_LAMBDA(const TeamMember &Team, + AccumType> &...Accums) { const int TeamId = Team.league_rank(); LinFunctor(TeamId, Team, Accums...); }, diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 64cdad749f25..9643fe264cd5 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -18,14 +18,15 @@ static std::string stripDefault(const std::string &Name) { // Constructor. Constructs the member auxiliary variables and registers their // fields with IOStreams AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, - Halo *MeshHalo, int NVertLayers, int NTracers) - : Mesh(Mesh), MeshHalo(MeshHalo), Name(stripDefault(Name)), - KineticAux(stripDefault(Name), Mesh, NVertLayers), - LayerThicknessAux(stripDefault(Name), Mesh, NVertLayers), - VorticityAux(stripDefault(Name), Mesh, NVertLayers), - VelocityDel2Aux(stripDefault(Name), Mesh, NVertLayers), + Halo *MeshHalo, const VertCoord *VCoord, + int NTracers) + : Mesh(Mesh), MeshHalo(MeshHalo), VCoord(VCoord), Name(stripDefault(Name)), + KineticAux(stripDefault(Name), Mesh, VCoord), + LayerThicknessAux(stripDefault(Name), Mesh, VCoord), + VorticityAux(stripDefault(Name), Mesh, VCoord), + VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), WindForcingAux(stripDefault(Name), Mesh), - TracerAux(stripDefault(Name), Mesh, NVertLayers, NTracers) { + TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -64,31 +65,53 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, State->getLayerThickness(LayerThickCell, ThickTimeLevel); State->getNormalVelocity(NormalVelEdge, VelTimeLevel); - const int NVertLayers = LayerThickCell.extent_int(1); - const int NChunks = NVertLayers / VecLength; - OMEGA_SCOPE(LocKineticAux, KineticAux); OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocVorticityAux, VorticityAux); OMEGA_SCOPE(LocVelocityDel2Aux, VelocityDel2Aux); OMEGA_SCOPE(LocWindForcingAux, WindForcingAux); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot); + OMEGA_SCOPE(MinLayerVertexTop, VCoord->MinLayerVertexTop); + OMEGA_SCOPE(MaxLayerVertexBot, VCoord->MaxLayerVertexBot); + OMEGA_SCOPE(MaxLayerVertexTop, VCoord->MaxLayerVertexTop); + OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + Pacer::start("AuxState:computeMomAux", 1); Pacer::start("AuxState:vertexAuxState1", 2); - parallelFor( - "vertexAuxState1", {Mesh->NVerticesAll, NChunks}, - KOKKOS_LAMBDA(int IVertex, int KChunk) { - LocVorticityAux.computeVarsOnVertex(IVertex, KChunk, LayerThickCell, - NormalVelEdge); + parallelForOuter( + "vertexAuxState1", {Mesh->NVerticesAll}, + KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) { + const int KMin = MinLayerVertexTop(IVertex); + const int KMax = MaxLayerVertexBot(IVertex); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVorticityAux.computeVarsOnVertex( + IVertex, KChunk, LayerThickCell, NormalVelEdge); + }); }); Pacer::stop("AuxState:vertexAuxState1", 2); Pacer::start("AuxState:cellAuxState1", 2); - parallelFor( - "cellAuxState1", {Mesh->NCellsAll, NChunks}, - KOKKOS_LAMBDA(int ICell, int KChunk) { - LocKineticAux.computeVarsOnCell(ICell, KChunk, NormalVelEdge); + parallelForOuter( + "cellAuxState1", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocKineticAux.computeVarsOnCell(ICell, KChunk, NormalVelEdge); + }); }); Pacer::stop("AuxState:cellAuxState1", 2); @@ -103,39 +126,79 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, Pacer::stop("AuxState:edgeAuxState1", 2); Pacer::start("AuxState:edgeAuxState2", 2); - parallelFor( - "edgeAuxState2", {Mesh->NEdgesAll, NChunks}, - KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocVorticityAux.computeVarsOnEdge(IEdge, KChunk); - LocLayerThicknessAux.computeVarsOnEdge(IEdge, KChunk, LayerThickCell, - NormalVelEdge); - LocVelocityDel2Aux.computeVarsOnEdge(IEdge, KChunk, VelocityDivCell, - RelVortVertex); + parallelForOuter( + "edgeAuxState2", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocLayerThicknessAux.computeVarsOnEdge( + IEdge, KChunk, LayerThickCell, NormalVelEdge); + LocVelocityDel2Aux.computeVarsOnEdge( + IEdge, KChunk, VelocityDivCell, RelVortVertex); + }); + }); + + parallelForOuter( + "edgeAuxState2", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeTop(IEdge); + const int KMax = MaxLayerEdgeBot(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVorticityAux.computeVarsOnEdge(IEdge, KChunk); + }); }); Pacer::stop("AuxState:edgeAuxState2", 2); Pacer::start("AuxState:vertexAuxState2", 2); - parallelFor( - "vertexAuxState2", {Mesh->NVerticesAll, NChunks}, - KOKKOS_LAMBDA(int IVertex, int KChunk) { - LocVelocityDel2Aux.computeVarsOnVertex(IVertex, KChunk); + parallelForOuter( + "vertexAuxState2", {Mesh->NVerticesAll}, + KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) { + const int KMin = MinLayerVertexBot(IVertex); + const int KMax = MaxLayerVertexTop(IVertex); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVelocityDel2Aux.computeVarsOnVertex(IVertex, KChunk); + }); }); Pacer::stop("AuxState:vertexAuxState2", 2); Pacer::start("AuxState:cellAuxState2", 2); - parallelFor( - "cellAuxState2", {Mesh->NCellsAll, NChunks}, - KOKKOS_LAMBDA(int ICell, int KChunk) { - LocVelocityDel2Aux.computeVarsOnCell(ICell, KChunk); + parallelForOuter( + "cellAuxState2", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVelocityDel2Aux.computeVarsOnCell(ICell, KChunk); + }); }); Pacer::stop("AuxState:cellAuxState2", 2); Pacer::start("AuxState:cellAuxState3", 2); - parallelFor( - "cellAuxState3", {Mesh->NCellsAll, NChunks}, - KOKKOS_LAMBDA(int ICell, int KChunk) { - LocLayerThicknessAux.computeVarsOnCells(ICell, KChunk, - LayerThickCell); + parallelForOuter( + "cellAuxState3", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocLayerThicknessAux.computeVarsOnCells(ICell, KChunk, + LayerThickCell); + }); }); Pacer::stop("AuxState:cellAuxState3", 2); @@ -151,33 +214,49 @@ void AuxiliaryState::computeAll(const OceanState *State, State->getLayerThickness(LayerThickCell, ThickTimeLevel); State->getNormalVelocity(NormalVelEdge, VelTimeLevel); - const int NVertLayers = LayerThickCell.extent_int(1); - const int NChunks = NVertLayers / VecLength; - const int NTracers = TracerArray.extent_int(0); + const int NTracers = TracerArray.extent_int(0); OMEGA_SCOPE(LocTracerAux, TracerAux); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); Pacer::start("AuxState:computeAll", 1); computeMomAux(State, ThickTimeLevel, VelTimeLevel); Pacer::start("AuxState:edgeAuxState4", 2); - parallelFor( - "edgeAuxState4", {NTracers, Mesh->NEdgesAll, NChunks}, - KOKKOS_LAMBDA(int LTracer, int IEdge, int KChunk) { - LocTracerAux.computeVarsOnEdge(LTracer, IEdge, KChunk, NormalVelEdge, - LayerThickCell, TracerArray); + parallelForOuter( + "edgeAuxState4", {NTracers, Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int LTracer, int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerAux.computeVarsOnEdge(LTracer, IEdge, KChunk, + NormalVelEdge, LayerThickCell, + TracerArray); + }); }); Pacer::stop("AuxState:edgeAuxState4", 2); const auto &MeanLayerThickEdge = LayerThicknessAux.MeanLayerThickEdge; Pacer::start("AuxState:cellAuxState4", 2); - parallelFor( - "cellAuxState4", {NTracers, Mesh->NCellsAll, NChunks}, - KOKKOS_LAMBDA(int LTracer, int ICell, int KChunk) { - LocTracerAux.computeVarsOnCells(LTracer, ICell, KChunk, - MeanLayerThickEdge, TracerArray); + parallelForOuter( + "cellAuxState4", {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int LTracer, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerAux.computeVarsOnCells( + LTracer, ICell, KChunk, MeanLayerThickEdge, TracerArray); + }); }); Pacer::stop("AuxState:cellAuxState4", 2); @@ -193,7 +272,8 @@ void AuxiliaryState::computeAll(const OceanState *State, // Create a non-default auxiliary state AuxiliaryState *AuxiliaryState::create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - int NVertLayers, const int NTracers) { + const VertCoord *VCoord, + const int NTracers) { if (AllAuxStates.find(Name) != AllAuxStates.end()) { LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it " "already exists", @@ -202,7 +282,7 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, } auto *NewAuxState = - new AuxiliaryState(Name, Mesh, MeshHalo, NVertLayers, NTracers); + new AuxiliaryState(Name, Mesh, MeshHalo, VCoord, NTracers); AllAuxStates.emplace(Name, NewAuxState); return NewAuxState; @@ -211,14 +291,14 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, // Create the default auxiliary state. Assumes that HorzMesh, VertCoord and // Halo have been initialized. void AuxiliaryState::init() { - const HorzMesh *DefMesh = HorzMesh::getDefault(); - Halo *DefHalo = Halo::getDefault(); + const HorzMesh *DefMesh = HorzMesh::getDefault(); + Halo *DefHalo = Halo::getDefault(); + const VertCoord *DefVCoord = VertCoord::getDefault(); - int NVertLayers = VertCoord::getDefault()->NVertLayers; - int NTracers = Tracers::getNumTracers(); + int NTracers = Tracers::getNumTracers(); - AuxiliaryState::DefaultAuxState = AuxiliaryState::create( - "Default", DefMesh, DefHalo, NVertLayers, NTracers); + AuxiliaryState::DefaultAuxState = + AuxiliaryState::create("Default", DefMesh, DefHalo, DefVCoord, NTracers); Config *OmegaConfig = Config::getOmegaConfig(); DefaultAuxState->readConfigOptions(OmegaConfig); diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 92eee116e6ca..313cc157f0f3 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -50,7 +50,8 @@ class AuxiliaryState { // Create a non-default auxiliary state static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh, - Halo *MeshHalo, int NVertLayers, int NTracers); + Halo *MeshHalo, const VertCoord *VCoord, + int NTracers); /// Get the default auxiliary state static AuxiliaryState *getDefault(); @@ -83,13 +84,14 @@ class AuxiliaryState { private: AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - int NVertLayers, int NTracers); + const VertCoord *VCoord, int NTracers); AuxiliaryState(const AuxiliaryState &) = delete; AuxiliaryState(AuxiliaryState &&) = delete; const HorzMesh *Mesh; Halo *MeshHalo; + const VertCoord *VCoord; static AuxiliaryState *DefaultAuxState; static std::map> AllAuxStates; }; diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index b8d1751c24c3..2530583e738c 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -14,26 +14,25 @@ namespace OMEGA { /// Constructor for Teos10Eos -Teos10Eos::Teos10Eos(int NVertLayers) : NVertLayers(NVertLayers) { +Teos10Eos::Teos10Eos(const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) { SpecVolPCoeffs = Array2DReal("SpecVolPCoeffs", 6, VecLength); } /// Constructor for LinearEos -LinearEos::LinearEos() {} +LinearEos::LinearEos(const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} /// Constructor for Eos -Eos::Eos(const std::string &Name_, ///< [in] Name for eos object - const HorzMesh *Mesh, ///< [in] Horizontal mesh - int NVertLayers ///< [in] Number of vertical layers +Eos::Eos(const std::string &Name, ///< [in] Name for eos object + const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord ///< [in] Vertical coordinate ) - : ComputeSpecVolTeos10(NVertLayers) { - SpecVol = Array2DReal("SpecVol", Mesh->NCellsAll, NVertLayers); + : ComputeSpecVolLinear(VCoord), ComputeSpecVolTeos10(VCoord), Name(Name), + Mesh(Mesh), VCoord(VCoord) { + SpecVol = Array2DReal("SpecVol", Mesh->NCellsAll, VCoord->NVertLayers); SpecVolDisplaced = - Array2DReal("SpecVolDisplaced", Mesh->NCellsAll, NVertLayers); - // Array dimension lengths - NCellsAll = Mesh->NCellsAll; - NChunks = NVertLayers / VecLength; - Name = Name_; + Array2DReal("SpecVolDisplaced", Mesh->NCellsAll, VCoord->NVertLayers); defineFields(); } @@ -60,8 +59,8 @@ void Eos::destroyInstance() { void Eos::init() { if (!Instance) { - Instance = new Eos("Default", HorzMesh::getDefault(), - VertCoord::getDefault()->NVertLayers); + Instance = + new Eos("Default", HorzMesh::getDefault(), VertCoord::getDefault()); } Error Err; // error code @@ -118,24 +117,42 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp, ComputeSpecVolLinear); /// Local view for linear EOS computation OMEGA_SCOPE(LocComputeSpecVolTeos10, ComputeSpecVolTeos10); /// Local view for TEOS-10 computation + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + deepCopy(LocSpecVol, 0); /// Initialize local specific volume to zero I4 KDisp = 0; /// No displacement in this case /// Dispatch to the correct EOS calculation if (EosChoice == EosType::LinearEos) { - parallelFor( - "eos-linear", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { - LocComputeSpecVolLinear(LocSpecVol, ICell, KChunk, ConservTemp, - AbsSalinity); + parallelForOuter( + "eos-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeSpecVolLinear(LocSpecVol, ICell, KChunk, + ConservTemp, AbsSalinity); + }); }); } else if (EosChoice == EosType::Teos10Eos) { - parallelFor( - "eos-teos10", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { - LocComputeSpecVolTeos10(LocSpecVol, ICell, KChunk, ConservTemp, - AbsSalinity, Pressure, KDisp); + parallelForOuter( + "eos-teos10", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeSpecVolTeos10(LocSpecVol, ICell, KChunk, + ConservTemp, AbsSalinity, Pressure, + KDisp); + }); }); } } @@ -150,6 +167,9 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, ComputeSpecVolLinear); /// Local view for linear EOS computation OMEGA_SCOPE(LocComputeSpecVolTeos10, ComputeSpecVolTeos10); /// Local view for TEOS-10 computation + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + deepCopy(LocSpecVolDisplaced, 0); /// Initialize local specific volume to zero @@ -160,18 +180,31 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, LOG_INFO("Eos::computeSpecVolDisp called with Linear EOS. " "SpecVol is independent of pressure/depth, so the " "displaced value will be the same as SpecVol."); - parallelFor( - "eos-linear", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { - LocComputeSpecVolLinear(LocSpecVolDisplaced, ICell, KChunk, - ConservTemp, AbsSalinity); + parallelForOuter( + "eos-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeSpecVolLinear(LocSpecVolDisplaced, ICell, KChunk, + ConservTemp, AbsSalinity); + }); }); } else if (EosChoice == EosType::Teos10Eos) { - parallelFor( - "eos-teos10", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { - LocComputeSpecVolTeos10(LocSpecVolDisplaced, ICell, KChunk, - ConservTemp, AbsSalinity, Pressure, KDisp); + parallelForOuter( + "eos-teos10", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeSpecVolTeos10(LocSpecVolDisplaced, ICell, KChunk, + ConservTemp, AbsSalinity, Pressure, + KDisp); + }); }); } } diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index b9eb0f949c31..c463a285156f 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -32,7 +32,7 @@ class Teos10Eos { Array2DReal SpecVolPCoeffs; /// constructor declaration - Teos10Eos(int NVertLayers); + Teos10Eos(const VertCoord *VCoord); // The functor takes the full arrays of specific volume (inout), // the indices ICell and KChunk, and the ocean tracers (conservative) @@ -45,8 +45,9 @@ class Teos10Eos { I4 KDisp) const { OMEGA_SCOPE(LocSpecVolPCoeffs, SpecVolPCoeffs); - const I4 KStart = KChunk * VecLength; - for (int KVec = 0; KVec < VecLength; ++KVec) { + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; /// Calculate the local specific volume polynomial pressure @@ -66,8 +67,8 @@ class Teos10Eos { calcDelta(LocSpecVolPCoeffs, KVec, Pressure(ICell, K)); } else { // Displacement, use the displaced pressure - I4 KTmp = Kokkos::min(K + KDisp, NVertLayers - 1); - KTmp = Kokkos::max(0, KTmp); + I4 KTmp = Kokkos::min(K + KDisp, MaxLayerCell(ICell)); + KTmp = Kokkos::max(MinLayerCell(ICell), KTmp); SpecVol(ICell, K) = calcRefProfile(Pressure(ICell, KTmp)) + calcDelta(LocSpecVolPCoeffs, KVec, Pressure(ICell, KTmp)); @@ -237,7 +238,8 @@ class Teos10Eos { } private: - const int NVertLayers; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; /// Linear Equation of State @@ -249,7 +251,7 @@ class LinearEos { Real RhoT0S0 = 1000.0; ///< Reference density (kg m^-3) at (T,S)=(0,0) /// constructor declaration - LinearEos(); + LinearEos(const VertCoord *VCoord); // The functor takes the full arrays of specific volume (inout), // the indices ICell and KChunk, and the ocean tracers (conservative) @@ -258,14 +260,21 @@ class LinearEos { KOKKOS_FUNCTION void operator()(Array2DReal SpecVol, I4 ICell, I4 KChunk, const Array2DReal &ConservTemp, const Array2DReal &AbsSalinity) const { - const I4 KStart = KChunk * VecLength; - for (int KVec = 0; KVec < VecLength; ++KVec) { + + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; SpecVol(ICell, K) = 1.0_Real / (RhoT0S0 + (DRhodT * ConservTemp(ICell, K) + DRhodS * AbsSalinity(ICell, K))); } } + + private: + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; /// Class for Equation of State (EOS) calculations @@ -302,7 +311,7 @@ class Eos { private: /// Private constructor - Eos(const std::string &Name, const HorzMesh *Mesh, int NVertLayers); + Eos(const std::string &Name, const HorzMesh *Mesh, const VertCoord *VCoord); /// Private destructor ~Eos(); @@ -316,8 +325,8 @@ class Eos { Eos(Eos &&) = delete; Eos &operator=(Eos &&) = delete; - I4 NCellsAll; ///< Number of horizontal cells - I4 NChunks; ///< Number of vertical chunks (for vectorization) + const HorzMesh *Mesh; ///< Horizontal mesh + const VertCoord *VCoord; ///< Vertical coordinate Teos10Eos ComputeSpecVolTeos10; ///< TEOS-10 specific volume calculator LinearEos ComputeSpecVolLinear; ///< Linear specific volume calculator diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 7d3efab243bd..a59adfe0c6a0 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -221,10 +221,12 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend) - : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh), KEGrad(Mesh), - SSHGrad(Mesh), VelocityDiffusion(Mesh), VelocityHyperDiff(Mesh), - WindForcing(Mesh), BottomDrag(Mesh, VCoord), TracerHorzAdv(Mesh), - TracerDiffusion(Mesh), TracerHyperDiff(Mesh), + : Mesh(Mesh), VCoord(VCoord), ThicknessFluxDiv(Mesh, VCoord), + PotientialVortHAdv(Mesh, VCoord), KEGrad(Mesh, VCoord), + SSHGrad(Mesh, VCoord), VelocityDiffusion(Mesh, VCoord), + VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), + BottomDrag(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), + TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend) { @@ -236,11 +238,7 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies TracerTend = Array3DReal("TracerTend", NTracersIn, Mesh->NCellsSize, VCoord->NVertLayers); - // Array dimension lengths - NCellsAll = Mesh->NCellsAll; - NEdgesAll = Mesh->NEdgesAll; - NTracers = NTracersIn; - NChunks = VCoord->NVertLayers / VecLength; + NTracers = NTracersIn; } // end constructor @@ -264,12 +262,26 @@ void Tendencies::computeThicknessTendenciesOnly( OMEGA_SCOPE(LocLayerThicknessTend, LayerThicknessTend); OMEGA_SCOPE(LocThicknessFluxDiv, ThicknessFluxDiv); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + Array2DReal NormalVelEdge; State->getNormalVelocity(NormalVelEdge, VelTimeLevel); Pacer::start("Tend:computeThicknessTendenciesOnly", 1); - deepCopy(LocLayerThicknessTend, 0); + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LocLayerThicknessTend(ICell, K) = 0; + }); + }); // Compute thickness flux divergence const Array2DReal &ThickFluxEdge = @@ -277,10 +289,17 @@ void Tendencies::computeThicknessTendenciesOnly( if (LocThicknessFluxDiv.Enabled) { Pacer::start("Tend:thicknessFluxDiv", 2); - parallelFor( - {NCellsAll, NChunks}, KOKKOS_LAMBDA(int ICell, int KChunk) { - LocThicknessFluxDiv(LocLayerThicknessTend, ICell, KChunk, - ThickFluxEdge, NormalVelEdge); + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocThicknessFluxDiv(LocLayerThicknessTend, ICell, KChunk, + ThickFluxEdge, NormalVelEdge); + }); }); Pacer::stop("Tend:thicknessFluxDiv", 2); } @@ -314,10 +333,23 @@ void Tendencies::computeVelocityTendenciesOnly( OMEGA_SCOPE(LocVelocityHyperDiff, VelocityHyperDiff); OMEGA_SCOPE(LocWindForcing, WindForcing); OMEGA_SCOPE(LocBottomDrag, BottomDrag); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); Pacer::start("Tend:computeVelocityTendenciesOnly", 1); - deepCopy(LocNormalVelocityTend, 0); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LocNormalVelocityTend(IEdge, K) = 0; + }); + }); const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; @@ -330,11 +362,18 @@ void Tendencies::computeVelocityTendenciesOnly( State->getNormalVelocity(NormVelEdge, VelTimeLevel); if (LocPotientialVortHAdv.Enabled) { Pacer::start("Tend:potientialVortHAdv", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocPotientialVortHAdv(LocNormalVelocityTend, IEdge, KChunk, - NormRVortEdge, NormFEdge, FluxLayerThickEdge, - NormVelEdge); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocPotientialVortHAdv(LocNormalVelocityTend, IEdge, KChunk, + NormRVortEdge, NormFEdge, + FluxLayerThickEdge, NormVelEdge); + }); }); Pacer::stop("Tend:potientialVortHAdv", 2); } @@ -343,9 +382,15 @@ void Tendencies::computeVelocityTendenciesOnly( const Array2DReal &KECell = AuxState->KineticAux.KineticEnergyCell; if (LocKEGrad.Enabled) { Pacer::start("Tend:KEGrad", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocKEGrad(LocNormalVelocityTend, IEdge, KChunk, KECell); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocKEGrad(LocNormalVelocityTend, IEdge, KChunk, KECell); + }); }); Pacer::stop("Tend:KEGrad", 2); } @@ -354,9 +399,15 @@ void Tendencies::computeVelocityTendenciesOnly( const Array2DReal &SSHCell = AuxState->LayerThicknessAux.SshCell; if (LocSSHGrad.Enabled) { Pacer::start("Tend:SSHGrad", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocSSHGrad(LocNormalVelocityTend, IEdge, KChunk, SSHCell); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocSSHGrad(LocNormalVelocityTend, IEdge, KChunk, SSHCell); + }); }); Pacer::stop("Tend:SSHGrad", 2); } @@ -366,10 +417,16 @@ void Tendencies::computeVelocityTendenciesOnly( const Array2DReal &RVortVertex = AuxState->VorticityAux.RelVortVertex; if (LocVelocityDiffusion.Enabled) { Pacer::start("Tend:velocityDiffusion", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocVelocityDiffusion(LocNormalVelocityTend, IEdge, KChunk, DivCell, - RVortVertex); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVelocityDiffusion(LocNormalVelocityTend, IEdge, KChunk, + DivCell, RVortVertex); + }); }); Pacer::stop("Tend:velocityDiffusion", 2); } @@ -380,10 +437,16 @@ void Tendencies::computeVelocityTendenciesOnly( AuxState->VelocityDel2Aux.Del2RelVortVertex; if (LocVelocityHyperDiff.Enabled) { Pacer::start("Tend:velocityHyperDiff", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocVelocityHyperDiff(LocNormalVelocityTend, IEdge, KChunk, - Del2DivCell, Del2RVortVertex); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocVelocityHyperDiff(LocNormalVelocityTend, IEdge, KChunk, + Del2DivCell, Del2RVortVertex); + }); }); Pacer::stop("Tend:velocityHyperDiff", 2); } @@ -394,10 +457,16 @@ void Tendencies::computeVelocityTendenciesOnly( AuxState->LayerThicknessAux.MeanLayerThickEdge; if (LocWindForcing.Enabled) { Pacer::start("Tend:windForcing", 2); - parallelFor( - {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { - LocWindForcing(LocNormalVelocityTend, IEdge, KChunk, - NormalStressEdge, MeanLayerThickEdge); + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocWindForcing(LocNormalVelocityTend, IEdge, KChunk, + NormalStressEdge, MeanLayerThickEdge); + }); }); Pacer::stop("Tend:windForcing", 2); } @@ -406,7 +475,7 @@ void Tendencies::computeVelocityTendenciesOnly( if (LocBottomDrag.Enabled) { Pacer::start("Tend:bottomDrag", 2); parallelFor( - {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { LocBottomDrag(LocNormalVelocityTend, IEdge, NormalVelEdge, KECell, MeanLayerThickEdge); }); @@ -436,21 +505,41 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerHorzAdv, TracerHorzAdv); OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); Pacer::start("Tend:computeTracerTendenciesOnly", 1); - deepCopy(LocTracerTend, 0); + parallelForOuter( + {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LocTracerTend(L, ICell, K) = 0; + }); + }); // compute tracer horizotal advection const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; const Array3DReal &HTracersEdge = AuxState->TracerAux.HTracersEdge; if (LocTracerHorzAdv.Enabled) { Pacer::start("Tend:tracerHorzAdv", 2); - parallelFor( - {NTracers, NCellsAll, NChunks}, - KOKKOS_LAMBDA(int L, int ICell, int KChunk) { - LocTracerHorzAdv(LocTracerTend, L, ICell, KChunk, NormalVelEdge, - HTracersEdge); + parallelForOuter( + {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv(LocTracerTend, L, ICell, KChunk, + NormalVelEdge, HTracersEdge); + }); }); Pacer::stop("Tend:tracerHorzAdv", 2); } @@ -460,11 +549,18 @@ void Tendencies::computeTracerTendenciesOnly( AuxState->LayerThicknessAux.MeanLayerThickEdge; if (LocTracerDiffusion.Enabled) { Pacer::start("Tend:tracerDiffusion", 2); - parallelFor( - {NTracers, NCellsAll, NChunks}, - KOKKOS_LAMBDA(int L, int ICell, int KChunk) { - LocTracerDiffusion(LocTracerTend, L, ICell, KChunk, TracerArray, - MeanLayerThickEdge); + parallelForOuter( + {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerDiffusion(LocTracerTend, L, ICell, KChunk, + TracerArray, MeanLayerThickEdge); + }); }); Pacer::stop("Tend:tracerDiffusion", 2); } @@ -473,11 +569,18 @@ void Tendencies::computeTracerTendenciesOnly( const Array3DReal &Del2TracersCell = AuxState->TracerAux.Del2TracersCell; if (LocTracerHyperDiff.Enabled) { Pacer::start("Tend:tracerHyperDiff", 2); - parallelFor( - {NTracers, NCellsAll, NChunks}, - KOKKOS_LAMBDA(int L, int ICell, int KChunk) { - LocTracerHyperDiff(LocTracerTend, L, ICell, KChunk, - Del2TracersCell); + parallelForOuter( + {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHyperDiff(LocTracerTend, L, ICell, KChunk, + Del2TracersCell); + }); }); Pacer::stop("Tend:tracerHyperDiff", 2); } @@ -500,15 +603,24 @@ void Tendencies::computeThicknessTendencies( OMEGA_SCOPE(LayerThicknessAux, AuxState->LayerThicknessAux); OMEGA_SCOPE(LayerThickCell, LayerThick); OMEGA_SCOPE(NormalVelEdge, NormVel); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); Pacer::start("Tend:computeThicknessTendencies", 1); Pacer::start("Tend:computeLayerThickAux", 2); - parallelFor( - "computeLayerThickAux", {NEdgesAll, NChunks}, - KOKKOS_LAMBDA(int IEdge, int KChunk) { - LayerThicknessAux.computeVarsOnEdge(IEdge, KChunk, LayerThickCell, - NormalVelEdge); + parallelForOuter( + "computeLayerThickAux", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LayerThicknessAux.computeVarsOnEdge( + IEdge, KChunk, LayerThickCell, NormalVelEdge); + }); }); Pacer::stop("Tend:computeLayerThickAux", 2); @@ -545,26 +657,44 @@ void Tendencies::computeTracerTendencies( OMEGA_SCOPE(TracerAux, AuxState->TracerAux); OMEGA_SCOPE(LayerThickCell, State->LayerThickness[ThickTimeLevel]); OMEGA_SCOPE(NormalVelEdge, State->NormalVelocity[VelTimeLevel]); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); Pacer::start("Tend:computeTracerTendencies", 1); Pacer::start("Tend:computeTracerAuxEdge", 2); - parallelFor( - "computeTracerAuxEdge", {NTracers, NEdgesAll, NChunks}, - KOKKOS_LAMBDA(int LTracer, int IEdge, int KChunk) { - TracerAux.computeVarsOnEdge(LTracer, IEdge, KChunk, NormalVelEdge, - LayerThickCell, TracerArray); + parallelForOuter( + "computeTracerAuxEdge", {NTracers, Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int LTracer, int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + TracerAux.computeVarsOnEdge(LTracer, IEdge, KChunk, + NormalVelEdge, LayerThickCell, + TracerArray); + }); }); Pacer::stop("Tend:computeTracerAuxEdge", 2); const auto &MeanLayerThickEdge = AuxState->LayerThicknessAux.MeanLayerThickEdge; Pacer::start("Tend:computeTracerAuxCell", 2); - parallelFor( - "computeTracerAuxCell", {NTracers, NCellsAll, NChunks}, - KOKKOS_LAMBDA(int LTracer, int ICell, int KChunk) { - TracerAux.computeVarsOnCells(LTracer, ICell, KChunk, - MeanLayerThickEdge, TracerArray); + parallelForOuter( + "computeTracerAuxCell", {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int LTracer, int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + TracerAux.computeVarsOnCells(LTracer, ICell, KChunk, + MeanLayerThickEdge, TracerArray); + }); }); Pacer::stop("Tend:computeTracerAuxCell", 2); diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 2544f0b77116..851ce1fc646c 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -167,11 +167,9 @@ class Tendencies { Tendencies(const Tendencies &) = delete; Tendencies(Tendencies &&) = delete; - // Mesh sizes - I4 NCellsAll; ///< Number of cells including full halo - I4 NEdgesAll; ///< Number of edges including full halo - I4 NTracers; ///< Number of tracers - I4 NChunks; ///< Number of vertical layer chunks + const HorzMesh *Mesh; ///< Pointer to horizontal mesh + const VertCoord *VCoord; ///< Pointer to vertical coordinate + I4 NTracers; ///< Number of tracers // Pointer to default tendencies static Tendencies *DefaultTendencies; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index aa4570d1f297..2ccaf6e5236b 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -14,62 +14,90 @@ #include "HorzMesh.h" #include "OceanState.h" #include "Tracers.h" -#include "VertCoord.h" namespace OMEGA { -ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(const HorzMesh *Mesh) +ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), - EdgeSignOnCell(Mesh->EdgeSignOnCell) {} + EdgeSignOnCell(Mesh->EdgeSignOnCell), MinLayerCell(VCoord->MinLayerCell), + MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -PotentialVortHAdvOnEdge::PotentialVortHAdvOnEdge(const HorzMesh *Mesh) +PotentialVortHAdvOnEdge::PotentialVortHAdvOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) : NEdgesOnEdge(Mesh->NEdgesOnEdge), EdgesOnEdge(Mesh->EdgesOnEdge), - WeightsOnEdge(Mesh->WeightsOnEdge), EdgeMask(Mesh->EdgeMask) {} + WeightsOnEdge(Mesh->WeightsOnEdge), EdgeMask(Mesh->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -KEGradOnEdge::KEGradOnEdge(const HorzMesh *Mesh) +KEGradOnEdge::KEGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord) : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge), - EdgeMask(Mesh->EdgeMask) {} + EdgeMask(Mesh->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -SSHGradOnEdge::SSHGradOnEdge(const HorzMesh *Mesh) +SSHGradOnEdge::SSHGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord) : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge), - EdgeMask(Mesh->EdgeMask) {} + EdgeMask(Mesh->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -VelocityDiffusionOnEdge::VelocityDiffusionOnEdge(const HorzMesh *Mesh) +VelocityDiffusionOnEdge::VelocityDiffusionOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), - MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask) {} + MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh) +VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), - MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {} + MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh) - : Enabled(false), EdgeMask(Mesh->EdgeMask) {} +WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) + : Enabled(false), EdgeMask(Mesh->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot) {} BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord) : Enabled(false), Coeff(0), CellsOnEdge(Mesh->CellsOnEdge), - NVertLayers(VCoord->NVertLayers), EdgeMask(Mesh->EdgeMask) {} + NVertLayers(VCoord->NVertLayers), EdgeMask(Mesh->EdgeMask), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh) +TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), - DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), EdgeMask(Mesh->EdgeMask) { -} + DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), EdgeMask(Mesh->EdgeMask), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh) +TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell), - MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask) {} + MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} -TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh) +TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell), - MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {} + MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} } // end namespace OMEGA diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index a6a07574a931..eead7467b7a5 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -28,7 +28,7 @@ class ThicknessFluxDivOnCell { bool Enabled; /// constructor declaration - ThicknessFluxDivOnCell(const HorzMesh *Mesh); + ThicknessFluxDivOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes cell index, vertical chunk index, and thickness flux /// array as inputs, outputs the tendency array @@ -36,23 +36,29 @@ class ThicknessFluxDivOnCell { const Array2DReal &ThicknessFlux, const Array2DReal &NormalVelEdge) const { - const I4 KStart = KChunk * VecLength; + const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const I4 KEndCell = KStartCell + KLenCell - 1; const Real InvAreaCell = 1._Real / AreaCell(ICell); Real DivTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnCell(ICell); ++J) { const I4 JEdge = EdgesOnCell(ICell, J); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + + const I4 KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const I4 KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); + + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const I4 KVec = K - KStartCell; DivTmp[KVec] -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * ThicknessFlux(JEdge, K) * NormalVelEdge(JEdge, K) * InvAreaCell; } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const I4 K = KStartCell + KVec; Tend(ICell, K) -= DivTmp[KVec]; } } @@ -63,6 +69,10 @@ class ThicknessFluxDivOnCell { Array1DReal DvEdge; Array1DReal AreaCell; Array2DReal EdgeSignOnCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Horizontal advection of potential vorticity defined on edges, for @@ -72,7 +82,7 @@ class PotentialVortHAdvOnEdge { bool Enabled; /// constructor declaration - PotentialVortHAdvOnEdge(const HorzMesh *Mesh); + PotentialVortHAdvOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes edge index, vertical chunk index, and arrays for /// normalized relative vorticity, normalized planetary vorticity, layer @@ -84,12 +94,13 @@ class PotentialVortHAdvOnEdge { const Array2DReal &FluxLayerThickEdge, const Array2DReal &NormVelEdge) const { - const I4 KStart = KChunk * VecLength; + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); Real VortTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnEdge(IEdge); ++J) { I4 JEdge = EdgesOnEdge(IEdge, J); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Real NormVort = (NormRVortEdge(IEdge, K) + NormFEdge(IEdge, K) + NormRVortEdge(JEdge, K) + NormFEdge(JEdge, K)) * @@ -101,7 +112,7 @@ class PotentialVortHAdvOnEdge { } } - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Tend(IEdge, K) += EdgeMask(IEdge, K) * VortTmp[KVec]; } @@ -112,6 +123,8 @@ class PotentialVortHAdvOnEdge { Array2DI4 EdgesOnEdge; Array2DReal WeightsOnEdge; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Gradient of kinetic energy defined on edges, for momentum equation @@ -120,19 +133,20 @@ class KEGradOnEdge { bool Enabled; /// constructor declaration - KEGradOnEdge(const HorzMesh *Mesh); + KEGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes edge index, vertical chunk index, and kinetic energy /// array as inputs, outputs the tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, const Array2DReal &KECell) const { - const I4 KStart = KChunk * VecLength; - const I4 JCell0 = CellsOnEdge(IEdge, 0); - const I4 JCell1 = CellsOnEdge(IEdge, 1); + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); const Real InvDcEdge = 1._Real / DcEdge(IEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Tend(IEdge, K) -= EdgeMask(IEdge, K) * (KECell(JCell1, K) - KECell(JCell0, K)) * InvDcEdge; @@ -143,6 +157,8 @@ class KEGradOnEdge { Array2DI4 CellsOnEdge; Array1DReal DcEdge; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Gradient of sea surface height defined on edges multipled by gravitational @@ -152,19 +168,20 @@ class SSHGradOnEdge { bool Enabled; /// constructor declaration - SSHGradOnEdge(const HorzMesh *Mesh); + SSHGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes edge index, vertical chunk index, and array of /// layer thickness/SSH, outputs tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, const Array2DReal &SshCell) const { - const I4 KStart = KChunk * VecLength; - const I4 ICell0 = CellsOnEdge(IEdge, 0); - const I4 ICell1 = CellsOnEdge(IEdge, 1); + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + const I4 ICell0 = CellsOnEdge(IEdge, 0); + const I4 ICell1 = CellsOnEdge(IEdge, 1); const Real InvDcEdge = 1._Real / DcEdge(IEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Tend(IEdge, K) -= EdgeMask(IEdge, K) * Grav * (SshCell(ICell1, K) - SshCell(ICell0, K)) * @@ -177,6 +194,8 @@ class SSHGradOnEdge { Array2DI4 CellsOnEdge; Array1DReal DcEdge; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Laplacian horizontal mixing, for momentum equation @@ -187,7 +206,7 @@ class VelocityDiffusionOnEdge { Real ViscDel2; /// constructor declaration - VelocityDiffusionOnEdge(const HorzMesh *Mesh); + VelocityDiffusionOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes edge index, vertical chunk index, and arrays for /// divergence of horizontal velocity (defined at cell centers) and relative @@ -196,7 +215,8 @@ class VelocityDiffusionOnEdge { const Array2DReal &DivCell, const Array2DReal &RVortVertex) const { - const I4 KStart = KChunk * VecLength; + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); const I4 ICell0 = CellsOnEdge(IEdge, 0); const I4 ICell1 = CellsOnEdge(IEdge, 1); @@ -206,7 +226,7 @@ class VelocityDiffusionOnEdge { const Real DcEdgeInv = 1._Real / DcEdge(IEdge); const Real DvEdgeInv = 1._Real / DvEdge(IEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; const Real Del2U = ((DivCell(ICell1, K) - DivCell(ICell0, K)) * DcEdgeInv - @@ -225,6 +245,8 @@ class VelocityDiffusionOnEdge { Array1DReal DvEdge; Array1DReal MeshScalingDel2; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Biharmonic horizontal mixing, for momentum equation @@ -236,7 +258,7 @@ class VelocityHyperDiffOnEdge { Real DivFactor; /// Constructor declaration - VelocityHyperDiffOnEdge(const HorzMesh *Mesh); + VelocityHyperDiffOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes the edge index, vertical chunk index, and arrays for /// the laplacian of divergence of horizontal velocity and the laplacian of @@ -245,7 +267,8 @@ class VelocityHyperDiffOnEdge { const Array2DReal &Del2DivCell, const Array2DReal &Del2RVortVertex) const { - const I4 KStart = KChunk * VecLength; + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); const I4 ICell0 = CellsOnEdge(IEdge, 0); const I4 ICell1 = CellsOnEdge(IEdge, 1); @@ -255,7 +278,7 @@ class VelocityHyperDiffOnEdge { const Real DcEdgeInv = 1._Real / DcEdge(IEdge); const Real DvEdgeInv = 1._Real / DvEdge(IEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; const Real Del2U = (DivFactor * (Del2DivCell(ICell1, K) - Del2DivCell(ICell0, K)) * @@ -275,6 +298,8 @@ class VelocityHyperDiffOnEdge { Array1DReal DvEdge; Array1DReal MeshScalingDel4; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; /// Wind forcing @@ -284,7 +309,7 @@ class WindForcingOnEdge { Real SaltWaterDensity; /// constructor declaration - WindForcingOnEdge(const HorzMesh *Mesh); + WindForcingOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); /// The functor takes the edge index, vertical chunk index, and arrays for /// normal wind stress and edge layer thickness, outputs tendency array @@ -292,7 +317,7 @@ class WindForcingOnEdge { const Array1DReal &NormalStressEdge, const Array2DReal &LayerThickEdge) const { if (KChunk == 0) { - const I4 K = 0; + const I4 K = MinLayerEdgeBot(IEdge); const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, K); Tend(IEdge, K) += EdgeMask(IEdge, K) * InvThickEdge * @@ -302,6 +327,7 @@ class WindForcingOnEdge { private: Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; }; /// Bottom drag @@ -320,7 +346,7 @@ class BottomDragOnEdge { const Array2DReal &NormalVelEdge, const Array2DReal &KECell, const Array2DReal &LayerThickEdge) const { - const I4 KBot = NVertLayers - 1; + const I4 KBot = MaxLayerEdgeTop(IEdge); const I4 JCell0 = CellsOnEdge(IEdge, 0); const I4 JCell1 = CellsOnEdge(IEdge, 1); @@ -337,6 +363,7 @@ class BottomDragOnEdge { I4 NVertLayers; Array2DI4 CellsOnEdge; Array2DReal EdgeMask; + Array1DI4 MaxLayerEdgeTop; }; // Tracer horizontal advection term @@ -344,30 +371,34 @@ class TracerHorzAdvOnCell { public: bool Enabled; - TracerHorzAdvOnCell(const HorzMesh *Mesh); + TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, const Array2DReal &NormVelEdge, const Array3DReal &HTracersOnEdge) const { - const I4 KStart = KChunk * VecLength; + const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const I4 KEndCell = KStartCell + KLenCell - 1; const Real InvAreaCell = 1._Real / AreaCell(ICell); Real HAdvTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnCell(ICell); ++J) { - const I4 JEdge = EdgesOnCell(ICell, J); + const I4 JEdge = EdgesOnCell(ICell, J); + const I4 KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const I4 KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const I4 KVec = K - KStartCell; HAdvTmp[KVec] -= EdgeMask(JEdge, K) * DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * HTracersOnEdge(L, JEdge, K) * NormVelEdge(JEdge, K) * InvAreaCell; } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const I4 K = KStartCell + KVec; Tend(L, ICell, K) -= HAdvTmp[KVec]; } } @@ -380,6 +411,10 @@ class TracerHorzAdvOnCell { Array1DReal DvEdge; Array1DReal AreaCell; Array2DReal EdgeMask; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; // Tracer horizontal diffusion term @@ -389,20 +424,24 @@ class TracerDiffOnCell { Real EddyDiff2; - TracerDiffOnCell(const HorzMesh *Mesh); + TracerDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, const Array3DReal &TracerCell, const Array2DReal &MeanLayerThickEdge) const { - const I4 KStart = KChunk * VecLength; + const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const I4 KEndCell = KStartCell + KLenCell - 1; const Real InvAreaCell = 1._Real / AreaCell(ICell); Real DiffTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnCell(ICell); ++J) { - const I4 JEdge = EdgesOnCell(ICell, J); + const I4 JEdge = EdgesOnCell(ICell, J); + const I4 KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const I4 KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); const I4 JCell0 = CellsOnEdge(JEdge, 0); const I4 JCell1 = CellsOnEdge(JEdge, 1); @@ -410,8 +449,8 @@ class TracerDiffOnCell { const Real RTemp = MeshScalingDel2(JEdge) * DvEdge(JEdge) / DcEdge(JEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const I4 KVec = K - KStartCell; const Real TracerGrad = (TracerCell(L, JCell1, K) - TracerCell(L, JCell0, K)); @@ -419,8 +458,8 @@ class TracerDiffOnCell { RTemp * MeanLayerThickEdge(JEdge, K) * TracerGrad; } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const I4 K = KStartCell + KVec; Tend(L, ICell, K) += EddyDiff2 * DiffTmp[KVec] * InvAreaCell; } } @@ -435,6 +474,10 @@ class TracerDiffOnCell { Array1DReal AreaCell; Array1DReal MeshScalingDel2; Array2DReal EdgeMask; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; // Tracer biharmonic horizontal mixing term @@ -444,19 +487,23 @@ class TracerHyperDiffOnCell { Real EddyDiff4; - TracerHyperDiffOnCell(const HorzMesh *Mesh); + TracerHyperDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, const Array3DReal &TrDel2Cell) const { - const I4 KStart = KChunk * VecLength; + const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const I4 KEndCell = KStartCell + KLenCell - 1; const Real InvAreaCell = 1._Real / AreaCell(ICell); Real HypTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnCell(ICell); ++J) { - const I4 JEdge = EdgesOnCell(ICell, J); + const I4 JEdge = EdgesOnCell(ICell, J); + const I4 KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const I4 KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); const I4 JCell0 = CellsOnEdge(JEdge, 0); const I4 JCell1 = CellsOnEdge(JEdge, 1); @@ -464,8 +511,8 @@ class TracerHyperDiffOnCell { const Real RTemp = MeshScalingDel4(JEdge) * DvEdge(JEdge) / DcEdge(JEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const I4 KVec = K - KStartCell; const Real Del2TrGrad = (TrDel2Cell(L, JCell1, K) - TrDel2Cell(L, JCell0, K)); @@ -473,8 +520,8 @@ class TracerHyperDiffOnCell { RTemp * Del2TrGrad; } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const I4 K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const I4 K = KStartCell + KVec; Tend(L, ICell, K) -= EddyDiff4 * HypTmp[KVec] * InvAreaCell; } } @@ -489,6 +536,10 @@ class TracerHyperDiffOnCell { Array1DReal AreaCell; Array1DReal MeshScalingDel4; Array2DReal EdgeMask; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 7f3133aa22fd..0bdffb908c79 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -40,7 +40,9 @@ void VertCoord::init2() { Config *OmegaConfig = Config::getOmegaConfig(); - DefaultVertCoord->completeSetup(OmegaConfig); + Halo *MeshHalo = Halo::getDefault(); + + DefaultVertCoord->completeSetup(OmegaConfig, MeshHalo); } // end init2 @@ -128,7 +130,8 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord //------------------------------------------------------------------------------ // Complete construction of new VertCoord instance -void VertCoord::completeSetup(Config *Options //< [in] configuration options +void VertCoord::completeSetup(Config *Options, //< [in] configuration options + Halo *MeshHalo //< [in] mesh halo ) { // Define field metadata @@ -204,6 +207,11 @@ void VertCoord::completeSetup(Config *Options //< [in] configuration options LocMaxLayerCell(ICell) -= 1; }); + // Exchange halos since IOStreams reads only owned cells, but + // we need halos filled for over-computation + MeshHalo->exchangeFullArrayHalo(MinLayerCell, OnCell); + MeshHalo->exchangeFullArrayHalo(MaxLayerCell, OnCell); + // Compute Edge and Vertex vertical ranges minMaxLayerEdge(); minMaxLayerVertex(); diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index 6423c5fa9cf1..7ed13fcd4a4d 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -15,6 +15,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Error.h" +#include "Halo.h" #include "HorzMesh.h" #include "Logging.h" #include "MachEnv.h" @@ -25,6 +26,26 @@ namespace OMEGA { +KOKKOS_INLINE_FUNCTION int vertRange(int KMin, int KMax) { + return KMax - KMin + 1; +} + +KOKKOS_INLINE_FUNCTION int vertRangeChunked(int KMin, int KMax) { + return (vertRange(KMin, KMax) + VecLength - 1) / VecLength; +} + +KOKKOS_INLINE_FUNCTION int chunkStart(int KChunk, int KMin) { + return KMin + KChunk * VecLength; +} + +KOKKOS_INLINE_FUNCTION int chunkLength(int KChunk, int KStart, int KMax) { + if constexpr (VecLength == 1) { + return 1; + } else { + return (KStart + VecLength - 1) > KMax ? (KMax - KStart + 1) : VecLength; + } +} + class VertCoord { private: @@ -152,7 +173,8 @@ class VertCoord { ); /// Read InitialVertCoord stream and complete initialization - void completeSetup(Config *Options /// [in] configuration options + void completeSetup(Config *Options, /// [in] configuration options + Halo *MeshHalo /// [in] mesh halo ); /// Copy member arrays from device to host diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp index 7e40aa2b7992..5e5bb08d4714 100644 --- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp @@ -7,14 +7,15 @@ namespace OMEGA { KineticAuxVars::KineticAuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, int NVertLayers) + const HorzMesh *Mesh, const VertCoord *VCoord) : KineticEnergyCell("KineticEnergyCell" + AuxStateSuffix, Mesh->NCellsSize, - NVertLayers), + VCoord->NVertLayers), VelocityDivCell("VelocityDivCell" + AuxStateSuffix, Mesh->NCellsSize, - NVertLayers), + VCoord->NVertLayers), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), DcEdge(Mesh->DcEdge), - DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {} + DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} void KineticAuxVars::registerFields( const std::string &AuxGroupName, // name of Auxiliary field group diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h index 3d23acceafc7..850807702f64 100644 --- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h @@ -4,6 +4,7 @@ #include "DataTypes.h" #include "HorzMesh.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include @@ -15,13 +16,16 @@ class KineticAuxVars { Array2DReal VelocityDivCell; KineticAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - int NVertLayers); + const VertCoord *VCoord); KOKKOS_FUNCTION void computeVarsOnCell(int ICell, int KChunk, const Array2DReal &NormalVelEdge) const { + + const int KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + const Real InvAreaCell = 1._Real / AreaCell(ICell); - const int KStart = KChunk * VecLength; Real KineticEnergyCellTmp[VecLength] = {0}; Real VelocityDivCellTmp[VecLength] = {0}; @@ -29,7 +33,7 @@ class KineticAuxVars { for (int J = 0; J < NEdgesOnCell(ICell); ++J) { const int JEdge = EdgesOnCell(ICell, J); const Real AreaEdge = 0.5_Real * DvEdge(JEdge) * DcEdge(JEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; KineticEnergyCellTmp[KVec] += AreaEdge * 0.5_Real * InvAreaCell * NormalVelEdge(JEdge, K) * @@ -39,7 +43,7 @@ class KineticAuxVars { NormalVelEdge(JEdge, K); } } - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; KineticEnergyCell(ICell, K) = KineticEnergyCellTmp[KVec]; VelocityDivCell(ICell, K) = VelocityDivCellTmp[KVec]; @@ -57,6 +61,8 @@ class KineticAuxVars { Array1DReal DcEdge; Array1DReal DvEdge; Array1DReal AreaCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index 6ed575df788a..cb937863528c 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -7,13 +7,17 @@ namespace OMEGA { LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - int NVertLayers) + const VertCoord *VCoord) : FluxLayerThickEdge("FluxLayerThickEdge" + AuxStateSuffix, - Mesh->NEdgesSize, NVertLayers), + Mesh->NEdgesSize, VCoord->NVertLayers), MeanLayerThickEdge("MeanLayerThickEdge" + AuxStateSuffix, - Mesh->NEdgesSize, NVertLayers), - SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize, NVertLayers), - CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth) {} + Mesh->NEdgesSize, VCoord->NVertLayers), + SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize, + VCoord->NVertLayers), + CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h index 44562aa6cc5c..0620d546fd8f 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h @@ -4,6 +4,7 @@ #include "DataTypes.h" #include "HorzMesh.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include @@ -20,16 +21,18 @@ class LayerThicknessAuxVars { FluxThickEdgeOption FluxThickEdgeChoice; LayerThicknessAuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, int NVertLayers); + const HorzMesh *Mesh, const VertCoord *VCoord); KOKKOS_FUNCTION void computeVarsOnEdge(int IEdge, int KChunk, const Array2DReal &LayerThickCell, const Array2DReal &NormalVelEdge) const { - const int KStart = KChunk * VecLength; + const int KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + const int JCell0 = CellsOnEdge(IEdge, 0); const int JCell1 = CellsOnEdge(IEdge, 1); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; MeanLayerThickEdge(IEdge, K) = 0.5_Real * (LayerThickCell(JCell0, K) + LayerThickCell(JCell1, K)); @@ -37,7 +40,7 @@ class LayerThicknessAuxVars { switch (FluxThickEdgeChoice) { case FluxThickEdgeOption::Center: - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; FluxLayerThickEdge(IEdge, K) = 0.5_Real * @@ -45,7 +48,7 @@ class LayerThicknessAuxVars { } break; case FluxThickEdgeOption::Upwind: - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; if (NormalVelEdge(IEdge, K) > 0) { FluxLayerThickEdge(IEdge, K) = LayerThickCell(JCell0, K); @@ -65,8 +68,10 @@ class LayerThicknessAuxVars { const Array2DReal &LayerThickCell) const { // Temporary for stacked shallow water - const int KStart = KChunk * VecLength; - for (int KVec = 0; KVec < VecLength; ++KVec) { + const int KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell); } @@ -88,6 +93,10 @@ class LayerThicknessAuxVars { private: Array2DI4 CellsOnEdge; Array1DReal BottomDepth; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp index 1ef7eefa6394..385cfea28401 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp @@ -6,16 +6,18 @@ namespace OMEGA { TracerAuxVars::TracerAuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, const I4 NVertLayers, + const HorzMesh *Mesh, const VertCoord *VCoord, const I4 NTracers) : HTracersEdge("ThickTracersEdge" + AuxStateSuffix, NTracers, - Mesh->NEdgesSize, NVertLayers), + Mesh->NEdgesSize, VCoord->NVertLayers), Del2TracersCell("Del2TracerCell" + AuxStateSuffix, NTracers, - Mesh->NCellsSize, NVertLayers), + Mesh->NCellsSize, VCoord->NVertLayers), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), - EdgeMask(Mesh->EdgeMask) {} + EdgeMask(Mesh->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} void TracerAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h index 39e6f5bf232a..1d12e55d9585 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h @@ -5,6 +5,7 @@ #include "Field.h" #include "HorzMesh.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include @@ -20,19 +21,21 @@ class TracerAuxVars { FluxTracerEdgeOption TracersOnEdgeChoice; TracerAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - const I4 NVertLayers, const I4 NTracers); + const VertCoord *VCoord, const I4 NTracers); KOKKOS_FUNCTION void computeVarsOnEdge(int L, int IEdge, int KChunk, const Array2DReal &NormalVelEdge, const Array2DReal &HCell, const Array3DReal &TrCell) const { - const int KStart = KChunk * VecLength; + const int KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + const int JCell0 = CellsOnEdge(IEdge, 0); const int JCell1 = CellsOnEdge(IEdge, 1); switch (TracersOnEdgeChoice) { case FluxTracerEdgeOption::Center: - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; HTracersEdge(L, IEdge, K) = 0.5_Real * (HCell(JCell0, K) * TrCell(L, JCell0, K) + @@ -40,7 +43,7 @@ class TracerAuxVars { } break; case FluxTracerEdgeOption::Upwind: - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; if (NormalVelEdge(IEdge, K) > 0) { HTracersEdge(L, IEdge, K) = @@ -63,7 +66,10 @@ class TracerAuxVars { const Array2DReal &LayerThickEdgeMean, const Array3DReal &TrCell) const { - const int KStart = KChunk * VecLength; + const int KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const int KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const int KEndCell = KStartCell + KLenCell - 1; + const Real InvAreaCell = 1._Real / AreaCell(ICell); Real Del2TrCellTmp[VecLength] = {0}; @@ -76,16 +82,19 @@ class TracerAuxVars { const Real DvDcEdge = DvEdge(JEdge) / DcEdge(JEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + const int KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const int KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); + + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const int KVec = K - KStartCell; const Real TracerGrad = TrCell(L, JCell1, K) - TrCell(L, JCell0, K); Del2TrCellTmp[KVec] -= EdgeMask(JEdge, K) * EdgeSignOnCell(ICell, J) * DvDcEdge * LayerThickEdgeMean(JEdge, K) * TracerGrad; } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const int K = KStartCell + KVec; Del2TracersCell(L, ICell, K) = Del2TrCellTmp[KVec] * InvAreaCell; } } @@ -103,6 +112,10 @@ class TracerAuxVars { Array1DReal DvEdge; Array1DReal AreaCell; Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp index 8140723b83f9..85b09b405b9d 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp @@ -7,19 +7,28 @@ namespace OMEGA { VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, int NVertLayers) - : Del2Edge("VelDel2Edge" + AuxStateSuffix, Mesh->NEdgesSize, NVertLayers), + const HorzMesh *Mesh, + const VertCoord *VCoord) + : Del2Edge("VelDel2Edge" + AuxStateSuffix, Mesh->NEdgesSize, + VCoord->NVertLayers), Del2DivCell("VelDel2DivCell" + AuxStateSuffix, Mesh->NCellsSize, - NVertLayers), + VCoord->NVertLayers), Del2RelVortVertex("VelDel2RelVortVertex" + AuxStateSuffix, - Mesh->NVerticesSize, NVertLayers), + Mesh->NVerticesSize, VCoord->NVertLayers), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), EdgesOnVertex(Mesh->EdgesOnVertex), CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge), EdgeMask(Mesh->EdgeMask), EdgeSignOnVertex(Mesh->EdgeSignOnVertex), - AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree) {} + AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), + MinLayerVertexBot(VCoord->MinLayerVertexBot), + MaxLayerVertexTop(VCoord->MaxLayerVertexTop), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) + +{} void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h index 468d5ea9279e..b6f05f9c0d8b 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h @@ -4,6 +4,7 @@ #include "DataTypes.h" #include "HorzMesh.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include @@ -16,12 +17,13 @@ class VelocityDel2AuxVars { Array2DReal Del2RelVortVertex; VelocityDel2AuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - int NVertLayers); + const VertCoord *VCoord); KOKKOS_FUNCTION void computeVarsOnEdge(int IEdge, int KChunk, const Array2DReal &VelocityDivCell, const Array2DReal &RelVortVertex) const { - const int KStart = KChunk * VecLength; + const int KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); const int JCell0 = CellsOnEdge(IEdge, 0); const int JCell1 = CellsOnEdge(IEdge, 1); @@ -32,7 +34,7 @@ class VelocityDel2AuxVars { const Real InvDvEdge = 1._Real / Kokkos::max(DvEdge(IEdge), 0.25_Real * DcEdge(IEdge)); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; const Real GradDiv = (VelocityDivCell(JCell1, K) - VelocityDivCell(JCell0, K)) * @@ -46,35 +48,43 @@ class VelocityDel2AuxVars { KOKKOS_FUNCTION void computeVarsOnCell(int ICell, int KChunk) const { const Real InvAreaCell = 1._Real / AreaCell(ICell); - const int KStart = KChunk * VecLength; + const int KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const int KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const int KEndCell = KStartCell + KLenCell - 1; Real Del2DivCellTmp[VecLength] = {0}; for (int J = 0; J < NEdgesOnCell(ICell); ++J) { const int JEdge = EdgesOnCell(ICell, J); const Real AreaEdge = 0.5_Real * DvEdge(JEdge) * DcEdge(JEdge); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + + const int KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); + const int KEndEdge = Kokkos::min(KEndCell, MaxLayerEdgeTop(JEdge)); + + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const int KVec = K - KStartCell; Del2DivCellTmp[KVec] -= DvEdge(JEdge) * InvAreaCell * EdgeSignOnCell(ICell, J) * Del2Edge(JEdge, K); } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + for (int KVec = 0; KVec < KLenCell; ++KVec) { + const int K = KStartCell + KVec; Del2DivCell(ICell, K) = Del2DivCellTmp[KVec]; } } KOKKOS_FUNCTION void computeVarsOnVertex(int IVertex, int KChunk) const { - const int KStart = KChunk * VecLength; + const int KStart = chunkStart(KChunk, MinLayerVertexBot(IVertex)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerVertexTop(IVertex)); + const Real InvAreaTriangle = 1._Real / AreaTriangle(IVertex); Real Del2RelVortVertexTmp[VecLength] = {0}; for (int J = 0; J < VertexDegree; ++J) { const int JEdge = EdgesOnVertex(IVertex, J); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; Del2RelVortVertexTmp[KVec] += InvAreaTriangle * DcEdge(JEdge) * EdgeSignOnVertex(IVertex, J) * @@ -82,7 +92,7 @@ class VelocityDel2AuxVars { } } - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; Del2RelVortVertex(IVertex, K) = Del2RelVortVertexTmp[KVec]; } @@ -106,6 +116,12 @@ class VelocityDel2AuxVars { Array1DReal AreaTriangle; Array2DReal EdgeMask; I4 VertexDegree; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + Array1DI4 MinLayerVertexBot; + Array1DI4 MaxLayerVertexTop; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp index 53fc77860d3e..4d20cac1f838 100644 --- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp @@ -7,23 +7,29 @@ namespace OMEGA { VorticityAuxVars::VorticityAuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, int NVertLayers) + const HorzMesh *Mesh, + const VertCoord *VCoord) : RelVortVertex("RelVortVertex" + AuxStateSuffix, Mesh->NVerticesSize, - NVertLayers), + VCoord->NVertLayers), NormRelVortVertex("NormRelVortVertex" + AuxStateSuffix, - Mesh->NVerticesSize, NVertLayers), + Mesh->NVerticesSize, VCoord->NVertLayers), NormPlanetVortVertex("NormPlanetVortVertex" + AuxStateSuffix, - Mesh->NVerticesSize, NVertLayers), + Mesh->NVerticesSize, VCoord->NVertLayers), NormRelVortEdge("NormRelVortEdge" + AuxStateSuffix, Mesh->NEdgesSize, - NVertLayers), + VCoord->NVertLayers), NormPlanetVortEdge("NormPlanetVortEdge" + AuxStateSuffix, - Mesh->NEdgesSize, NVertLayers), + Mesh->NEdgesSize, VCoord->NVertLayers), VertexDegree(Mesh->VertexDegree), CellsOnVertex(Mesh->CellsOnVertex), EdgesOnVertex(Mesh->EdgesOnVertex), EdgeSignOnVertex(Mesh->EdgeSignOnVertex), DcEdge(Mesh->DcEdge), KiteAreasOnVertex(Mesh->KiteAreasOnVertex), AreaTriangle(Mesh->AreaTriangle), FVertex(Mesh->FVertex), - VerticesOnEdge(Mesh->VerticesOnEdge) {} + VerticesOnEdge(Mesh->VerticesOnEdge), + MinLayerVertexTop(VCoord->MinLayerVertexTop), + MaxLayerVertexBot(VCoord->MaxLayerVertexBot), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeTop(VCoord->MinLayerEdgeTop), + MaxLayerEdgeBot(VCoord->MaxLayerEdgeBot) {} void VorticityAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h index 64feb94f81cd..eecfcf5fd993 100644 --- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h @@ -4,6 +4,7 @@ #include "DataTypes.h" #include "HorzMesh.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include @@ -19,13 +20,18 @@ class VorticityAuxVars { Array2DReal NormPlanetVortEdge; VorticityAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - int NVertLayers); + const VertCoord *VCoord); KOKKOS_FUNCTION void computeVarsOnVertex(int IVertex, int KChunk, const Array2DReal &LayerThickCell, const Array2DReal &NormalVelEdge) const { - const int KStart = KChunk * VecLength; + + const int KStartVertex = chunkStart(KChunk, MinLayerVertexTop(IVertex)); + const int KLenVertex = + chunkLength(KChunk, KStartVertex, MaxLayerVertexBot(IVertex)); + const int KEndVertex = KStartVertex + KLenVertex - 1; + const Real InvAreaTriangle = 1._Real / AreaTriangle(IVertex); Real LayerThickVertex[VecLength] = {0}; @@ -33,21 +39,32 @@ class VorticityAuxVars { for (int J = 0; J < VertexDegree; ++J) { const int JCell = CellsOnVertex(IVertex, J); - const int JEdge = EdgesOnVertex(IVertex, J); - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + const int KStartCell = Kokkos::max(KStartVertex, MinLayerCell(JCell)); + const int KEndCell = Kokkos::min(KEndVertex, MaxLayerCell(JCell)); + + for (int K = KStartCell; K <= KEndCell; ++K) { + const int KVec = K - KStartVertex; LayerThickVertex[KVec] += InvAreaTriangle * KiteAreasOnVertex(IVertex, J) * LayerThickCell(JCell, K); + } + + const int JEdge = EdgesOnVertex(IVertex, J); + const int KStartEdge = + Kokkos::max(KStartVertex, MinLayerEdgeTop(JEdge)); + const int KEndEdge = Kokkos::min(KEndVertex, MaxLayerEdgeBot(JEdge)); + + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const int KVec = K - KStartVertex; RelVortVertexTmp[KVec] += InvAreaTriangle * DcEdge(JEdge) * EdgeSignOnVertex(IVertex, J) * NormalVelEdge(JEdge, K); } } - for (int KVec = 0; KVec < VecLength; ++KVec) { - const int K = KStart + KVec; + for (int KVec = 0; KVec < KLenVertex; ++KVec) { + const int K = KStartVertex + KVec; const Real InvLayerThickVertex = 1._Real / LayerThickVertex[KVec]; RelVortVertex(IVertex, K) = RelVortVertexTmp[KVec]; @@ -59,11 +76,12 @@ class VorticityAuxVars { } KOKKOS_FUNCTION void computeVarsOnEdge(int IEdge, int KChunk) const { - const int KStart = KChunk * VecLength; + const int KStart = chunkStart(KChunk, MinLayerEdgeTop(IEdge)); + const int KLen = chunkLength(KChunk, KStart, MaxLayerEdgeBot(IEdge)); const int JVertex0 = VerticesOnEdge(IEdge, 0); const int JVertex1 = VerticesOnEdge(IEdge, 1); - for (int KVec = 0; KVec < VecLength; ++KVec) { + for (int KVec = 0; KVec < KLen; ++KVec) { const int K = KStart + KVec; NormRelVortEdge(IEdge, K) = 0.5_Real * @@ -89,6 +107,13 @@ class VorticityAuxVars { Array1DReal AreaTriangle; Array2DI4 VerticesOnEdge; Array1DReal FVertex; + + Array1DI4 MinLayerVertexTop; + Array1DI4 MaxLayerVertexBot; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeTop; + Array1DI4 MaxLayerEdgeBot; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h index b1546f441216..11d33816d32f 100644 --- a/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/WindForcingAuxVars.h @@ -5,6 +5,7 @@ #include "HorzMesh.h" #include "HorzOperators.h" #include "OmegaKokkos.h" +#include "VertCoord.h" #include diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index 217c21de46df..c95b8f3f70d3 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -7,6 +7,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -50,8 +51,9 @@ int initState() { int Err = 0; TestSetup Setup; - auto *Mesh = HorzMesh::getDefault(); - auto *State = OceanState::getDefault(); + auto *Mesh = HorzMesh::getDefault(); + auto *State = OceanState::getDefault(); + auto *VCoord = VertCoord::getDefault(); Array2DReal LayerThickCell; Array2DReal NormalVelEdge; @@ -64,18 +66,21 @@ int initState() { Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickCell, Geom, Mesh, OnCell); + LayerThickCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, + VCoord->MaxLayerCell); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracerArray, Geom, Mesh, OnCell); + TracerArray, Geom, Mesh, OnCell, VCoord->MinLayerCell, + VCoord->MaxLayerCell); 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, ExchangeHalos::Yes, + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, ExchangeHalos::Yes, CartProjection::No); return Err; @@ -102,6 +107,9 @@ int initAuxStateTest(const std::string &mesh) { IO::init(DefComm); Decomp::init(mesh); + // Initialize streams + IOStream::init(); + int HaloErr = Halo::init(); if (HaloErr != 0) { Err++; @@ -109,8 +117,11 @@ int initAuxStateTest(const std::string &mesh) { } VertCoord::init1(); + HorzMesh::init(); + VertCoord::init2(); + Tracers::init(); int StateErr = OceanState::init(); if (StateErr != 0) { @@ -138,10 +149,11 @@ int testAuxState() { return -1; } - const auto *Mesh = HorzMesh::getDefault(); - auto *MeshHalo = Halo::getDefault(); + const auto *Mesh = HorzMesh::getDefault(); + auto *MeshHalo = Halo::getDefault(); + const auto *VCoord = VertCoord::getDefault(); // test creation of another auxiliary state - AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, 12, 3); + AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, VCoord, 3); // test retrievel of another if (AuxiliaryState::get("AnotherAuxState")) { @@ -194,77 +206,112 @@ int testAuxState() { int NTracers = Tracers::getNumTracers(); const Real KineticEnergySum = - sum(DefAuxState->KineticAux.KineticEnergyCell, NCellsOwned); + sum(DefAuxState->KineticAux.KineticEnergyCell, NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); if (!Kokkos::isfinite(KineticEnergySum)) { Err++; LOG_ERROR("AuxStateTest: KineticEnergy FAIL"); } const Real VelocityDivSum = - sum(DefAuxState->KineticAux.VelocityDivCell, NCellsOwned); + sum(DefAuxState->KineticAux.VelocityDivCell, NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); if (!Kokkos::isfinite(VelocityDivSum)) { Err++; LOG_ERROR("AuxStateTest: VelocityDivCell FAIL"); } const Real FluxLayerThickSum = - sum(DefAuxState->LayerThicknessAux.FluxLayerThickEdge, NEdgesOwned); + sum(DefAuxState->LayerThicknessAux.FluxLayerThickEdge, NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); if (!Kokkos::isfinite(FluxLayerThickSum)) { Err++; LOG_ERROR("AuxStateTest: FluxLayerThickEdge FAIL"); } const Real MeanLayerThickSum = - sum(DefAuxState->LayerThicknessAux.MeanLayerThickEdge, NEdgesOwned); + sum(DefAuxState->LayerThicknessAux.MeanLayerThickEdge, NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); if (!Kokkos::isfinite(MeanLayerThickSum)) { Err++; LOG_ERROR("AuxStateTest: MeanLayerThickEdge FAIL"); } const Real RelVortVSum = - sum(DefAuxState->VorticityAux.RelVortVertex, NVerticesOwned); + sum(DefAuxState->VorticityAux.RelVortVertex, NVerticesOwned, + VCoord->MinLayerVertexTop, VCoord->MaxLayerVertexBot); if (!Kokkos::isfinite(RelVortVSum)) { Err++; LOG_ERROR("AuxStateTest: RelVortVertex FAIL"); } const Real NormRelVortVSum = - sum(DefAuxState->VorticityAux.NormRelVortVertex, NVerticesOwned); + sum(DefAuxState->VorticityAux.NormRelVortVertex, NVerticesOwned, + VCoord->MinLayerVertexTop, VCoord->MaxLayerVertexBot); if (!Kokkos::isfinite(NormRelVortVSum)) { Err++; LOG_ERROR("AuxStateTest: NormRelVortVertex FAIL"); } const Real NormPlanetVortVSum = - sum(DefAuxState->VorticityAux.NormPlanetVortVertex, NVerticesOwned); + sum(DefAuxState->VorticityAux.NormPlanetVortVertex, NVerticesOwned, + VCoord->MinLayerVertexTop, VCoord->MaxLayerVertexBot); if (!Kokkos::isfinite(NormPlanetVortVSum)) { Err++; LOG_ERROR("AuxStateTest: NormPlanetVortVertex FAIL"); } const Real NormRelVortESum = - sum(DefAuxState->VorticityAux.NormRelVortEdge, NEdgesOwned); + sum(DefAuxState->VorticityAux.NormRelVortEdge, NEdgesOwned, + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot); if (!Kokkos::isfinite(NormRelVortESum)) { Err++; LOG_ERROR("AuxStateTest: NormRelVortEdge FAIL"); } const Real NormPlanetVortESum = - sum(DefAuxState->VorticityAux.NormPlanetVortEdge, NEdgesOwned); + sum(DefAuxState->VorticityAux.NormPlanetVortEdge, NEdgesOwned, + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot); if (!Kokkos::isfinite(NormPlanetVortESum)) { Err++; LOG_ERROR("AuxStateTest: NormPlanetVortEdge FAIL"); } + const Real Del2EdgeSum = + sum(DefAuxState->VelocityDel2Aux.Del2Edge, NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); + if (!Kokkos::isfinite(Del2EdgeSum)) { + Err++; + LOG_ERROR("AuxStateTest: Del2Edge FAIL"); + } + + const Real Del2DivCellSum = + sum(DefAuxState->VelocityDel2Aux.Del2DivCell, NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + if (!Kokkos::isfinite(Del2DivCellSum)) { + Err++; + LOG_ERROR("AuxStateTest: Del2DivCell FAIL"); + } + + const Real Del2RelVortVertexSum = + sum(DefAuxState->VelocityDel2Aux.Del2RelVortVertex, NVerticesOwned, + VCoord->MinLayerVertexBot, VCoord->MaxLayerVertexTop); + if (!Kokkos::isfinite(Del2RelVortVertexSum)) { + Err++; + LOG_ERROR("AuxStateTest: Del2RelVortVertex FAIL"); + } + const Real HTracersESum = - sum(DefAuxState->TracerAux.HTracersEdge, NTracers, NEdgesOwned); + sum(DefAuxState->TracerAux.HTracersEdge, NTracers, NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); if (!Kokkos::isfinite(HTracersESum)) { Err++; LOG_ERROR("AuxStateTest: HTracersOnEdge FAIL"); } const Real Del2TracersCSum = - sum(DefAuxState->TracerAux.Del2TracersCell, NTracers, NCellsOwned); + sum(DefAuxState->TracerAux.Del2TracersCell, NTracers, NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); if (!Kokkos::isfinite(Del2TracersCSum)) { Err++; LOG_ERROR("AuxStateTest: Del2TracersOnCell FAIL"); diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 00505498439c..264a80fca5bc 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -345,7 +345,8 @@ int testKineticAuxVars(const Array2DReal &LayerThicknessCell, int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result @@ -363,7 +364,7 @@ int testKineticAuxVars(const Array2DReal &LayerThicknessCell, // Compute numerical result - KineticAuxVars KineticAux("", Mesh, NVertLayers); + KineticAuxVars KineticAux("", Mesh, VCoord); parallelFor( {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int KLayer) { @@ -449,7 +450,8 @@ int testLayerThicknessAuxVars(const Array2DReal &LayerThickCell, int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result @@ -460,7 +462,7 @@ int testLayerThicknessAuxVars(const Array2DReal &LayerThickCell, // Compute numerical result - LayerThicknessAuxVars LayerThicknessAux("", Mesh, NVertLayers); + LayerThicknessAuxVars LayerThicknessAux("", Mesh, VCoord); LayerThicknessAux.FluxThickEdgeChoice = FluxThickEdgeOption::Upwind; parallelFor( {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) { @@ -499,7 +501,8 @@ int testVorticityAuxVars(const Array2DReal &LayerThickCell, const auto Decomp = Decomp::getDefault(); const auto Mesh = HorzMesh::getDefault(); - VorticityAuxVars VorticityAux("", Mesh, NVertLayers); + const auto VCoord = VertCoord::getDefault(); + VorticityAuxVars VorticityAux("", Mesh, VCoord); // Compute exact results for vertex variables @@ -615,7 +618,8 @@ int testVelocityDel2AuxVars(Real RTol) { const auto Decomp = Decomp::getDefault(); const auto Mesh = HorzMesh::getDefault(); - VelocityDel2AuxVars VelocityDel2Aux("", Mesh, NVertLayers); + const auto VCoord = VertCoord::getDefault(); + VelocityDel2AuxVars VelocityDel2Aux("", Mesh, VCoord); // Use analytical expressions to compute inputs @@ -721,9 +725,10 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, TestSetup Setup; int Err = 0; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); - TracerAuxVars TracerAux("", Mesh, NVertLayers, NTracers); + TracerAuxVars TracerAux("", Mesh, VCoord, NTracers); TracerAux.TracersOnEdgeChoice = FluxTracerEdgeOption::Upwind; // Set input arrays @@ -825,11 +830,7 @@ int initAuxVarsTest(const std::string &mesh) { VertCoord::init1(); // Reset NVertLayers to the test value - auto *DefVertCoord = VertCoord::getDefault(); - DefVertCoord->NVertLayers = NVertLayers; - Dimension::destroy("NVertLayers"); - std::shared_ptr VertDim = - Dimension::create("NVertLayers", NVertLayers); + resetVertCoord(VertCoord::getDefault(), NVertLayers); HorzMesh::init(); diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 028d136e9c16..08535f59cc7e 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -98,8 +98,9 @@ I4 initEosTest(const std::string &mesh) { /// Test Linear EOS calculation for all cells/layers int testEosLinear() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto *Mesh = HorzMesh::getDefault(); + const auto *VCoord = VertCoord::getDefault(); /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::LinearEos; @@ -117,24 +118,40 @@ int testEosLinear() { /// Compute specific volume TestEos->computeSpecVol(TArray, SArray, PArray); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + /// Check all array values against expected value - int numMismatches = 0; + int NumMismatches = 0; Array2DReal SpecVol = TestEos->SpecVol; - parallelReduce( - "CheckSpecVolMatrix-linear", {Mesh->NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVol(i, j), LinearExpValue, RTol)) { - localCount++; - } + parallelReduceOuter( + "CheckSpecVolMatrix-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVol(ICell, K), LinearExpValue, RTol)) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); }, - numMismatches); + NumMismatches); auto SpecVolH = createHostMirrorCopy(SpecVol); - if (numMismatches != 0) { + if (NumMismatches != 0) { Err++; LOG_ERROR("EosTest: SpecVol Linear isApprox FAIL, " "expected {}, got {} with {} mismatches", - LinearExpValue, SpecVolH(1, 1), numMismatches); + LinearExpValue, SpecVolH(1, 1), NumMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalc Linear: PASS"); @@ -145,8 +162,9 @@ int testEosLinear() { /// Test Linear EOS calculation with vertical displacement int testEosLinearDisplaced() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto *Mesh = HorzMesh::getDefault(); + const auto *VCoord = VertCoord::getDefault(); /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::LinearEos; @@ -164,24 +182,41 @@ int testEosLinearDisplaced() { /// Compute displaced specific volume TestEos->computeSpecVolDisp(TArray, SArray, PArray, KDisp); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + /// Check all array values against expected value - int numMismatches = 0; + int NumMismatches = 0; Array2DReal SpecVolDisplaced = TestEos->SpecVolDisplaced; - parallelReduce( - "CheckSpecVolDispMatrix-linear", {Mesh->NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVolDisplaced(i, j), LinearExpValue, RTol)) { - localCount++; - } + parallelReduceOuter( + "CheckSpecVolDispMatrix-linear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVolDisplaced(ICell, K), LinearExpValue, + RTol)) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); }, - numMismatches); + NumMismatches); auto SpecVolDisplacedH = createHostMirrorCopy(SpecVolDisplaced); - if (numMismatches != 0) { + if (NumMismatches != 0) { Err++; LOG_ERROR("EosTest: Linear SpecVolDisp isApprox FAIL, " "expected {}, got {} with {} mismatches", - LinearExpValue, SpecVolDisplacedH(1, 1), numMismatches); + LinearExpValue, SpecVolDisplacedH(1, 1), NumMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalcDisp Linear: PASS"); @@ -192,8 +227,9 @@ int testEosLinearDisplaced() { /// Test TEOS-10 EOS calculation for all cells/layers int testEosTeos10() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto *Mesh = HorzMesh::getDefault(); + const auto *VCoord = VertCoord::getDefault(); /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; @@ -211,24 +247,40 @@ int testEosTeos10() { /// Compute specific volume TestEos->computeSpecVol(TArray, SArray, PArray); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + /// Check all array values against expected value - int numMismatches = 0; + int NumMismatches = 0; Array2DReal SpecVol = TestEos->SpecVol; - parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVol(i, j), TeosExpValue, RTol)) { - localCount++; - } + parallelReduceOuter( + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVol(ICell, K), TeosExpValue, RTol)) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); }, - numMismatches); + NumMismatches); auto SpecVolH = createHostMirrorCopy(SpecVol); - if (numMismatches != 0) { + if (NumMismatches != 0) { Err++; LOG_ERROR("EosTest: TEOS SpecVol isApprox FAIL, " "expected {}, got {} with {} mismatches", - TeosExpValue, SpecVolH(1, 1), numMismatches); + TeosExpValue, SpecVolH(1, 1), NumMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalc TEOS-10: PASS"); @@ -239,8 +291,9 @@ int testEosTeos10() { /// Test TEOS-10 EOS calculation with vertical displacement int testEosTeos10Displaced() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto *Mesh = HorzMesh::getDefault(); + const auto *VCoord = VertCoord::getDefault(); /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; @@ -258,24 +311,41 @@ int testEosTeos10Displaced() { /// Compute displaced specific volume TestEos->computeSpecVolDisp(TArray, SArray, PArray, KDisp); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + /// Check all array values against expected value - int numMismatches = 0; + int NumMismatches = 0; Array2DReal SpecVolDisplaced = TestEos->SpecVolDisplaced; - parallelReduce( - "CheckSpecVolDispMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVolDisplaced(i, j), TeosExpValue, RTol)) { - localCount++; - } + parallelReduceOuter( + "CheckSpecVolDispMatrix-Teos", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(SpecVolDisplaced(ICell, K), TeosExpValue, + RTol)) { + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); }, - numMismatches); + NumMismatches); auto SpecVolDisplacedH = createHostMirrorCopy(SpecVolDisplaced); - if (numMismatches != 0) { + if (NumMismatches != 0) { Err++; LOG_ERROR("EosTest: TEOS SpecVolDisp isApprox FAIL, " "expected {}, got {} with {} mismatches", - TeosExpValue, SpecVolDisplacedH(1, 1), numMismatches); + TeosExpValue, SpecVolDisplacedH(1, 1), NumMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalcDisp TEOS-10: PASS"); diff --git a/components/omega/test/ocn/OceanTestCommon.h b/components/omega/test/ocn/OceanTestCommon.h index 337119f73fbc..a73fff1de369 100644 --- a/components/omega/test/ocn/OceanTestCommon.h +++ b/components/omega/test/ocn/OceanTestCommon.h @@ -2,11 +2,13 @@ #define OMEGA_OCEAN_TEST_COMMON_H #include "DataTypes.h" +#include "Dimension.h" #include "Halo.h" #include "HorzMesh.h" #include "Logging.h" #include "MachEnv.h" #include "OmegaKokkos.h" +#include "VertCoord.h" namespace OMEGA { @@ -66,23 +68,58 @@ KOKKOS_INLINE_FUNCTION void tangentVector(Real (&TanVec)[3], enum class EdgeComponent { Normal, Tangential }; enum class Geometry { Planar, Spherical }; enum class ExchangeHalos { Yes, No }; +enum class SetBoundary { Yes = 1, No = 0 }; + +// change the number of layers in a vertical coordinate to a constant value +void resetVertCoord(VertCoord *VCoord, int NVertLayers) { + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayers + 1; + + deepCopy(VCoord->MinLayerCell, 0); + deepCopy(VCoord->MaxLayerCell, NVertLayers - 1); + + VCoord->minMaxLayerEdge(); + VCoord->minMaxLayerVertex(); + + Dimension::destroy("NVertLayers"); + Dimension::destroy("NVertLayersP1"); + std::shared_ptr VertDim = + Dimension::create("NVertLayers", NVertLayers); + std::shared_ptr VertDimP1 = + Dimension::create("NVertLayersP1", NVertLayers + 1); +} + +// helper to get vertical iteration bounds that can be provided +// either as an integer or an array +template KOKKOS_FUNCTION int getVertBound(const T &VertBound, int I) { + if constexpr (std::is_integral_v) { + return VertBound; + } else { + static_assert(Kokkos::is_view_v); + return VertBound(I); + } +} // set scalar field on chosen elements (cells/vertices/edges) based on // analytical formula and optionally exchange halos -template +template int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, - const HorzMesh *Mesh, MeshElement Element, - ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes) { + const HorzMesh *Mesh, MeshElement Element, const VertMin &VMin, + const VertMax &VMax, + ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, + SetBoundary SetBndOpt = SetBoundary::No) { int Err = 0; int NElementsOwned; + int NElementsSize; Array1DReal XElement, YElement; Array1DReal LonElement, LatElement; switch (Element) { case OnCell: NElementsOwned = Mesh->NCellsOwned; + NElementsSize = Mesh->NCellsSize; XElement = createDeviceMirrorCopy(Mesh->XCellH); YElement = createDeviceMirrorCopy(Mesh->YCellH); LonElement = createDeviceMirrorCopy(Mesh->LonCellH); @@ -90,6 +127,7 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, break; case OnVertex: NElementsOwned = Mesh->NVerticesOwned; + NElementsSize = Mesh->NVerticesSize; XElement = createDeviceMirrorCopy(Mesh->XVertexH); YElement = createDeviceMirrorCopy(Mesh->YVertexH); LonElement = createDeviceMirrorCopy(Mesh->LonVertexH); @@ -97,6 +135,7 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, break; case OnEdge: NElementsOwned = Mesh->NEdgesOwned; + NElementsSize = Mesh->NEdgesSize; XElement = createDeviceMirrorCopy(Mesh->XEdgeH); YElement = createDeviceMirrorCopy(Mesh->YEdgeH); LonElement = createDeviceMirrorCopy(Mesh->LonEdgeH); @@ -108,9 +147,14 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, return 1; } + const int NElementsSet = NElementsOwned + static_cast(SetBndOpt); + if constexpr (Array::rank == 1) { parallelFor( - {NElementsOwned}, KOKKOS_LAMBDA(int IElement) { + {NElementsSet}, KOKKOS_LAMBDA(int IElement) { + if (SetBndOpt == SetBoundary::Yes) { + IElement = IElement == 0 ? (NElementsSize - 1) : (IElement - 1); + } if (Geom == Geometry::Planar) { const Real X = XElement(IElement); const Real Y = YElement(IElement); @@ -126,35 +170,54 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, if constexpr (Array::rank == 2) { const int NVertLayers = ScalarElement.extent_int(1); - parallelFor( - {NElementsOwned, NVertLayers}, 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); + parallelForOuter( + {NElementsSet}, KOKKOS_LAMBDA(int IElement, const TeamMember &Team) { + if (SetBndOpt == SetBoundary::Yes) { + IElement = IElement == 0 ? (NElementsSize - 1) : (IElement - 1); } + const int KMin = getVertBound(VMin, IElement); + const int KMax = getVertBound(VMax, IElement); + const int KRange = KMax - KMin + 1; + parallelForInner( + Team, KRange, INNER_LAMBDA(int KOff) { + const int K = KMin + KOff; + 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 constexpr (Array::rank == 3) { - const int NTracers = ScalarElement.extent_int(0); - const int NVertLayers = ScalarElement.extent_int(2); - parallelFor( - {NTracers, NElementsOwned, NVertLayers}, - 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); + const int NTracers = ScalarElement.extent_int(0); + parallelForOuter( + {NTracers, NElementsSet}, + KOKKOS_LAMBDA(int L, int IElement, const TeamMember &Team) { + if (SetBndOpt == SetBoundary::Yes) { + IElement = IElement == 0 ? (NElementsSize - 1) : (IElement - 1); } + const int KMin = getVertBound(VMin, IElement); + const int KMax = getVertBound(VMax, IElement); + const int KRange = KMax - KMin + 1; + parallelForInner( + Team, KRange, INNER_LAMBDA(int KOff) { + const int K = KMin + KOff; + 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); + } + }); }); } @@ -167,15 +230,28 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, return Err; } +// This overload calls setScalar with vertical bounds based on the array size +template +int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, + const HorzMesh *Mesh, MeshElement Element, + ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes) { + const int VMin = 0; + const int VMax = ScalarElement.extent_int(Array::rank - 1) - 1; + return setScalar(Fun, ScalarElement, Geom, Mesh, Element, VMin, VMax, + ExchangeHalosOpt); +} + enum class CartProjection { Yes, No }; // set vector field on edges based on analytical formula and optionally // exchange halos -template +template int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, EdgeComponent EdgeComp, Geometry Geom, const HorzMesh *Mesh, + const VertMin &VMin, const VertMax &VMax, ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, - CartProjection CartProjectionOpt = CartProjection::Yes) { + CartProjection CartProjectionOpt = CartProjection::Yes, + SetBoundary SetBndOpt = SetBoundary::No) { int Err = 0; @@ -197,6 +273,8 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, auto &AngleEdge = Mesh->AngleEdge; auto &CellsOnEdge = Mesh->CellsOnEdge; auto &VerticesOnEdge = Mesh->VerticesOnEdge; + const int NEdgesSize = Mesh->NEdgesSize; + const int NEdgesSet = Mesh->NEdgesOwned + static_cast(SetBndOpt); auto ProjectVector = KOKKOS_LAMBDA(int IEdge) { Real VecFieldEdge; @@ -278,16 +356,28 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, if constexpr (Array::rank == 1) { parallelFor( - {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { + {NEdgesSet}, KOKKOS_LAMBDA(int IEdge) { + if (SetBndOpt == SetBoundary::Yes) { + IEdge = IEdge == 0 ? (NEdgesSize - 1) : (IEdge - 1); + } VectorFieldEdge(IEdge) = ProjectVector(IEdge); }); } if constexpr (Array::rank == 2) { - const int NVertLayers = VectorFieldEdge.extent_int(1); - parallelFor( - {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int K) { - VectorFieldEdge(IEdge, K) = ProjectVector(IEdge); + parallelForOuter( + {NEdgesSet}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + if (SetBndOpt == SetBoundary::Yes) { + IEdge = IEdge == 0 ? (NEdgesSize - 1) : (IEdge - 1); + } + const int KMin = getVertBound(VMin, IEdge); + const int KMax = getVertBound(VMax, IEdge); + const int KRange = KMax - KMin + 1; + parallelForInner( + Team, KRange, INNER_LAMBDA(int KOff) { + const int K = KMin + KOff; + VectorFieldEdge(IEdge, K) = ProjectVector(IEdge); + }); }); } @@ -300,97 +390,128 @@ int setVectorEdge(const Functor &Fun, const Array &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)); +// This overload calls setVectorEdge with vertical bounds based on the array +// size +template +int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, + EdgeComponent EdgeComp, Geometry Geom, const HorzMesh *Mesh, + ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, + CartProjection CartProjectionOpt = CartProjection::Yes) { - return MaxVal; + const int VMin = 0; + const int VMax = VectorFieldEdge.extent_int(Array::rank - 1) - 1; + return setVectorEdge(Fun, VectorFieldEdge, EdgeComp, Geom, Mesh, VMin, VMax, + ExchangeHalosOpt, CartProjectionOpt); } -inline Real maxVal(const Array2DReal &Arr) { - Real MaxVal; - +template Real reduceArray(const Array1DReal &Arr, int Extent0) { + Real Res; + Reducer R(Res); parallelReduce( - {Arr.extent_int(0), Arr.extent_int(1)}, - KOKKOS_LAMBDA(int I, int J, Real &Accum) { - Accum = Kokkos::max(Arr(I, J), Accum); - }, - Kokkos::Max(MaxVal)); + {Extent0}, KOKKOS_LAMBDA(int I, Real &Accum) { R.join(Accum, Arr(I)); }, + R); - return MaxVal; + return Res; } -inline Real maxVal(const Array3DReal &Arr) { - Real MaxVal; +template Real reduceArray(const Array1DReal &Arr) { + return reduceArray(Arr, Arr.extent_int(0)); +} - parallelReduce( - {Arr.extent_int(0), Arr.extent_int(1), Arr.extent_int(2)}, - KOKKOS_LAMBDA(int I, int J, int K, Real &Accum) { - Accum = Kokkos::max(Arr(I, J, K), Accum); +template +Real reduceArray(const Array2DReal &Arr, int Extent0, const VertMin &VMin, + const VertMax &VMax) { + Real Res; + Reducer ROuter(Res); + + parallelReduceOuter( + {Extent0}, + KOKKOS_LAMBDA(int I, const TeamMember &Team, Real &AccumOuter) { + Real ResInner; + Reducer RInner(ResInner); + + const int KMin = getVertBound(VMin, I); + const int KMax = getVertBound(VMax, I); + const int KRange = KMax - KMin + 1; + + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, Real &AccumInner) { + const int K = KMin + KOff; + RInner.join(AccumInner, Arr(I, K)); + }, + RInner); + + Kokkos::single(PerTeam(Team), + [&]() { ROuter.join(AccumOuter, ResInner); }); }, - Kokkos::Max(MaxVal)); + ROuter); - return MaxVal; + return Res; } -inline Real sum(const Array1DReal &Arr, int Extent0) { - Real Sum; - - parallelReduce( - {Extent0}, KOKKOS_LAMBDA(int I, Real &Accum) { Accum += Arr(I); }, Sum); - - return Sum; +template Real reduceArray(const Array2DReal &Arr, int Extent0) { + return reduceArray(Arr, Extent0, 0, Arr.extent_int(1) - 1); } -inline Real sum(const Array1DReal &Arr) { return sum(Arr, Arr.extent_int(0)); } +template Real reduceArray(const Array2DReal &Arr) { + return reduceArray(Arr, Arr.extent_int(0), 0, + Arr.extent_int(1) - 1); +} -inline Real sum(const Array2DReal &Arr, int Extent0, int Extent1) { - Real Sum; +template +Real reduceArray(const Array3DReal &Arr, int Extent0, int Extent1, + const VertMin &VMin, const VertMax &VMax) { + Real Res; + Reducer ROuter(Res); - parallelReduce( + parallelReduceOuter( {Extent0, Extent1}, - KOKKOS_LAMBDA(int I, int J, Real &Accum) { Accum += Arr(I, J); }, Sum); - - return Sum; -} + KOKKOS_LAMBDA(int L, int I, const TeamMember &Team, Real &AccumOuter) { + Real ResInner; + Reducer RInner(ResInner); + + const int KMin = getVertBound(VMin, I); + const int KMax = getVertBound(VMax, I); + const int KRange = KMax - KMin + 1; + + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, Real &AccumInner) { + const int K = KMin + KOff; + RInner.join(AccumInner, Arr(L, I, K)); + }, + RInner); + + Kokkos::single(PerTeam(Team), + [&]() { ROuter.join(AccumOuter, ResInner); }); + }, + ROuter); -inline Real sum(const Array2DReal &Arr) { - return sum(Arr, Arr.extent_int(0), Arr.extent_int(1)); + return Res; } -inline Real sum(const Array2DReal &Arr, int Extent0) { - return sum(Arr, Extent0, Arr.extent_int(1)); +template +Real reduceArray(const Array3DReal &Arr, int Extent0, int Extent1) { + return reduceArray(Arr, Extent0, Extent1, 0, Arr.extent_int(2) - 1); } -inline Real sum(const Array3DReal &Arr, int Extent0, int Extent1, int Extent2) { - Real Sum; - - parallelReduce( - {Extent0, Extent1, Extent2}, - KOKKOS_LAMBDA(int I, int J, int K, Real &Accum) { - Accum += Arr(I, J, K); - }, - Sum); - - return Sum; +template Real reduceArray(const Array3DReal &Arr, int Extent0) { + return reduceArray(Arr, Extent0, Arr.extent_int(1), 0, + Arr.extent_int(2) - 1); } -inline Real sum(const Array3DReal &Arr) { - return sum(Arr, Arr.extent_int(0), Arr.extent_int(1), Arr.extent_int(2)); +template Real reduceArray(const Array3DReal &Arr) { + return reduceArray(Arr, Arr.extent_int(0), Arr.extent_int(1), 0, + Arr.extent_int(2) - 1); } -inline Real sum(const Array3DReal &Arr, int Extent0) { - return sum(Arr, Extent0, Arr.extent_int(1), Arr.extent_int(2)); +template Real maxVal(ArgTypes &&...Args) { + return reduceArray>(std::forward(Args)...); } -inline Real sum(const Array3DReal &Arr, int Extent0, int Extent1) { - return sum(Arr, Extent0, Extent1, Arr.extent_int(2)); +template Real sum(ArgTypes &&...Args) { + return reduceArray>(std::forward(Args)...); } struct ErrorMeasures { diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 2785876987ae..8d32df9ced66 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -9,6 +9,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -51,8 +52,9 @@ int initState() { int Err = 0; TestSetup Setup; - auto *Mesh = HorzMesh::getDefault(); - auto *State = OceanState::getDefault(); + auto *Mesh = HorzMesh::getDefault(); + auto *VCoord = VertCoord::getDefault(); + auto *State = OceanState::getDefault(); Array2DReal LayerThickCell; Array2DReal NormalVelEdge; @@ -65,21 +67,28 @@ int initState() { int NTracers = Tracers::getNumTracers(); + deepCopy(LayerThickCell, NAN); + deepCopy(NormalVelEdge, NAN); + deepCopy(TracersCell, NAN); + Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThickness(X, Y); }, - LayerThickCell, Geom, Mesh, OnCell); + LayerThickCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, + VCoord->MaxLayerCell, ExchangeHalos::Yes, SetBoundary::Yes); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracersCell, Geom, Mesh, OnCell); + TracersCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, + VCoord->MaxLayerCell, ExchangeHalos::Yes, SetBoundary::Yes); 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, ExchangeHalos::Yes, - CartProjection::No); + NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, ExchangeHalos::Yes, + CartProjection::No, SetBoundary::Yes); return Err; } @@ -104,6 +113,9 @@ int initTendenciesTest(const std::string &mesh) { IO::init(DefComm); Decomp::init(mesh); + // Initialize streams + IOStream::init(); + int HaloErr = Halo::init(); if (HaloErr != 0) { Err++; @@ -111,9 +123,12 @@ int initTendenciesTest(const std::string &mesh) { } VertCoord::init1(); + HorzMesh::init(); Tracers::init(); + VertCoord::init2(); + int StateErr = OceanState::init(); if (StateErr != 0) { Err++; @@ -191,21 +206,24 @@ int testTendencies() { int NTracers = Tracers::getNumTracers(); const Real LayerThickTendSum = - sum(DefTendencies->LayerThicknessTend, NCellsOwned); + sum(DefTendencies->LayerThicknessTend, NCellsOwned, VCoord->MinLayerCell, + VCoord->MaxLayerCell); if (!Kokkos::isfinite(LayerThickTendSum) || LayerThickTendSum == 0) { Err++; LOG_ERROR("TendenciesTest: LayerThickTend FAIL"); } const Real NormVelTendSum = - sum(DefTendencies->NormalVelocityTend, NEdgesOwned); + sum(DefTendencies->NormalVelocityTend, NEdgesOwned, + VCoord->MinLayerEdgeBot, VCoord->MaxLayerEdgeTop); if (!Kokkos::isfinite(NormVelTendSum) || NormVelTendSum == 0) { Err++; LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); } const Real TraceTendSum = - sum(DefTendencies->TracerTend, NTracers, NCellsOwned); + sum(DefTendencies->TracerTend, NTracers, NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); if (!Kokkos::isfinite(TraceTendSum) || TraceTendSum == 0) { Err++; LOG_ERROR("TendenciesTest: TraceTendSum FAIL"); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index e537e9e6341a..529a2aa9ce23 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -338,7 +338,8 @@ int testThickFluxDiv(int NVertLayers, Real RTol) { int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array2DReal ExactThickFluxDiv("ExactThickFluxDiv", Mesh->NCellsOwned, @@ -365,7 +366,7 @@ int testThickFluxDiv(int NVertLayers, Real RTol) { // Compute numerical result Array2DReal NumThickFluxDiv("NumThickFluxDiv", Mesh->NCellsOwned, NVertLayers); - ThicknessFluxDivOnCell ThickFluxDivOnC(Mesh); + ThicknessFluxDivOnCell ThickFluxDivOnC(Mesh, VCoord); parallelFor( {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int KLayer) { ThickFluxDivOnC(NumThickFluxDiv, ICell, KLayer, OnesEdge, @@ -393,7 +394,8 @@ int testPotVortHAdv(int NVertLayers, Real RTol) { int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array2DReal ExactPotVortHAdv("ExactPotVortHAdv", Mesh->NEdgesOwned, @@ -441,7 +443,7 @@ int testPotVortHAdv(int NVertLayers, Real RTol) { // Compute numerical result Array2DReal NumPotVortHAdv("NumPotVortHAdv", Mesh->NEdgesOwned, NVertLayers); - PotentialVortHAdvOnEdge PotVortHAdvOnE(Mesh); + PotentialVortHAdvOnEdge PotVortHAdvOnE(Mesh, VCoord); parallelFor( {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) { PotVortHAdvOnE(NumPotVortHAdv, IEdge, KLayer, NormRelVortEdge, @@ -469,7 +471,8 @@ int testKEGrad(int NVertLayers, Real RTol) { int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array2DReal ExactKEGrad("ExactKEGrad", Mesh->NEdgesOwned, NVertLayers); @@ -491,7 +494,7 @@ int testKEGrad(int NVertLayers, Real RTol) { // Compute numerical result Array2DReal NumKEGrad("NumKEGrad", Mesh->NEdgesOwned, NVertLayers); - KEGradOnEdge KEGradOnE(Mesh); + KEGradOnEdge KEGradOnE(Mesh, VCoord); parallelFor( {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) { KEGradOnE(NumKEGrad, IEdge, KLayer, KECell); @@ -517,7 +520,8 @@ int testSSHGrad(int NVertLayers, Real RTol) { int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array2DReal ExactSSHGrad("ExactSSHGrad", Mesh->NEdgesOwned, NVertLayers); @@ -539,7 +543,7 @@ int testSSHGrad(int NVertLayers, Real RTol) { // Compute numerical result Array2DReal NumSSHGrad("NumSSHGrad", Mesh->NEdgesOwned, NVertLayers); - SSHGradOnEdge SSHGradOnE(Mesh); + SSHGradOnEdge SSHGradOnE(Mesh, VCoord); parallelFor( {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) { SSHGradOnE(NumSSHGrad, IEdge, KLayer, SSHCell); @@ -562,14 +566,15 @@ int testVelDiff(int NVertLayers, Real RTol) { Error Err1; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); Config *OmegaConfig = Config::getOmegaConfig(); Config TendConfig("Tendencies"); Err1 += OmegaConfig->get(TendConfig); CHECK_ERROR_ABORT(Err1, "Tendencies: Tendencies group not found in Config"); - VelocityDiffusionOnEdge VelDiffOnE(Mesh); + VelocityDiffusionOnEdge VelDiffOnE(Mesh, VCoord); Err1 += TendConfig.get("ViscDel2", VelDiffOnE.ViscDel2); CHECK_ERROR_ABORT(Err1, "Tendencies: ViscDel2 not found in TendConfig"); @@ -627,14 +632,15 @@ int testVelHyperDiff(int NVertLayers, Real RTol) { Error Err1; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); Config *OmegaConfig = Config::getOmegaConfig(); Config TendConfig("Tendencies"); Err1 += OmegaConfig->get(TendConfig); CHECK_ERROR_ABORT(Err1, "Tendencies: Tendencies group not found in Config"); - VelocityHyperDiffOnEdge VelHyperDiffOnE(Mesh); + VelocityHyperDiffOnEdge VelHyperDiffOnE(Mesh, VCoord); Err1 += TendConfig.get("ViscDel4", VelHyperDiffOnE.ViscDel4); CHECK_ERROR_ABORT(Err1, "Tendencies: ViscDel4 not found in TendConfig"); @@ -697,7 +703,8 @@ int testWindForcing(int NVertLayers) { int Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); const Real SaltWaterDensity = 0.987654321; @@ -737,7 +744,7 @@ int testWindForcing(int NVertLayers) { // Compute numerical result Array2DReal NumWindForcing("NumWindForcing", Mesh->NEdgesOwned, NVertLayers); - WindForcingOnEdge WindForcingOnE(Mesh); + WindForcingOnEdge WindForcingOnE(Mesh, VCoord); WindForcingOnE.SaltWaterDensity = SaltWaterDensity; parallelFor( @@ -840,7 +847,8 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array3DReal ExactTrFluxDiv("ExactTrFluxDiv", NTracers, Mesh->NCellsOwned, @@ -869,7 +877,7 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { // Compute numerical result Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, NVertLayers); - TracerHorzAdvOnCell TrHorzAdvOnC(Mesh); + TracerHorzAdvOnCell TrHorzAdvOnC(Mesh, VCoord); parallelFor( {NTracers, Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int L, int ICell, int KLayer) { @@ -896,7 +904,8 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array3DReal ExactTracerDiff("ExactTracerDiff", NTracers, Mesh->NCellsOwned, @@ -923,7 +932,7 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) { // Compute numerical result Array3DReal NumTracerDiff("NumTracerDiff", NTracers, Mesh->NCellsOwned, NVertLayers); - TracerDiffOnCell TrDiffOnC(Mesh); + TracerDiffOnCell TrDiffOnC(Mesh, VCoord); TrDiffOnC.EddyDiff2 = 1._Real; parallelFor( @@ -952,7 +961,8 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); // Compute exact result Array3DReal ExactTracerHyperDiff("ExactTracerHyperDiff", NTracers, @@ -973,7 +983,7 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { // Compute numerical result Array3DReal NumTracerHyperDiff("NumTracerHyperDiff", NTracers, Mesh->NCellsOwned, NVertLayers); - TracerHyperDiffOnCell TrHypDiffOnC(Mesh); + TracerHyperDiffOnCell TrHypDiffOnC(Mesh, VCoord); TrHypDiffOnC.EddyDiff4 = 1._Real; parallelFor( {NTracers, Mesh->NCellsOwned, NVertLayers}, @@ -1020,13 +1030,8 @@ void initTendTest(const std::string &MeshFile, int NVertLayers) { } VertCoord::init1(); - // Reset NVertLayers to the test value - auto *DefVertCoord = VertCoord::getDefault(); - DefVertCoord->NVertLayers = NVertLayers; - Dimension::destroy("NVertLayers"); - std::shared_ptr VertDim = - Dimension::create("NVertLayers", NVertLayers); + resetVertCoord(VertCoord::getDefault(), NVertLayers); HorzMesh::init(); @@ -1048,8 +1053,7 @@ int tendencyTermsTest(const std::string &MeshFile = DefaultMeshFile) { initTendTest(MeshFile, NVertLayers); - const auto &Mesh = HorzMesh::getDefault(); - int NTracers = 3; + int NTracers = 3; const Real RTol = sizeof(Real) == 4 ? 2e-2 : 1e-5; diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 7a280d0bcf1e..4ce769d46254 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -172,11 +172,8 @@ int initTimeStepperTest(const std::string &mesh) { // Initialize the vertical coordinate and reset NVertLayers to 1 VertCoord::init1(); - auto *DefVertCoord = VertCoord::getDefault(); - DefVertCoord->NVertLayers = 1; - Dimension::destroy("NVertLayers"); - std::shared_ptr VertDim = - Dimension::create("NVertLayers", NVertLayers); + auto *DefVertCoord = VertCoord::getDefault(); + resetVertCoord(DefVertCoord, NVertLayers); HorzMesh::init(); Tracers::init(); @@ -208,7 +205,7 @@ int initTimeStepperTest(const std::string &mesh) { } auto *TestAuxState = AuxiliaryState::create("TestAuxState", DefMesh, DefHalo, - NVertLayers, NTracers); + DefVertCoord, NTracers); Config *OmegaConfig = Config::getOmegaConfig(); TestAuxState->readConfigOptions(OmegaConfig);