diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index d592acd1359a..7718ac3675e3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -23,8 +23,11 @@ Omega: State: NTimeLevels: 2 Advection: + Coef3rdOrder: 0.25 FluxThicknessType: Center FluxTracerType: Center + VerticalTracerFluxLimiterEnabled: true + VerticalTracerFluxOrder: 3 WindStress: InterpType: Isotropic VertCoord: @@ -50,6 +53,9 @@ Omega: EddyDiff4: 0.0 UseCustomTendency: false ManufacturedSolutionTendency: false + ThicknessVertAdvTendencyEnable: true + VelocityVertAdvTendencyEnable: true + TracerVertAdvTendencyEnable: true Tracers: Base: [Temperature, Salinity] Debug: [Debug1, Debug2, Debug3] diff --git a/components/omega/src/infra/OmegaKokkos.h b/components/omega/src/infra/OmegaKokkos.h index 5a3eb0b82499..a7ab8bca89da 100644 --- a/components/omega/src/infra/OmegaKokkos.h +++ b/components/omega/src/infra/OmegaKokkos.h @@ -78,6 +78,8 @@ using ScratchMemSpace = ExecSpace::scratch_memory_space; using Kokkos::MemoryUnmanaged; using Kokkos::PerTeam; using Kokkos::TeamThreadRange; +using RealScratchArray = + Kokkos::View; /// team_size for hierarchical parallelism #ifdef OMEGA_TARGET_DEVICE @@ -321,7 +323,8 @@ KOKKOS_INLINE_FUNCTION void teamBarrier(const TeamMember &Team) { // parallelForOuter: with label template inline void parallelForOuter(const std::string &Label, - const int (&UpperBounds)[N], F &&Functor) { + const int (&UpperBounds)[N], F &&Functor, + int ScratchValsPerTeam = 0) { auto LinFunctor = LinearIdxWrapper{std::forward(Functor), UpperBounds}; int LinBound = 1; @@ -330,6 +333,12 @@ inline void parallelForOuter(const std::string &Label, } auto Policy = TeamPolicy(LinBound, OMEGA_TEAMSIZE); + + if (ScratchValsPerTeam > 0) { + Policy.set_scratch_size( + 0, Kokkos::PerTeam(ScratchValsPerTeam * sizeof(Real))); + } + Kokkos::parallel_for( Label, Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int TeamId = Team.league_rank(); @@ -339,8 +348,10 @@ inline void parallelForOuter(const std::string &Label, // parallelForOuter: without label template -inline void parallelForOuter(const int (&UpperBounds)[N], F &&Functor) { - parallelForOuter("", UpperBounds, std::forward(Functor)); +inline void parallelForOuter(const int (&UpperBounds)[N], F &&Functor, + int ScratchValsPerTeam = 0) { + parallelForOuter("", UpperBounds, std::forward(Functor), + ScratchValsPerTeam); } // parallelForInner diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 9643fe264cd5..ef3cea652d64 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -196,8 +196,9 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocLayerThicknessAux.computeVarsOnCells(ICell, KChunk, - LayerThickCell); + LocLayerThicknessAux.computeVarsOnCells( + ICell, KChunk, LayerThickCell, NormalVelEdge, 0._Real); + // TODO: make timestep available to this call }); }); Pacer::stop("AuxState:cellAuxState3", 2); diff --git a/components/omega/src/ocn/CustomTendencyTerms.cpp b/components/omega/src/ocn/CustomTendencyTerms.cpp index ff9c17da6435..79ab03630912 100644 --- a/components/omega/src/ocn/CustomTendencyTerms.cpp +++ b/components/omega/src/ocn/CustomTendencyTerms.cpp @@ -10,6 +10,7 @@ #include "Config.h" #include "GlobalConstants.h" #include "TimeStepper.h" +#include "VertCoord.h" namespace OMEGA { @@ -76,8 +77,8 @@ void ManufacturedSolution::init() { /// This test case assumes that the restingThickness is horizontally uniform /// and that only one vertical level is used so only one set of indices is /// used here. - HorzMesh *DefHorzMesh = HorzMesh::getDefault(); - R8 H0 = DefHorzMesh->BottomDepthH(0); + VertCoord *DefVCoord = VertCoord::getDefault(); + R8 H0 = DefVCoord->BottomDepthH(0); // Define and compute common constants R8 Kx = TwoPi / WavelengthX; // Wave in X-dir diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 0c626b49905e..0494273df9cf 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -111,9 +111,6 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh // Read x/y/z and lon/lat coordinates for cells, edges, and vertices readCoordinates(); - // Read the cell-centered bottom depth - readBottomDepth(); - // Read the mesh areas, lengths, and angles readMeasurements(); @@ -439,12 +436,6 @@ void HorzMesh::readCoordinates() { } // end readCoordinates -//------------------------------------------------------------------------------ -// Read the cell-centered bottom depth -void HorzMesh::readBottomDepth() { - readCellArray(BottomDepthH, "bottomDepth"); -} // end readDepth - //------------------------------------------------------------------------------ // Read the mesh areas (cell, triangle, and kite), // lengths (between centers and vertices), and edge angles @@ -600,7 +591,6 @@ void HorzMesh::copyToDevice() { AngleEdge = createDeviceMirrorCopy(AngleEdgeH); WeightsOnEdge = createDeviceMirrorCopy(WeightsOnEdgeH); FVertex = createDeviceMirrorCopy(FVertexH); - BottomDepth = createDeviceMirrorCopy(BottomDepthH); FEdge = createDeviceMirrorCopy(FEdgeH); XCell = createDeviceMirrorCopy(XCellH); YCell = createDeviceMirrorCopy(YCellH); diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index 0b244c0e17a9..424c34e3a5b9 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -42,16 +42,12 @@ class HorzMesh { void readCoordinates(); - void readBottomDepth(); - void readMeasurements(); void readWeights(); void readCoriolis(); - // void computeEdgeSign(); - void copyToDevice(); // int computeMesh(); @@ -232,11 +228,6 @@ class HorzMesh { Array1DReal FVertex; ///< Coriolis parameter at vertices (radians s^-1) HostArray1DReal FVertexH; ///< Coriolis parameter at vertices (radians s^-1) - // Depth - - Array1DReal BottomDepth; ///< Depth of the bottom of the ocean (m) - HostArray1DReal BottomDepthH; ///< Depth of the bottom of the ocean (m) - // Edge sign Array2DReal EdgeSignOnCell; ///< Sign of vector connecting cells diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp new file mode 100644 index 000000000000..886b839025e1 --- /dev/null +++ b/components/omega/src/ocn/VertAdv.cpp @@ -0,0 +1,940 @@ +//===-- ocn/VertAdv.cpp - vertical advection --------------------*- C++ -*-===// +// +// The VertAdv class contains methods for calculating the vertical velocity and +// the tendencies of thickness, horizontal velocity, and tracers due to vertical +// advective transport, and member variables for storing the vertical velocity +// and the fluxes of tracers at vertical interfaces in each column. +// +//===----------------------------------------------------------------------===// + +#include "VertAdv.h" +#include "Field.h" +#include "Tracers.h" + +#include + +namespace OMEGA { + +// create static class members +VertAdv *VertAdv::DefaultVertAdv = nullptr; +std::map> VertAdv::AllVertAdvs; + +//------------------------------------------------------------------------------ +// create the default VertAdv, requires prior initialization of default +// HorzMesh, VertCoord, and Tracers +void VertAdv::init() { + + auto Mesh = HorzMesh::getDefault(); + auto VCoord = VertCoord::getDefault(); + + Config *OmegaConfig = Config::getOmegaConfig(); + + VertAdv::DefaultVertAdv = create("Default", Mesh, VCoord, OmegaConfig); + +} // end init + +//------------------------------------------------------------------------------ +// constructor +VertAdv::VertAdv(const std::string &Name_, //< [in] name for new VertAdv + const HorzMesh *InMesh, //< [in] associated HorzMesh + const VertCoord *InVCoord, //< [in] associated VertCoord + Config *Options //< [in] configuration options +) { + + // Read options from Config object + readConfigOptions(Options); + + // Store name suffix + Name = Name_; + + // Store pointers to Mesh and VertCoord objects + Mesh = InMesh; + VCoord = InVCoord; + + // Store dimension sizes + NVertLayers = VCoord->NVertLayers; + NVertLayersP1 = VCoord->NVertLayersP1; + NCellsOwned = Mesh->NCellsOwned; + NCellsAll = Mesh->NCellsAll; + NCellsSize = Mesh->NCellsSize; + NEdgesOwned = Mesh->NEdgesOwned; + NEdgesAll = Mesh->NEdgesAll; + NEdgesSize = Mesh->NEdgesSize; + NTracers = Tracers::getNumTracers(); + + // Allocate member arrays + VerticalVelocity = + Array2DReal("VerticalVelocity", NCellsSize, NVertLayersP1); + TotalVerticalVelocity = + Array2DReal("TotalVerticalVelocity", NCellsSize, NVertLayersP1); + VertFlux = Array3DReal("VertFlux", NTracers, NCellsSize, NVertLayersP1); + + // Allocate host copies + VerticalVelocityH = createHostMirrorCopy(VerticalVelocity); + TotalVerticalVelocityH = createHostMirrorCopy(TotalVerticalVelocity); + VertFluxH = createHostMirrorCopy(VertFlux); + + // Low-order flux array only needed for flux-corrected transport + if (VertAdvChoice == VertAdvOption::FCT) { + LowOrderVertFlux = + Array3DReal("LowOrderVertFlux", NTracers, NCellsSize, NVertLayersP1); + LowOrderVertFluxH = createHostMirrorCopy(LowOrderVertFlux); + } + + defineFields(); + +} // end constructor + +//------------------------------------------------------------------------------ +// create a new VertAdv instance +VertAdv *VertAdv::create(const std::string &Name, //< [in] name for new VertAdv + const HorzMesh *Mesh, //< [in] associated HorzMesh + const VertCoord *VCoord, //< [in] associated VertCoord + Config *Options //< [in] configuration options +) { + // Check to see if a VertAdv of the same name already exists and, if so, + // exit with an error + if (AllVertAdvs.find(Name) != AllVertAdvs.end()) { + LOG_ERROR("Attempted to create a VertAdv with name {} but a VertAdv " + "of that name already exists", + Name); + return nullptr; + } + + // create a new VertAdv on the heap and put it in a map of unique_ptrs, + // which will manage its lifetime + auto *NewVertAdv = new VertAdv(Name, Mesh, VCoord, Options); + AllVertAdvs.emplace(Name, NewVertAdv); + + return NewVertAdv; + +} // end create + +void VertAdv::defineFields() { + + // set field names (append Name if not default) + VerticalVelocityFldName = "VerticalVelocity"; + TotalVertVelocityFldName = "TotalVerticalVelocity"; + VertFluxFldName = "VertFlux"; + LowOrderVertFluxFldName = "LowOrderVertFlux"; + + if (Name != "Default") { + VerticalVelocityFldName.append(Name); + TotalVertVelocityFldName.append(Name); + VertFluxFldName.append(Name); + LowOrderVertFluxFldName.append(Name); + } + + const Real FillValueReal = -9.99e30; + I4 NDims = 2; + std::vector DimNames(NDims); + DimNames[0] = "NCells"; + DimNames[1] = "NVertLayersP1"; + + auto VerticalVelocityField = Field::create( + VerticalVelocityFldName, // field name + "Vertical pseudovelocity across a pseudoheight surface", // long name or + // description + "m s^-1", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + auto TotalVertVelocityField = + Field::create(TotalVertVelocityFldName, // field name + "Total vertical pseudovelocity across a moving, tilted " + "pseudoheight surface", // long name or description + "m s^-1", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + NDims = 3; + DimNames.resize(NDims); + DimNames[2] = "NTracers"; + + auto VertFluxField = Field::create( + VertFluxFldName, // field name + "Vertical flux of tracers across a pseudoheight surface", // long name or + // description + "", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + auto LowOrderVertFluxField = + Field::create(LowOrderVertFluxFldName, // field name + "Low-order vertical flux of tracers across a pseudoheight " + "surface", // long name or description + "", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + GroupName = "VertAdv"; + if (Name != "Default") { + GroupName.append(Name); + } + + // Create a field group for VertAdv fields + auto VertAdvGroup = FieldGroup::create(GroupName); + + VertAdvGroup->addField(VerticalVelocityFldName); + VertAdvGroup->addField(TotalVertVelocityFldName); + VertAdvGroup->addField(VertFluxFldName); + VertAdvGroup->addField(LowOrderVertFluxFldName); + + // Associate Fields with data + VerticalVelocityField->attachData(VerticalVelocity); + TotalVertVelocityField->attachData(TotalVerticalVelocity); + VertFluxField->attachData(VertFlux); + LowOrderVertFluxField->attachData(LowOrderVertFlux); + +} // end defineFields + +//------------------------------------------------------------------------------ +// destructor +VertAdv::~VertAdv() { + + if (FieldGroup::exists(GroupName)) { + Field::destroy(VerticalVelocityFldName); + Field::destroy(TotalVertVelocityFldName); + Field::destroy(VertFluxFldName); + Field::destroy(LowOrderVertFluxFldName); + FieldGroup::destroy(GroupName); + } + +} // end destructor + +//------------------------------------------------------------------------------ +// Removes all VertAdvs to clean up before exit +void VertAdv::clear() { + + AllVertAdvs.clear(); // removes all VertAdvs from the list and in the + // process, calls the destructors for each + +} // end clear + +//------------------------------------------------------------------------------ +// Removes a VertAdv from map by name +void VertAdv::erase(std::string Name) { + AllVertAdvs.erase(Name); // removes the VertAdv from the list and in the + // process, calls the destructor +} // end erase + +//------------------------------------------------------------------------------ +// Perform deepCopy for each variable array from device to host +void VertAdv::copyToHost() { + + deepCopy(VerticalVelocityH, VerticalVelocity); + deepCopy(TotalVerticalVelocityH, TotalVerticalVelocity); + deepCopy(VertFluxH, VertFlux); + if (VertAdvChoice == VertAdvOption::FCT) { + deepCopy(LowOrderVertFluxH, LowOrderVertFlux); + } +} + +//------------------------------------------------------------------------------ +// Perform deepCopy for each variable array from host to device +void VertAdv::copyToDevice() { + + deepCopy(VerticalVelocity, VerticalVelocityH); + deepCopy(TotalVerticalVelocity, TotalVerticalVelocityH); + deepCopy(VertFlux, VertFluxH); + if (VertAdvChoice == VertAdvOption::FCT) { + deepCopy(LowOrderVertFlux, LowOrderVertFluxH); + } +} + +//------------------------------------------------------------------------------ +// Get default VertAdv +VertAdv *VertAdv::getDefault() { return VertAdv::DefaultVertAdv; } + +//------------------------------------------------------------------------------ +// Get VertAdv by name +VertAdv *VertAdv::get(const std::string Name //< [in] Name of VertAdv +) { + + // look for an instance of this name + auto it = AllVertAdvs.find(Name); + + // if found, return the VertAdv pointer + if (it != AllVertAdvs.end()) { + return it->second.get(); + + // otherwise print error and return null pointer + } else { + LOG_ERROR("VertAdv::get: Attempt to retrieve non-existant VertAdv:"); + LOG_ERROR("{} has not been defined or has been removed", Name); + return nullptr; + } + +} // end get VertAdv + +//------------------------------------------------------------------------------ +// Read and set config options +void VertAdv::readConfigOptions(Config *OmegaConfig) { + + Error Err; // Error code + + Config TendConfig("Tendencies"); + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, "VertAdv: Tendencies group not in Config"); + + Err += TendConfig.get("ThicknessVertAdvTendencyEnable", + this->ThickVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Err += + TendConfig.get("VelocityVertAdvTendencyEnable", this->VelVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("TracerVertAdvTendencyEnable", + this->TracerVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Config AdvectConfig("Advection"); + Err += OmegaConfig->get(AdvectConfig); + CHECK_ERROR_ABORT(Err, "VertAdv: Advection group not in Config"); + + bool FluxLimiterOn; + Err += AdvectConfig.get("VerticalTracerFluxLimiterEnabled", FluxLimiterOn); + CHECK_ERROR_ABORT( + Err, + "VertAdv: VerticalTracerFluxLimiterEnabled not found in AdvectConfig"); + if (FluxLimiterOn) { + VertAdvChoice = VertAdvOption::FCT; + } else { + VertAdvChoice = VertAdvOption::Standard; + } + + I4 VertFluxOrder; + Err += AdvectConfig.get("VerticalTracerFluxOrder", VertFluxOrder); + CHECK_ERROR_ABORT( + Err, "VertAdv: VerticalTracerFluxOrder not found in AdvectConfig"); + + switch (VertFluxOrder) { + case (2): + VertFluxChoice = VertFluxOption::Second; + break; + case (3): + VertFluxChoice = VertFluxOption::Third; + break; + case (4): + VertFluxChoice = VertFluxOption::Fourth; + break; + default: + ABORT_ERROR("VertAdv: Invalid option for VerticalTracerFluxOrder found " + "in AdvectConfig. Must be 2, 3, or 4"); + } + + Err += AdvectConfig.get("Coef3rdOrder", Coef3rdOrder); + CHECK_ERROR_ABORT(Err, "VertAdv: Coef3rdOrder not found in AdvectConfig"); + +} // end readConfigOptions + +//------------------------------------------------------------------------------ +// Compute VerticalVelocity and TotalVerticalVelocity from the horizontal +// velocity (NormalVelocity) and the layer thickness used for fluxes through +// edges (FluxLayerThickEdge) +void VertAdv::computeVerticalVelocity( + const Array2DReal &NormalVelocity, //< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge //< [in] layer thickness at edges +) { + + OMEGA_SCOPE(LocVertVel, VerticalVelocity); + OMEGA_SCOPE(LocNVertLayers, NVertLayers); + OMEGA_SCOPE(LocAreaCell, Mesh->AreaCell); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocNEOnC, Mesh->NEdgesOnCell); + OMEGA_SCOPE(LocEOnC, Mesh->EdgesOnCell); + OMEGA_SCOPE(LocDvE, Mesh->DvEdge); + OMEGA_SCOPE(LocESOnC, Mesh->EdgeSignOnCell); + + // Loop over all cells owned by the task + parallelForOuter( + "computeVerticalVelocity", {NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + RealScratchArray DivHU(Team.team_scratch(0), LocNVertLayers); + + const Real InvAreaCell = 1._Real / LocAreaCell(ICell); + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + I4 KRange = vertRangeChunked(KMin, KMax); + + // Compute thickness-weighted divergence of horizontal velocity + // in each layer + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + Real DivHUTmp[VecLength] = {0}; + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int J = 0; J < LocNEOnC(ICell); ++J) { + const I4 JEdge = LocEOnC(ICell, J); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) * + FluxLayerThickEdge(JEdge, K) * + NormalVelocity(JEdge, K) * InvAreaCell; + } + } + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + DivHU(K) = DivHUTmp[KVec]; + } + }); + + Team.team_barrier(); + + // Set velocity through top and bottom interfaces to zero + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + LocVertVel(ICell, KMin) = 0.; + LocVertVel(ICell, KMax + 1) = 0.; + }); + KRange = vertRange(KMin + 1, KMax); + + // Prefix sum of divergence to determine velocity through each + // interface + parallelScanInner( + Team, KRange, INNER_LAMBDA(int K, Real &Accum, bool IsFinal) { + const I4 KRev = KMax - K; + Accum -= DivHU(KRev); + + if (IsFinal) { + LocVertVel(ICell, KRev) = Accum; + } + }); + }, + NVertLayers); + + // TODO: currently assuming TotalVerticalVelocity = VerticalVelocity, i.e. + // purely from divergence of horizontal velocity. Need to add optional + // corrections to transport velocity from other contributions, e.g. + // movement of vertical interfaces, contribution of horizontal velocity + // through tilted interface. + deepCopy(TotalVerticalVelocity, VerticalVelocity); + +} // end computeVerticalVelocity + +//------------------------------------------------------------------------------ +// Compute thickness tendency due to vertical advection +void VertAdv::computeThicknessVAdvTend( + const Array2DReal &ThickTend //< [inout] thickness tendency +) { + + // Return if vertical advection thickness tendency not enabled + if (!ThickVertAdvEnabled) + return; + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVelocity, TotalVerticalVelocity); + + // Loop over every owned cell, pseudo thickness tendency is simply + // difference in pseudo velocity between bottom and top interface for + // each layer + parallelForOuter( + "computeThicknessVAdvTend", {NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + ThickTend(ICell, K) += LocTotVertVelocity(ICell, K + 1) - + LocTotVertVelocity(ICell, K); + } + }); + }); + +} // end computeThicknessVAdvTend + +//------------------------------------------------------------------------------ +// Compute velocity tendency due to vertical advection +void VertAdv::computeVelocityVAdvTend( + const Array2DReal &VelTend, //< [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, //< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge //< [in] layer thickness at edges +) { + + // Return if vertical advection velocity tendency not enabled + if (!VelVertAdvEnabled) + return; + + OMEGA_SCOPE(LocCOnE, Mesh->CellsOnEdge); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocTotVertVelocity, TotalVerticalVelocity); + OMEGA_SCOPE(EdgeMask, VCoord->EdgeMask); + OMEGA_SCOPE(LocNVertLayersP1, NVertLayersP1); + + // Loop over every owned edge + parallelForOuter( + "computeVelocityVAdvTend", {NEdgesOwned}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const I4 Cell1 = LocCOnE(IEdge, 0); + const I4 Cell2 = LocCOnE(IEdge, 1); + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + I4 KRange = vertRangeChunked(KMin + 1, KMax); + + // Allocate scratch space for W times Du/Dz at vertical interfaces + // between edges + RealScratchArray WDuDzEdge(Team.team_scratch(0), LocNVertLayersP1); + + // Flux is zero at top and bottom + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + WDuDzEdge(KMin) = 0._Real; + WDuDzEdge(KMax + 1) = 0._Real; + }); + + // Average vertical velocities from cell centers to edges and multiply + // by derivative of horizontal velocity to obtain flux of horizontal + // velocity at the vertical interfaces between edges + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + const Real WAvg = 0.5_Real * (LocTotVertVelocity(Cell1, K) + + LocTotVertVelocity(Cell2, K)); + WDuDzEdge(K) = + WAvg * + (NormalVelocity(IEdge, K - 1) - + NormalVelocity(IEdge, K)) / + (0.5_Real * (FluxLayerThickEdge(IEdge, K - 1) + + FluxLayerThickEdge(IEdge, K))); + } + }); + + Team.team_barrier(); + + KRange = vertRangeChunked(KMin, KMax); + // Average W*Du/Dz from interfaces to layer midpoints + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + VelTend(IEdge, K) -= EdgeMask(IEdge, K) * 0.5_Real * + (WDuDzEdge(K) + WDuDzEdge(K + 1)); + } + }); + }, + NVertLayersP1); + +} // end computeVelocityVAdvTend + +//------------------------------------------------------------------------------ +// Compute tracer tendency due to vertical advection, TimeStep is only needed +// as an arugement for flux-corrected transport +void VertAdv::computeTracerVAdvTend( + const Array3DReal &TracerTend, //< [inout] tracer tendencies + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &LayerThickness, //< [in] layer thickness + const TimeInterval TimeStep //< [in] (optional) time step +) { + + // Compute tracer fluxes at the interfaces + computeVerticalFluxes(Tracers, LayerThickness); + + // Dispatch to appropriate algorithm based on configuration settings + switch (VertAdvChoice) { + case VertAdvOption::Standard: + computeStdVAdvTend(TracerTend); + break; + case VertAdvOption::FCT: + R8 Dt; + TimeStep.get(Dt, TimeUnits::Seconds); + computeFCTVAdvTend(TracerTend, Tracers, LayerThickness, Dt); + break; + } + +} // end computeTracerVAdvTend + +//------------------------------------------------------------------------------ +// Compute tracer fluxes due to vertical advection, the particular scheme used +// is chosen via configuration settings +void VertAdv::computeVerticalFluxes( + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &LayerThickness //< [in] layer thickness +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVel, TotalVerticalVelocity); + OMEGA_SCOPE(LocVertFlux, VertFlux); + OMEGA_SCOPE(LocLowOrderVertFlux, LowOrderVertFlux); + OMEGA_SCOPE(LocCoef3rdOrder, Coef3rdOrder); + + // Compute the fluxes used for the standard VAdv scheme, or the high-order + // fluxes used for the FCT VAdv scheme, store in VertFlux member array + switch (VertFluxChoice) { + // 2nd-order centered fluxes + case VertFluxOption::Second: + parallelForOuter( + "computeVerticalFluxes-Second", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + const Real InvLayerThickSum = + 1._Real / (LayerThickness(ICell, K - 1) + + LayerThickness(ICell, K)); + const Real VerticalWeightK = + LayerThickness(ICell, K - 1) * InvLayerThickSum; + const Real VerticalWeightKm1 = + LayerThickness(ICell, K) * InvLayerThickSum; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (VerticalWeightK * Tracers(L, ICell, K) + + VerticalWeightKm1 * Tracers(L, ICell, K - 1)); + } + }); + }); + break; + // 3rd-order upwind fluxes + case VertFluxOption::Third: + parallelForOuter( + "computeVerticalFluxes-Third", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocVertFlux(L, ICell, K) = + (LocTotVertVel(ICell, K) * + (7._Real * (Tracers(L, ICell, K) + + Tracers(L, ICell, K - 1)) - + (Tracers(L, ICell, K + 1) + + Tracers(L, ICell, K - 2))) - + LocCoef3rdOrder * + std::abs(LocTotVertVel(ICell, K)) * + ((Tracers(L, ICell, K + 1) - + Tracers(L, ICell, K - 2)) - + 3._Real * (Tracers(L, ICell, K) - + Tracers(L, ICell, K - 1)))) / + 12._Real; + } + }); + }); + break; + // 4th-order centered fluxes + case VertFluxOption::Fourth: + parallelForOuter( + "computeVerticalFluxes-Fourth", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (7._Real * (Tracers(L, ICell, K) + + Tracers(L, ICell, K - 1)) - + (Tracers(L, ICell, K + 1) + + Tracers(L, ICell, K - 2))) / + 12._Real; + } + }); + }); + break; + } + + // Handle fluxes on interfaces near top and bottom layers. Fluxes are 0 on + // top-most (KMin) and bottom-most (KMax + 1) interfaces, use second order + // for fluxes on next-to-top (KMin + 1) and next-to-bottom (KMax) interfaces. + parallelFor( + "computeVerticalFluxes-TopBot", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + for (int K : {KMin, KMax + 1}) { + LocVertFlux(L, ICell, K) = 0._Real; + } + for (int K : {KMin + 1, KMax}) { + const Real InvLayerThickSum = + 1._Real / + (LayerThickness(ICell, K - 1) + LayerThickness(ICell, K)); + const Real VerticalWeightK = + LayerThickness(ICell, K - 1) * InvLayerThickSum; + const Real VerticalWeightKm1 = + LayerThickness(ICell, K) * InvLayerThickSum; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (VerticalWeightK * Tracers(L, ICell, K) + + VerticalWeightKm1 * Tracers(L, ICell, K - 1)); + } + }); + + // If using FCT scheme, compute 1st-order upwind fluxes for low-order and + // remove low-order flux from high-order flux + if (VertAdvChoice == VertAdvOption::FCT) { + parallelForOuter( + "computeVerticalFluxes-LowOrder", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 1, KMax); + + for (int K : {KMin, KMax + 1}) { + LocLowOrderVertFlux(L, ICell, K) = 0._Real; + } + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocLowOrderVertFlux(L, ICell, K) = + Kokkos::min(0._Real, LocTotVertVel(ICell, K)) * + Tracers(L, ICell, K - 1) + + Kokkos::max(0._Real, LocTotVertVel(ICell, K)) * + Tracers(L, ICell, K); + + LocVertFlux(L, ICell, K) -= + LocLowOrderVertFlux(L, ICell, K); + } + }); + }); + } + +} // end computeVerticalFluxes + +//------------------------------------------------------------------------------ +// Compute tracer tendencies due to vertical advection using standard advection +// scheme +void VertAdv::computeStdVAdvTend( + const Array3DReal &TracerTend //< [inout] tracer tendencies +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocVertFlux, VertFlux); + + // Loop over owned cells, tracer tendency in each layer is computed from + // difference between fluxes through bottom and top interfaces + parallelForOuter( + "computeStdVAdvTend", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + TracerTend(L, ICell, K) += + LocVertFlux(L, ICell, K + 1) - LocVertFlux(L, ICell, K); + } + }); + }); + +} // end computeStdVAdvTend + +//------------------------------------------------------------------------------ +// Compute tracer tendencies due to vertical advection using flux-corrected +// transport scheme. ProvThickness input is provisional layer thickness after +// horizontal thickness flux +void VertAdv::computeFCTVAdvTend( + const Array3DReal &TracerTend, //< [inout] tracer tendencies + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &ProvThickness, //< [in] provisional layer thickness + const Real Dt //< [in] time step +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVel, TotalVerticalVelocity); + OMEGA_SCOPE(LocVertFlux, VertFlux); + OMEGA_SCOPE(LocLowOrderVertFlux, LowOrderVertFlux); + OMEGA_SCOPE(LocNVertLayers, NVertLayers); + OMEGA_SCOPE(LocEps, Eps); + + parallelForOuter( + "computeFCTVAdvTend", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + I4 KRange = vertRangeChunked(KMin, KMax); + + RealScratchArray InvNewProvThick(Team.team_scratch(0), + LocNVertLayers); + RealScratchArray WorkTend(Team.team_scratch(0), LocNVertLayers); + RealScratchArray FlxIn(Team.team_scratch(0), LocNVertLayers); + RealScratchArray FlxOut(Team.team_scratch(0), LocNVertLayers); + RealScratchArray RescaledFlux(Team.team_scratch(0), + LocNVertLayers + 1); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + InvNewProvThick(K) = + 1._Real / (ProvThickness(ICell, K) + + Dt * (LocTotVertVel(ICell, K + 1) - + LocTotVertVel(ICell, K))); + Real TracerMax; + Real TracerMin; + // Determine bounds on tracer from neighbor values for + // limiting + if (K == KMin) { + TracerMax = Kokkos::max(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1)); + TracerMin = Kokkos::min(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1)); + } else if (K == KMax) { + TracerMax = Kokkos::max(Tracers(L, ICell, K - 1), + Tracers(L, ICell, K)); + TracerMin = Kokkos::min(Tracers(L, ICell, K - 1), + Tracers(L, ICell, K)); + } else { + TracerMax = + Kokkos::max(Tracers(L, ICell, K - 1), + Kokkos::max(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1))); + TracerMin = + Kokkos::min(Tracers(L, ICell, K - 1), + Kokkos::min(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1))); + } + + // Accumulate upwind flux in WorkTend + WorkTend(K) = LocLowOrderVertFlux(L, ICell, K + 1) - + LocLowOrderVertFlux(L, ICell, K); + // Accumulate remaining high-order flux into layer + FlxIn(K) = + Kokkos::max(0._Real, LocVertFlux(L, ICell, K + 1)) - + Kokkos::min(0._Real, LocVertFlux(L, ICell, K)); + // Accumulate remaining high-order flux out of layer + FlxOut(K) = + Kokkos::min(0._Real, LocVertFlux(L, ICell, K + 1)) - + Kokkos::max(0._Real, LocVertFlux(L, ICell, K)); + // Build scale factors to limit flux for FCT using the + // bounds determined above and bounds on newly updated + // values. Factors are stored in FlxIn and FlxOut + // scratch space. + Real TracerMinNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * (WorkTend(K) + FlxOut(K))) * + InvNewProvThick(K); + Real TracerMaxNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * (WorkTend(K) + FlxIn(K))) * + InvNewProvThick(K); + Real TracerUpwindNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * WorkTend(K)) * + InvNewProvThick(K); + Real ScaleFactor = + (TracerMax - TracerUpwindNew) / + (TracerMaxNew - TracerUpwindNew + LocEps); + FlxIn(K) = + Kokkos::min(1._Real, Kokkos::max(0._Real, ScaleFactor)); + ScaleFactor = (TracerUpwindNew - TracerMin) / + (TracerUpwindNew - TracerMinNew + LocEps); + FlxOut(K) = + Kokkos::min(1._Real, Kokkos::max(0._Real, ScaleFactor)); + } + }); + + Team.team_barrier(); + + KRange = vertRangeChunked(KMin + 1, KMax); + + // Rescale the high-order vertical flux + RescaledFlux(KMin) = 0._Real; + RescaledFlux(KMax + 1) = 0._Real; + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + + RescaledFlux(K) = + Kokkos::max(0._Real, LocVertFlux(L, ICell, K)) * + Kokkos::min(FlxOut(K), FlxIn(K - 1)) + + Kokkos::min(0._Real, LocVertFlux(L, ICell, K)) * + Kokkos::min(FlxOut(K - 1), FlxIn(K)); + } + }); + + Team.team_barrier(); + + // Accumulate total FCT vertical advection tendency + KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + WorkTend(K) += RescaledFlux(K + 1) - RescaledFlux(K); + TracerTend(L, ICell, K) += WorkTend(K); + } + }); + // TODO: Monotonicity and diagnostic checks + }, + 5 * NVertLayers + 1); + +} // end computeFTCVAdvTend + +} // end namespace OMEGA +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h new file mode 100644 index 000000000000..95bd30968913 --- /dev/null +++ b/components/omega/src/ocn/VertAdv.h @@ -0,0 +1,224 @@ +#ifndef OMEGA_VERTADV_H +#define OMEGA_VERTADV_H +//===-- ocn/VertAdv.h - vertical advection definitions ----------*- C++ -*-===// +// +/// \file +/// \brief Contains methods and variables for vertical advection +/// +/// This header defines the VertAdv class which contains methods and variables +/// for calculating the vertical velocity and for vertical transport of +/// thickness, velocity, and tracers. +// +//===----------------------------------------------------------------------===// + +#include "AuxiliaryState.h" +#include "Config.h" +#include "DataTypes.h" +#include "Error.h" +#include "HorzMesh.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "TimeMgr.h" +#include "VertCoord.h" + +#include +#include +#include + +namespace OMEGA { + +// Enum for flux types +enum class VertFluxOption { Second, Third, Fourth }; +// Enum for advection scheme type +enum class VertAdvOption { Standard, FCT }; + +class VertAdv { + + public: + // Logicals to enable/disable advection of specific fields + bool ThickVertAdvEnabled = false; + bool VelVertAdvEnabled = false; + bool TracerVertAdvEnabled = false; + + // enums storing options chosen by user + VertFluxOption VertFluxChoice; + VertAdvOption VertAdvChoice; + + // Mesh dimensions + I4 NVertLayers; + I4 NVertLayersP1; + I4 NCellsOwned; + I4 NCellsAll; + I4 NCellsSize; + I4 NEdgesOwned; + I4 NEdgesAll; + I4 NEdgesSize; + + // Number of tracers + I4 NTracers; + + // Coefficient for blending 3rd-order and 4th-order reconstruction of + // tracers with VertFluxOption::Third. Coef3rdOrder = 1 gives purely + // 3rd-order, Coef3rdOrder = 0 gives purely 4th-order. + Real Coef3rdOrder; + + Array2DReal VerticalVelocity; ///< pseudovelocity through top of cell + Array2DReal + TotalVerticalVelocity; ///< transport velocity through top of Cell + Array3DReal VertFlux; ///< fluxes at vertical interfaces + Array3DReal LowOrderVertFlux; ///< low-order fluxes for FCT + + // Arrays on host + HostArray2DReal VerticalVelocityH; + HostArray2DReal TotalVerticalVelocityH; + HostArray3DReal VertFluxH; + HostArray3DReal LowOrderVertFluxH; + + // VertAdv instance name and group name + std::string Name; + std::string GroupName; + + // Field names + std::string VerticalVelocityFldName; + std::string TotalVertVelocityFldName; + std::string VertFluxFldName; + std::string LowOrderVertFluxFldName; + + // public methods + + // Initialize the default VertAdv instance + static void init(); + + /// Creates a new vertical advection object by calling the constructor and + /// puts it in the AllVertAdvs map. + static VertAdv * + create(const std::string &Name, ///< [in] name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options + ); + + /// Copy member arrays from device to host + void copyToHost(); + + /// Copy member arrays from host to device + void copyToDevice(); + + /// Destructor - deallocates all memory and deletes a VertAdv + ~VertAdv(); + + /// Deallocates arrays + static void clear(); + + /// Remove a VertAdv by name + static void erase(std::string InName); + + /// Retrieve the default VertAdv + static VertAdv *getDefault(); + + /// Retreive a VertAdv by name + static VertAdv *get(std::string name); + + /// Read and set config options + void readConfigOptions(Config *Options); + + /// Determine transport due to vertical advection from divergence of + /// horizontal advection. + void computeVerticalVelocity( + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges + ); + + /// Compute pseudo thickness tendency due to vertical advection + void computeThicknessVAdvTend( + const Array2DReal &ThickTend ///< [inout] thickness tendency + ); + + /// Compute velocity tendency due to vertical advection + void computeVelocityVAdvTend( + const Array2DReal &VelTend, ///< [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges + ); + + /// Compute tracer tendency due to vertical advection, The LayerThickness at + /// the beginning of the time step is passed for standard advection, while + /// the provisional thickness including horizontal thickness flux is passed + /// for flux-corrected transport. The TimeStep is only needed as an argument + /// for flux-corrected transport + void computeTracerVAdvTend( + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness, ///< [in] layer thickness + const TimeInterval TimeStep = + TimeInterval() ///< [in] (optional) time step + ); + + /// Compute tracer fluxes due to vertical advection + void computeVerticalFluxes( + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness ///< [in] layer thickness + ); + + /// Compute tracer tendencies due to vertical advection using standard + /// advection scheme + void computeStdVAdvTend( + const Array3DReal &TracerTend ///< [inout] tracer tendencies + ); + + /// Compute tracer tendencies due to vertical advection using flux-corrected + /// transport scheme + void computeFCTVAdvTend( + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &ProvThickness, ///< [in] provisional layer thickness + const Real Dt ///< [in] time step + ); + + private: + // Vertical loop bounds from VertCoord + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + + // Arrays from HorzMesh + Array2DI4 CellsOnEdge; + Array2DI4 EdgesOnCell; + Array1DI4 NEdgesOnCell; + Array1DReal AreaCell; + Array1DReal DvEdge; + Array2DReal EdgeSignOnCell; + + // small number to avoid numerical difficulties in computing + // scale factors for FCT. + Real Eps = 1e-10; + + const HorzMesh *Mesh; + const VertCoord *VCoord; + + static VertAdv *DefaultVertAdv; + static std::map> AllVertAdvs; + + // private methods + + /// construct a new vertical coordinate object + VertAdv(const std::string &Name, ///< [in] Name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options + ); + + /// define field metadata + void defineFields(); + + // Forbid copy and move construction + VertAdv(const VertAdv &) = delete; + VertAdv(VertAdv &&) = delete; + +}; // end class VertAdv + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // defined OMEGA_VERTICALADV_H diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index cb937863528c..8a69aa7c8269 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -14,7 +14,12 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, Mesh->NEdgesSize, VCoord->NVertLayers), SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize, VCoord->NVertLayers), - CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth), + ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize, + VCoord->NVertLayers), + AreaCell(Mesh->AreaCell), DvEdge(Mesh->DvEdge), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), + BottomDepth(VCoord->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} @@ -78,21 +83,37 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName, DimNames // dimension names ); + // Provisional Thickness + auto ProvThicknessField = Field::create( + ProvThickness.label(), // field name + "layer thickness after horizontal flux", // long Name or description + "m", // units + "", // CF standard Name + 0, // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + // Add fields to Aux field group FieldGroup::addFieldToGroup(FluxLayerThickEdge.label(), AuxGroupName); FieldGroup::addFieldToGroup(MeanLayerThickEdge.label(), AuxGroupName); FieldGroup::addFieldToGroup(SshCell.label(), AuxGroupName); + FieldGroup::addFieldToGroup(ProvThickness.label(), AuxGroupName); // Attach field data FluxLayerThickEdgeField->attachData(FluxLayerThickEdge); MeanLayerThickEdgeField->attachData(MeanLayerThickEdge); SshCellField->attachData(SshCell); + ProvThicknessField->attachData(ProvThickness); } void LayerThicknessAuxVars::unregisterFields() const { Field::destroy(FluxLayerThickEdge.label()); Field::destroy(MeanLayerThickEdge.label()); Field::destroy(SshCell.label()); + Field::destroy(ProvThickness.label()); } } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h index 0620d546fd8f..a415bfa02ba3 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h @@ -17,6 +17,7 @@ class LayerThicknessAuxVars { Array2DReal FluxLayerThickEdge; Array2DReal MeanLayerThickEdge; Array2DReal SshCell; + Array2DReal ProvThickness; FluxThickEdgeOption FluxThickEdgeChoice; @@ -63,9 +64,10 @@ class LayerThicknessAuxVars { } } - KOKKOS_FUNCTION void - computeVarsOnCells(int ICell, int KChunk, - const Array2DReal &LayerThickCell) const { + KOKKOS_FUNCTION void computeVarsOnCells(int ICell, int KChunk, + const Array2DReal &LayerThickCell, + const Array2DReal &NormalVelEdge, + const Real Dt) const { // Temporary for stacked shallow water const int KStart = chunkStart(KChunk, MinLayerCell(ICell)); @@ -76,6 +78,25 @@ class LayerThicknessAuxVars { SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell); } + Real TmpProv[VecLength] = {0.}; + + Real DtInvAreaCell = Dt / AreaCell(ICell); + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + const I4 JEdge = EdgesOnCell(ICell, J); + const Real Factor = + DtInvAreaCell * DvEdge(JEdge) * EdgeSignOnCell(ICell, J); + for (int KVec = 0; KVec < KLen; ++KVec) { + const int K = KStart + KVec; + TmpProv[KVec] += + Factor * FluxLayerThickEdge(JEdge, K) * NormalVelEdge(JEdge, K); + } + } + + for (int KVec = 0; KVec < KLen; ++KVec) { + const int K = KStart + KVec; + ProvThickness(ICell, K) = LayerThickCell(ICell, K) + TmpProv[KVec]; + } + /* Real TotalThickness = 0.0; for (int K = 0; K < NVertLayers; K++) { @@ -91,6 +112,11 @@ class LayerThicknessAuxVars { void unregisterFields() const; private: + Array1DReal AreaCell; + Array1DReal DvEdge; + Array1DI4 NEdgesOnCell; + Array2DI4 EdgesOnCell; + Array2DReal EdgeSignOnCell; Array2DI4 CellsOnEdge; Array1DReal BottomDepth; Array1DI4 MinLayerEdgeBot; diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 2f4b182f3ed0..3aeb12df79e9 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -476,3 +476,14 @@ add_omega_test( ocn/VertMixTest.cpp "-n;4" ) + +################## +# VAdv test +################## + +add_omega_test( + VERTADV_TEST + testVertAdv.exe + ocn/VertAdvTest.cpp + "-n;1" +) diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 0c218325b100..f178b3e8f171 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -368,28 +368,6 @@ int main(int argc, char *argv[]) { LOG_INFO("HorzMeshTest: Vertex lon/lat test PASS"); } - // Test bounds of bathymetry - // Find minimum and maximum values of the bottom depth - // and compares to reasonable values - // Tests that the bottom depth has been read in correctly - R8 MaxBathy = -1e10; - R8 MinBathy = 1e10; - for (int Cell = 0; Cell < LocCells; Cell++) { - if (Mesh->BottomDepthH(Cell) < MinBathy) { - MinBathy = Mesh->BottomDepthH(Cell); - } - if (Mesh->BottomDepthH(Cell) > MaxBathy) { - MaxBathy = Mesh->BottomDepthH(Cell); - } - } - - if ((MinBathy > 0) && (MaxBathy < 11000.0)) { - LOG_INFO("HorzMeshTest: Bathy min/max test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Bathy min/max test FAIL"); - } - // Test cell areas // Find the global sum of all the local cell areas // and compares to reasonable value for Earth's ocean area diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp new file mode 100644 index 000000000000..18b1bea12f9f --- /dev/null +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -0,0 +1,441 @@ +//===-- Test driver for OMEGA Vertical Advection ----------------*- C++ -*-===/ +// +/// \file +/// \brief Test driver for OMEGA VertAdv class +/// +/// +// +//===-----------------------------------------------------------------------===/ + +#include "VertAdv.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "GlobalConstants.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "TimeStepper.h" +#include "VertCoord.h" + +#include "Pacer.h" +#include +#include +#include + +using namespace OMEGA; + +constexpr Real Alpha = 0.197_Real; // random phase angle +constexpr Real Length = 1024._Real; + +// Test function used for tendTest +KOKKOS_FUNCTION Real testFunc(const Real X) { + Real Arg = Pi * (X + Alpha) / Length; + return std::sin(Arg) * std::sin(Arg) * std::cos(Arg) * std::cos(Arg); +} + +// Function used for analytical derivative in tendTest +KOKKOS_FUNCTION Real testDeriv(const Real X) { + Real C = Pi / Length; + return (C / 2._Real) * std::sin(4._Real * C * (X + Alpha)); +} + +// For a given resolution, use computeStdVAdvTend to compute tracer +// tendencies, then calculate L2 error +Real tendTest(const int NLayers, VertAdv *VAdv) { + + // Setup test for this resolution + Real Delta = Length / NLayers; + + VAdv->NVertLayers = NLayers; + VAdv->TotalVerticalVelocity = + Array2DReal("TotalVerticaVelocity", 1, NLayers + 1); + Array1DReal XLayer("XLayer", NLayers); + OMEGA_SCOPE(LocTotVertVel, VAdv->TotalVerticalVelocity); + Array2DReal LayerThick("LayerThickness", 1, NLayers); + Array3DReal TracerArray("TracerArray", 1, 1, NLayers); + Array3DReal Tend("TracerArray", 1, 1, NLayers); + + // Set uniform velocity throughout domain + deepCopy(VAdv->TotalVerticalVelocity, 1._Real); + // Set velocities at top and bottom interfaces to 0. + parallelFor( + {1}, KOKKOS_LAMBDA(const int &) { + LocTotVertVel(0, 0) = 0._Real; + LocTotVertVel(0, NLayers) = 0._Real; + }); + // Set X values, LayerThickness, and Tracer values throughout domain + parallelFor( + {NLayers}, KOKKOS_LAMBDA(int K) { + XLayer(K) = Delta * (K + 0.5_Real); + LayerThick(0, K) = Delta; + TracerArray(0, 0, K) = testFunc(XLayer(K)); + }); + + // Compute tracer tendencies + VAdv->computeTracerVAdvTend(Tend, TracerArray, LayerThick); + + HostArray3DReal TendH = createHostMirrorCopy(Tend); + HostArray1DReal XLayerH = createHostMirrorCopy(XLayer); + Real SumSq = 0._Real; + int Count = 0; + + // Compute errors at interior of range, ignoring layers near top and bottom + // where lower order is used + for (int K = 2; K < NLayers - 2; ++K) { + Real Expected = testDeriv(XLayerH(K)); + Real Diff = TendH(0, 0, K) / Delta - Expected; + SumSq += Diff * Diff; + Count++; + } + Real L2 = std::sqrt(SumSq / Count); + + return L2; +} + +// Initialize needed modules +void initVertAdvTest() { + + I4 Err; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + // Initialize the Logging system + initLogging(DefEnv); + + // Open config file + Config("Omega"); + Config::readAll("omega.yml"); + + // First step of time stepper initialization needed for IOstream + TimeStepper::init1(); + + // Initialize the IO system + IO::init(DefComm); + + // Create the default decomposition (initializes the decomposition) + Decomp::init(); + + // Initialize streams + IOStream::init(); + + // Initialize the default halo + Err = Halo::init(); + if (Err != 0) + ABORT_ERROR("VertAdvTest: error initializing default halo"); + + // Initialize the default mesh + HorzMesh::init(); + + // Initialize the default vertical coordinate + VertCoord::init(false); + + // Initialize tracers + Tracers::init(); + + // Initialize the default vertical advection + VertAdv::init(); +} + +// Clean-up modules +void finalizeVertAdvTest() { + + IOStream::finalize(); + Tracers::clear(); + VertAdv::clear(); + VertCoord::clear(); + TimeStepper::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); +} + +int main(int argc, char *argv[]) { + + // Initialize error code + Error ErrAll; + + // Initialize the global MPI environment + MPI_Init(&argc, &argv); + Kokkos::initialize(); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + { + initVertAdvTest(); + + auto DefMesh = HorzMesh::getDefault(); + auto DefDecomp = Decomp::getDefault(); + auto DefVertAdv = VertAdv::getDefault(); + auto DefVCoord = VertCoord::getDefault(); + + // Only need to test within one column or along one edge and with one + // tracer + DefVertAdv->NCellsOwned = 1; + DefVertAdv->NCellsAll = 1; + DefVertAdv->NEdgesOwned = 1; + DefVertAdv->NEdgesAll = 1; + DefVertAdv->NTracers = 1; + + I4 NVertLayers = DefVertAdv->NVertLayers; + I4 NVertLayersP1 = DefVertAdv->NVertLayersP1; + + Array2DReal FluxLayerThickEdge("FluxLayerThickEdge", + DefDecomp->NEdgesSize, NVertLayers); + Array2DReal NormalVelEdge("NormalVelEdge", DefDecomp->NEdgesSize, + NVertLayers); + + const I4 ICell0 = 0; + const I4 NEdgesOnCell0 = DefDecomp->NEdgesOnCellH(ICell0); + + // Test for computeVerticalVelocity + I4 Err = 0; + + OMEGA_SCOPE(LocEOnC, DefDecomp->EdgesOnCell); + OMEGA_SCOPE(LocESOnC, DefMesh->EdgeSignOnCell); + + // Set unifrom values for layer thickness and velocity on edges of column + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(int K) { + for (int J = 0; J < NEdgesOnCell0; ++J) { + const I4 JEdge = LocEOnC(ICell0, J); + FluxLayerThickEdge(JEdge, K) = 1._Real; + NormalVelEdge(JEdge, K) = 1._Real * LocESOnC(ICell0, J); + } + }); + + // Compute vertical velocity + DefVertAdv->computeVerticalVelocity(NormalVelEdge, FluxLayerThickEdge); + + HostArray2DReal VertVelH = + createHostMirrorCopy(DefVertAdv->VerticalVelocity); + + Real Tol = 1e-10; + + // Divergence is equal in each layer. Expected vertical velocity through + // each interface is sum over negative divergence starting from bottom + // layer + Real Perim0 = 0.; + for (int J = 0; J < NEdgesOnCell0; ++J) { + Perim0 += DefMesh->DvEdgeH(DefMesh->EdgesOnCellH(ICell0, J)); + } + Real InvAreaCell0 = 1._Real / DefMesh->AreaCellH(ICell0); + + for (int K = 1; K < NVertLayersP1; ++K) { + Real Expected = (NVertLayers - K) * Perim0 * InvAreaCell0; + Real Diff = std::abs(VertVelH(ICell0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeVerticalVelocity FAIL"); + } + + // Test for computeThicknessVAdvTend + Err = 0; + + Array2DReal TendCell2D("TendCell2D", DefMesh->NCellsSize, NVertLayers); + // Compute thickness tendencies for each layer using velocities from + // previous tests + DefVertAdv->computeThicknessVAdvTend(TendCell2D); + HostArray2DReal TendCell2DH = createHostMirrorCopy(TendCell2D); + + // Expected thickness tendency for each layer is difference between + // velocity at bottom and top interface, i.e. the divergence + for (int K = 1; K < NVertLayers; ++K) { + Real Expected = -Perim0 * InvAreaCell0; + Real Diff = std::abs(TendCell2DH(ICell0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeThicknessVAdvTend FAIL"); + } + + // Test for computeVelocityVAdvTend + Err = 0; + + // Set uniform vertical velocity in each column along the edge and + // uniform layer thickness at the edge. Horizontal velocity is set + // to a periodic function + Array2DReal TendEdge2D("TendEdge2D", DefMesh->NEdgesSize, NVertLayers); + const I4 IEdge0 = 0; + OMEGA_SCOPE(LocCOnE, DefDecomp->CellsOnEdge); + OMEGA_SCOPE(LocVertVel, DefVertAdv->TotalVerticalVelocity); + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(int K) { + NormalVelEdge(IEdge0, K) = sin(Pi * K / NVertLayers); + FluxLayerThickEdge(IEdge0, K) = 1._Real; + const I4 Cell1 = LocCOnE(IEdge0, 0); + const I4 Cell2 = LocCOnE(IEdge0, 1); + LocVertVel(Cell1, K) = 1._Real; + LocVertVel(Cell2, K) = 1._Real; + }); + + // Compute velocity tendencies for each layer + DefVertAdv->computeVelocityVAdvTend(TendEdge2D, NormalVelEdge, + FluxLayerThickEdge); + + Tol = 1._Real / (NVertLayers * NVertLayers); + + HostArray2DReal TendEdge2DH = createHostMirrorCopy(TendEdge2D); + + // Expected velocity tendency is derivative of initial distribution + for (int K = 0; K < NVertLayers; ++K) { + Real Expected = (Pi / NVertLayers) * cos(Pi * K / NVertLayers); + if (K == 0 or K == NVertLayers - 1) + Expected /= 2._Real; + Real Diff = std::abs(TendEdge2DH(IEdge0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeVelocityVAdvTend FAIL"); + } + + // Tests for computeTracerVAdvTend + + DefVertAdv->VertAdvChoice = VertAdvOption::Standard; + + // Setup for convergence tests with computeStdVAdvTend + constexpr I4 NLayersArray[] = {128, 256, 512, 1024}; + constexpr VertFluxOption AllOrders[] = {VertFluxOption::Second, + VertFluxOption::Third, + VertFluxOption::Fourth}; + static const std::map OrderOfAcc = { + {VertFluxOption::Second, 2}, + {VertFluxOption::Third, 3}, + {VertFluxOption::Fourth, 4}}; + static const std::map OrderStr = { + {VertFluxOption::Second, "2nd"}, + {VertFluxOption::Third, "3rd"}, + {VertFluxOption::Fourth, "4th"}}; + + // Set Coef3rdOrder to 1 so VertFluxOption::Third is calculated with + // purely 3rd-order and not a blend of 3rd- and 4th-order fluxes + DefVertAdv->Coef3rdOrder = 1.; + + Tol = 0.1_Real; + OMEGA_SCOPE(LocMaxLyrCell, DefVCoord->MaxLayerCell); + + // Compute L2 error for computeStdVAdvTend with each VertFluxOption over a + // set of resolutions. Compare successive errors for each Order to confirm + // expected order of accuracy. + for (VertFluxOption Order : AllOrders) { + Err = 0; + DefVertAdv->VertFluxChoice = Order; + std::vector L2Errors; + for (I4 NLayers : NLayersArray) { + parallelFor( + {1}, + KOKKOS_LAMBDA(const int &) { LocMaxLyrCell(0) = NLayers - 1; }); + Real L2Err = tendTest(NLayers, DefVertAdv); + L2Errors.push_back(L2Err); + } + Real ExpectedRat = std::pow(2._Real, OrderOfAcc.at(Order)); + for (I4 I = 0; I < L2Errors.size() - 1; ++I) { + Real Rat = L2Errors[I] / L2Errors[I + 1]; + Real RelErr = std::abs(Rat - ExpectedRat) / ExpectedRat; + if (RelErr > Tol) { + ++Err; + } + } + if (Err != 0) { + ErrAll += + Error(ErrorCode::Fail, + "VertAdvTest: computeStdVAdvTend with {} order FAIL", + OrderStr.at(Order)); + } + } + + // Monotonicity check for computeFCTVAdvTend + + NVertLayers = 256; + DefVertAdv->VertAdvChoice = VertAdvOption::FCT; + DefVertAdv->VertFluxChoice = VertFluxOption::Third; + DefVertAdv->Coef3rdOrder = 0.25; + std::vector TestVel = {1._Real, -1._Real}; + std::vector StartIdx = {NVertLayers / 2, NVertLayers / 4}; + parallelFor( + {1}, + KOKKOS_LAMBDA(const int &) { LocMaxLyrCell(0) = NVertLayers - 1; }); + + // Test a top-hat function for a uniform velocity in both directions + for (int IMono = 0; IMono < 2; ++IMono) { + Err = 0; + Array3DReal TracerArray("TracerArray", 1, 1, NVertLayers); + Array2DReal LayerThick("LayerThick", 1, NVertLayers); + Array3DReal TendCell3D("TendCell3D", 1, 1, NVertLayers); + deepCopy(LayerThick, 1._Real); + DefVertAdv->NVertLayers = NVertLayers; + LocVertVel = Array2DReal("TotalVerticaVelocity", 1, NVertLayers + 1); + deepCopy(LocVertVel, TestVel[IMono]); + parallelFor( + {1}, KOKKOS_LAMBDA(const int &) { + LocVertVel(0, 0) = 0._Real; + LocVertVel(0, NVertLayers) = 1._Real; + }); + deepCopy(TracerArray, 0._Real); + const I4 KStart = StartIdx[IMono]; + const I4 KRange = NVertLayers / 4; + parallelFor( + {KRange}, KOKKOS_LAMBDA(const int K) { + TracerArray(0, 0, KStart + K) = 1._Real; + }); + Real Dt = 0.5_Real; + TimeInterval TimeStep(Dt, TimeUnits::Seconds); + + const I4 NTSteps = 100; + + // Integrate 100 time steps and confirm monotonicity + for (int N = 0; N < NTSteps; ++N) { + deepCopy(TendCell3D, 0._Real); + DefVertAdv->computeTracerVAdvTend(TendCell3D, TracerArray, + LayerThick, TimeStep); + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(const int K) { + TracerArray(0, 0, K) += Dt * TendCell3D(0, 0, K); + }); + } + HostArray3DReal TracerH = createHostMirrorCopy(TracerArray); + HostArray3DReal TendH = createHostMirrorCopy(TendCell3D); + for (int K = 0; K < NVertLayers; ++K) { + if (TracerH(0, 0, K) < 0._Real or TracerH(0, 0, K) > 1._Real) { + ++Err; + } + } + if (Err != 0) { + ErrAll += + Error(ErrorCode::Fail, + "VertAdvTest: computeFCTVAdvTend with velocity = {} " + "monotonicity FAIL", + TestVel[IMono]); + } + } + + finalizeVertAdvTest(); + } + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + CHECK_ERROR_ABORT(ErrAll, "VertAdv unit tests FAIL"); + + return 0; +}