diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..a62bbe2754b7 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -98,6 +98,12 @@ Omega: RiCrit: 0.7 Exponent: 3.0 RiSmoothLoops: 2 + Submeso: + Enable: false + Tau: 172800.0 + Ce: 0.08 + LfMin: 1.0e3 + DsMax: 100.0e3 IOStreams: HorzMeshIn: UsePointerFile: false diff --git a/components/omega/doc/devGuide/SubmesoEddies.md b/components/omega/doc/devGuide/SubmesoEddies.md new file mode 100644 index 000000000000..685d5756dfcf --- /dev/null +++ b/components/omega/doc/devGuide/SubmesoEddies.md @@ -0,0 +1,73 @@ +(omega-dev-submeso-eddies)= + +# Submesoscale Eddy Parameterization + +Omega includes a `SubmesoEddies` class (`src/ocn/SubmesoEddies.{h,cpp}`) that implements the +submesoscale mixed-layer instability closure used to produce an eddy-induced +normal transport velocity. It provides methods to compute mixed-layer depth from +a density threshold, buoyancy gradients, and eddy velocity. +The current implementation follows the Fox-Kemper et al. (2011) (FK11) closure, +as described in the +{ref}`omega-design-submesoscale-eddies` design document. + +## Initialization + +The `SubmesoEddies` class is implemented as a singleton. Before creating it, +[`HorzMesh`](#omega-dev-horz-mesh) and [`VertCoord`](#omega-dev-vert-coord) +must be initialized. Create the instance with the static method +```c++ +SubmesoEddies::init(); +``` +Retrieve the pointer at any time with: +```c++ +SubmesoEddies* DefSubEddies = SubmesoEddies::getInstance(); +``` + +## Data and algorithms + +Core public fields include: + +- `DenMixLayerIndex`, `DenMixLayerDepth` +- `GradBuoyEdgeInterface` +- `EddyVelocity` + +Key implementation details: + +- `computeDenMixLayerDepth` uses a density-threshold criterion referenced to a + fixed near-surface depth (`ReferenceDepth = 10 m`) and linear interpolation to + estimate the crossing depth. +- `computeBuoyGrad` computes horizontal buoyancy gradients at edges and adds + the tilted-coordinate correction using `BruntVaisalaFreqSq`. +- `computeEddyVelocity` forms mixed layer averaged buoyancy and stratification + terms, evaluates frontal-width limits (`LfMin`, `DsMax`), computes a + streamfunction with `shapeFunction()`, and takes its vertical divergence to + produce edge-normal eddy velocity. + +## Computation of mixed layer depth +To compute the mixed layer depth `DenMixedLayerDepth` from the specific volume `SpecVol`, do +```c++ +SubEddies.computeDenMixLayerDepth(SpecVol); +``` + +## Computation of buoyancy gradient +To compute buoyancy gradient `GradBuoyEdgeInterface` from specific volume +`SpecVol`, mean pseudo-thickness on edges `MeanPseudoThickEdge`, mid-layer +`z` coordinate `GeomZMid`, and squared Brunt-Vaisala frequency +`BruntVaisalaFreqSq`, use +```c++ +SubEddies.computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid, BruntVaisalaFreqSq); +``` + +## Computation of eddy velocity +To compute eddy velocity array `EddyVelocity` from squared Brunt-Vaisala +frequency field `BruntVaisalaFreqSq` and mean pseudo-thickness on edges +`MeanPseudoThickEdge`, use +```c++ +SubEddies.computeEddyVelocity(BruntVaisalaFreqSq, MeanPseudoThickEdge); +``` + +## Finalization +To clear the singleton instance, use the static method +```c++ +SubmesoEddies::destroyInstance(); +``` diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index aff5c7a1b8a9..03c888487815 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -53,6 +53,7 @@ userGuide/VertCoord userGuide/PGrad userGuide/Timing userGuide/VerticalMixingCoeff +userGuide/SubmesoEddies userGuide/VertAdv userGuide/Forcing userGuide/SfcCoupling @@ -102,6 +103,7 @@ devGuide/VertCoord devGuide/PGrad devGuide/Timing devGuide/VerticalMixingCoeff +devGuide/SubmesoEddies devGuide/VertAdv devGuide/Forcing devGuide/SfcCoupling diff --git a/components/omega/doc/userGuide/SubmesoEddies.md b/components/omega/doc/userGuide/SubmesoEddies.md new file mode 100644 index 000000000000..0ac2d12c7aed --- /dev/null +++ b/components/omega/doc/userGuide/SubmesoEddies.md @@ -0,0 +1,39 @@ +(omega-user-submeso-eddies)= + +# Submesoscale Eddy Parameterization + +Omega includes an optional submesoscale mixed layer instability (MLI) +parameterization through the `SubmesoEddies` class. When enabled, the model +computes an eddy-induced transport velocity on edges and adds it to the normal +transport velocity used by thickness and tracer advection. + +The current implementation follows the Fox-Kemper et al. (2011) (FK11) closure +described in the +{ref}`omega-design-submesoscale-eddies` design document. + +## Configuration + +Configure the parameterization in the `Submeso` section of the YAML input: + +```yaml +Submeso: + Enable: false + Tau: 172800.0 + Ce: 0.08 + LfMin: 1.0e3 + DsMax: 100.0e3 +``` + +- `Enable`: turns the parameterization on/off. +- `Tau`: MLI timescale parameter (s). +- `Ce`: nondimensional efficiency coefficient. +- `LfMin`: minimum frontal width limiter (m). +- `DsMax`: maximum grid-length limiter used in the closure (m). + +## Diagnostics + +When enabled, the following fields are available in the `Submeso` field group: + +- `DenMixLayerDepth` (m): density-threshold mixed-layer depth. +- `GradBuoyEdgeInterface` (s^-2): buoyancy gradient on edge interfaces. +- `EddyVelocity` (m/s): eddy-induced transport velocity. diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 9a37a0739d98..2db0071ca189 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -4,6 +4,7 @@ #include "Field.h" #include "Logging.h" #include "Pacer.h" +#include "SubmesoEddies.h" #include "Tendencies.h" #include "TimeStepper.h" @@ -30,7 +31,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), SurfTracerRestAux(stripDefault(Name), Mesh, NTracers), TracerAux(stripDefault(Name), Mesh, VCoord, NTracers), - TimeStep(TimeStep) { + TransportAux(stripDefault(Name), Mesh, VCoord), TimeStep(TimeStep) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -46,6 +47,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VelocityDel2Aux.registerFields(GroupName, AuxMeshName); SurfTracerRestAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); + TransportAux.registerFields(GroupName, AuxMeshName); } // Destructor. Unregisters the fields with IOStreams and destroys this auxiliary @@ -57,6 +59,7 @@ AuxiliaryState::~AuxiliaryState() { VelocityDel2Aux.unregisterFields(); SurfTracerRestAux.unregisterFields(); TracerAux.unregisterFields(); + TransportAux.unregisterFields(); FieldGroup::destroy(GroupName); } @@ -98,12 +101,104 @@ void AuxiliaryState::computeMomVertAux(const OceanState *State, // compute geometric height VCoord->computeGeomZHeight(PseudoThickCell, EosInstance->SpecVol); + // compute Brunt-Vaisala freqency squared + EosInstance->computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, PressureMid, + EosInstance->SpecVol); + // compute target thickness VCoord->computeTargetThickness(); Pacer::stop("AuxState:computeMomVertAux", 2); } +// Compute transport velocity for pseudo-thickness and tracers +void AuxiliaryState::computeTransportVelocity(const OceanState *State, + const Array3DReal &TracerArray, + int ThickTimeLevel, + int VelTimeLevel) const { + Pacer::start("AuxState:computeTransportVelocity", 2); + + Array2DReal NormalVel = State->getNormalVelocity(VelTimeLevel); + + const auto &NormalTransportVelocity = TransportAux.NormalTransportVelocity; + + deepCopy(NormalTransportVelocity, NormalVel); + + auto *SubEddies = SubmesoEddies::getInstance(); + + if (SubEddies && SubEddies->Enable) { + + Eos *EosInstance = Eos::getInstance(); + + const auto &MeanPseudoThickEdge = PseudoThicknessAux.MeanPseudoThickEdge; + const auto &SpecVol = EosInstance->SpecVol; + const auto &BVFreqSq = EosInstance->BruntVaisalaFreqSq; + const auto &GeomZMid = VCoord->GeomZMid; + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + + SubEddies->computeDenMixLayerDepth(SpecVol); + SubEddies->computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid, + BVFreqSq); + SubEddies->computeEddyVelocity(BVFreqSq, MeanPseudoThickEdge); + + const auto &EddyVelocity = SubEddies->EddyVelocity; + + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int MinLyrEdgeBot = MinLayerEdgeBot(IEdge); + const int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge); + + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop}, + INNER_LAMBDA(int K) { + NormalTransportVelocity(IEdge, K) += EddyVelocity(IEdge, K); + }); + }); + } + + Pacer::stop("AuxState:computeTransportVelocity", 2); +} + +// Compute the auxiliary variables needed for pseudo-thickness equation +void AuxiliaryState::computePseudoThicknessAux(const OceanState *State, + const Array3DReal &TracerArray, + int ThickTimeLevel, + int VelTimeLevel) const { + + Array2DReal PseudoThick = State->getPseudoThickness(ThickTimeLevel); + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); + OMEGA_SCOPE(LocPseudoThicknessAux, PseudoThicknessAux); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + Pacer::start("Tend:computePseudoThickAux", 2); + + parallelForOuter( + "computePseudoThickAux", {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) { + LocPseudoThicknessAux.computeVarsOnEdge( + IEdge, KChunk, PseudoThick, NormalVelEdge); + }); + }); + + auto *SubEddies = SubmesoEddies::getInstance(); + + if (SubEddies && SubEddies->Enable) { + computeMomVertAux(State, TracerArray, ThickTimeLevel, VelTimeLevel); + } + + computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel); + + Pacer::stop("Tend:computePseudoThickAux", 2); +} + // Compute the auxiliary variables needed for momentum equation void AuxiliaryState::computeMomAux(const OceanState *State, const Array3DReal &TracerArray, @@ -196,6 +291,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State, }); }); + computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel); + parallelForOuter( "edgeAuxState2", {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { @@ -246,15 +343,57 @@ void AuxiliaryState::computeMomAux(const OceanState *State, Pacer::start("AuxState:computeVerticalPseudoVelocity", 2); - const auto &FluxPseudoThickEdge = PseudoThicknessAux.FluxPseudoThickEdge; - VAdv->computeVerticalPseudoVelocity(NormalVelEdge, FluxPseudoThickEdge, - PseudoThickCell, ProjDtSeconds); + const auto &FluxPseudoThickEdge = PseudoThicknessAux.FluxPseudoThickEdge; + const auto &NormalTransportVelocity = TransportAux.NormalTransportVelocity; + VAdv->computeVerticalPseudoVelocity(NormalTransportVelocity, + FluxPseudoThickEdge, PseudoThickCell, + ProjDtSeconds); Pacer::stop("AuxState:computeVerticalPseudoVelocity", 2); Pacer::stop("AuxState:computeMomAux", 1); } +// Compute the auxiliary variables needed for tracer equation +void AuxiliaryState::computeTracerAux(const OceanState *State, + const Array3DReal &TracerArray, + int ThickTimeLevel, + int VelTimeLevel) const { + + OMEGA_SCOPE(LocTracerAux, TracerAux); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + const auto &MeanPseudoThickEdge = PseudoThicknessAux.MeanPseudoThickEdge; + + const int NTracers = Tracers::getNumTracers(); + + Pacer::start("Tend:computeTracerAuxCell", 2); + + auto *SubEddies = SubmesoEddies::getInstance(); + + if (SubEddies && SubEddies->Enable) { + computeMomVertAux(State, TracerArray, ThickTimeLevel, VelTimeLevel); + } + + computeTransportVelocity(State, TracerArray, ThickTimeLevel, VelTimeLevel); + + 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) { + LocTracerAux.computeVarsOnCells( + LTracer, ICell, KChunk, MeanPseudoThickEdge, TracerArray); + }); + }); + Pacer::stop("Tend:computeTracerAuxCell", 2); +} + // Compute the auxiliary variables void AuxiliaryState::computeAll(const OceanState *State, const Array3DReal &TracerArray, diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 91047d6e160e..6cca9ed2db30 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -14,6 +14,7 @@ #include "auxiliaryVars/PseudoThicknessAuxVars.h" #include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" +#include "auxiliaryVars/TransportAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" @@ -42,6 +43,7 @@ class AuxiliaryState { VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; SurfTracerRestAuxVars SurfTracerRestAux; + TransportAuxVars TransportAux; ~AuxiliaryState(); @@ -74,16 +76,31 @@ class AuxiliaryState { /// Exchange halo I4 exchangeHalo(); + // Compute all auxiliary variables needed for pseudo-thickness equation + void computePseudoThicknessAux(const OceanState *State, + const Array3DReal &TracerArray, + int ThickTimeLevel, int VelTimeLevel) const; + // Compute auxiliary variables for vertical dynamics void computeMomVertAux(const OceanState *State, const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel) const; + // Compute transport velocity for pseudo-thickness and tracers + void computeTransportVelocity(const OceanState *State, + const Array3DReal &TracerArray, + int ThickTimeLevel, int VelTimeLevel) const; + // Compute all auxiliary variables needed for momentum equation void computeMomAux(const OceanState *State, const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel, const TimeInterval ProjDt) const; + // Compute all auxiliary variables needed for tracer equation + void computeTracerAux(const OceanState *State, + const Array3DReal &TracerArray, int ThickTimeLevel, + int VelTimeLevel) const; + /// Compute all auxiliary variables based on an ocean state at a given time /// level void computeAll(const OceanState *State, const Array3DReal &TracerArray, diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index 3dfbcf2bd6f7..12c59f8f5362 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -20,6 +20,7 @@ #include "OceanState.h" #include "PGrad.h" #include "SfcCoupling.h" +#include "SubmesoEddies.h" #include "Tendencies.h" #include "TimeMgr.h" #include "TimeStepper.h" @@ -46,6 +47,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time SfcCoupling::clear(); Tracers::clear(); TimeStepper::clear(); + SubmesoEddies::destroyInstance(); PressureGrad::clear(); Eos::destroyInstance(); Tendencies::clear(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 472725f07dc4..179641f2a0c2 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -27,6 +27,7 @@ #include "PGrad.h" #include "Pacer.h" #include "SfcCoupling.h" +#include "SubmesoEddies.h" #include "Tendencies.h" #include "TimeMgr.h" #include "TimeStepper.h" @@ -292,6 +293,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { Eos::init(); PressureGrad::init(); VertMix::init(); + SubmesoEddies::init(); Tendencies::init(); // Validate SurfaceTracerRestoring configuration diff --git a/components/omega/src/ocn/SubmesoEddies.cpp b/components/omega/src/ocn/SubmesoEddies.cpp new file mode 100644 index 000000000000..274851f2b7de --- /dev/null +++ b/components/omega/src/ocn/SubmesoEddies.cpp @@ -0,0 +1,513 @@ +#include "SubmesoEddies.h" +#include "Field.h" +#include "GlobalConstants.h" + +namespace OMEGA { + +static KOKKOS_FUNCTION Real shapeFunction(Real Z, Real H) { + const Real Tmp = (2 * Z) / H + 1; + return Kokkos::max(0._Real, (1 - Tmp * Tmp) * (1 + 5 * Tmp * Tmp / 21)); +} + +/// Instance management +SubmesoEddies *SubmesoEddies::Instance = nullptr; + +/// Get instance of SubmesoEddies +SubmesoEddies *SubmesoEddies::getInstance() { return Instance; } + +/// Destroy instance of SubmesoEddies +void SubmesoEddies::destroyInstance() { + delete Instance; + Instance = nullptr; +} + +/// Initializes the SubmesoEddies class and its options. +/// It assumes that HorzMesh and VertCoord were initialized and +/// initializes the SubmesoEddies class by using the default mesh and vertical +/// coordinate, reading the config file, and setting the parametrization +/// parameters. +void SubmesoEddies::init() { + + if (!Instance) { + Instance = + new SubmesoEddies(HorzMesh::getDefault(), VertCoord::getDefault()); + } + + Error Err; // error code + + /// Retrieve default eos + SubmesoEddies *SubEddies = SubmesoEddies::getInstance(); + + /// Get Submeso group from Omega config + Config *OmegaConfig = Config::getOmegaConfig(); + Config SubmesoConfig("Submeso"); + Err += OmegaConfig->get(SubmesoConfig); + CHECK_ERROR_ABORT(Err, + "SubmesoEddies::init: Submeso group not found in Config"); + + Err += SubmesoConfig.get("Enable", SubEddies->Enable); + CHECK_ERROR_ABORT(Err, + "SubmesoEddies::init: Enable not found in SubmesoConfig"); + + Err += SubmesoConfig.get("Tau", SubEddies->Tau); + CHECK_ERROR_ABORT(Err, + "SubmesoEddies::init: Tau not found in SubmesoConfig"); + + Err += SubmesoConfig.get("Ce", SubEddies->Ce); + CHECK_ERROR_ABORT(Err, "SubmesoEddies::init: Ce not found in SubmesoConfig"); + + Err += SubmesoConfig.get("LfMin", SubEddies->LfMin); + CHECK_ERROR_ABORT(Err, + "SubmesoEddies::init: LfMin not found in SubmesoConfig"); + + Err += SubmesoConfig.get("DsMax", SubEddies->DsMax); + CHECK_ERROR_ABORT(Err, + "SubmesoEddies::init: DsMax not found in SubmesoConfig"); + + // Precompute time scale + SubEddies->computeTimeScale(); + +} // end init + +SubmesoEddies::SubmesoEddies(const HorzMesh *Mesh, const VertCoord *VCoord) + : Mesh(Mesh), VCoord(VCoord), TimeScale("TimeScale", Mesh->NEdgesSize), + DenMixLayerDepth("DenMixLayerDepth", Mesh->NCellsSize), + DenMixLayerIndex("DenMixLayerIndex", Mesh->NCellsSize), + GradBuoyEdgeInterface("GradBuoyEdgeInterface", Mesh->NEdgesSize, + VCoord->NVertLayersP1), + EddyVelocity("EddyVelocity", Mesh->NEdgesSize, VCoord->NVertLayers) { + + // define fields for IO + defineFields(); +} + +void SubmesoEddies::defineFields() { + + // Create a group for the submesoscale eddy parametrization fields + auto SubmesoGroup = FieldGroup::create("Submeso"); + + // Create and add mixed layer depth field + { + int NDims = 1; + std::vector DimNames(NDims); + DimNames[0] = "NCells"; + + auto DenMixLayerDepthField = + Field::create(DenMixLayerDepth.label(), // Field name + "Mixed Layer Depth", // Long Name + "m", // Units + "", // CF-ish Name + 0.0, // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Number of dimensions + DimNames // Dimension names + ); + + DenMixLayerDepthField->attachData(DenMixLayerDepth); + + SubmesoGroup->addField(DenMixLayerDepth.label()); + } + + // Create and add buoyancy gradient field + { + int NDims = 2; + std::vector DimNames(NDims); + DimNames[0] = "NEdges"; + DimNames[1] = "NVertLayersP1"; + + auto BuoyancyGradientInterfaceField = + Field::create(GradBuoyEdgeInterface.label(), // Field name + "Buoyancy Gradient", // Long Name + "1/s^2", // Units + "", // CF-ish Name + std::numeric_limits::lowest(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Number of dimensions + DimNames // Dimension names + ); + + BuoyancyGradientInterfaceField->attachData( + GradBuoyEdgeInterface); + + SubmesoGroup->addField(GradBuoyEdgeInterface.label()); + } + + // Create and add eddy velocity field + { + int NDims = 2; + std::vector DimNames(NDims); + DimNames[0] = "NEdges"; + DimNames[1] = "NVertLayers"; + + auto EddyVelocityField = + Field::create(EddyVelocity.label(), // Field name + "Eddy Velocity", // Long Name + "m/s", // Units + "", // CF-ish Name + std::numeric_limits::lowest(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Number of dimensions + DimNames // Dimension names + ); + + EddyVelocityField->attachData(EddyVelocity); + + SubmesoGroup->addField(EddyVelocity.label()); + } +} + +void SubmesoEddies::computeTimeScale() { + OMEGA_SCOPE(TimeScale, this->TimeScale); + OMEGA_SCOPE(Tau, this->Tau); + + const auto &FEdge = Mesh->FEdge; + + parallelFor( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + TimeScale(IEdge) = + Kokkos::sqrt(FEdge(IEdge) * FEdge(IEdge) + 1._Real / (Tau * Tau)); + }); +} + +void SubmesoEddies::computeDenMixLayerDepth(const Array2DReal &SpecVol) { + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + const auto &GeomZInterface = VCoord->GeomZInterface; + const auto &GeomZMid = VCoord->GeomZMid; + + OMEGA_SCOPE(ReferenceDepth, this->ReferenceDepth); + OMEGA_SCOPE(DenThreshold, this->DenThreshold); + OMEGA_SCOPE(DenMixLayerDepth, this->DenMixLayerDepth); + OMEGA_SCOPE(DenMixLayerIndex, this->DenMixLayerIndex); + + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const Real SSH = GeomZInterface(ICell, MinLayerCell(ICell)); + + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + + // Find first interface where depth >= reference depth + int KRef; + parallelSearchInner( + Team, Range{KMin + 1, KMax}, + INNER_LAMBDA(int K) { + const Real Depth = SSH - GeomZInterface(ICell, K); + return Depth >= ReferenceDepth; + }, + KRef); + + // Not found, setting to KMax + if (KRef == -1) { + KRef = KMax; + } + + const int KRefM1 = Kokkos::max(KRef - 1, MinLayerCell(ICell)); + + const Real DepthKRef = SSH - GeomZMid(ICell, KRef); + const Real DepthKRefM1 = SSH - GeomZMid(ICell, KRefM1); + + const Real ReferenceSpecVol = + linearInterp(ReferenceDepth, SpecVol(ICell, KRef), DepthKRef, + SpecVol(ICell, KRefM1), DepthKRefM1); + + // Start searching from reference level - 1 + int KDen; + parallelSearchInner( + Team, Range{KRefM1, KMax}, + INNER_LAMBDA(int K) { + return (ReferenceSpecVol / SpecVol(ICell, K) - 1) >= + DenThreshold * ReferenceSpecVol; + }, + KDen); + + // Not found. Setting to the depth of the deepest layer + if (KDen == -1) { + DenMixLayerIndex(ICell) = KMax; + DenMixLayerDepth(ICell) = SSH - GeomZMid(ICell, KMax); + } else { // Found + const int KDenM1 = Kokkos::max(KDen - 1, MinLayerCell(ICell)); + + const Real DepthKDen = SSH - GeomZMid(ICell, KDen); + const Real DepthKDenM1 = SSH - GeomZMid(ICell, KDenM1); + + const Real FactorKDen = + ReferenceSpecVol / SpecVol(ICell, KDen) - 1; + const Real FactorKDenM1 = + ReferenceSpecVol / SpecVol(ICell, KDenM1) - 1; + + const Real MixedLayerDepth = + linearInterp(DenThreshold * ReferenceSpecVol, DepthKDen, + FactorKDen, DepthKDenM1, FactorKDenM1); + + DenMixLayerIndex(ICell) = KDen; + DenMixLayerDepth(ICell) = MixedLayerDepth; + } + }); +} + +void SubmesoEddies::computeBuoyGrad(const Array2DReal &SpecVol, + const Array2DReal &MeanPseudoThickEdge, + const Array2DReal &GeomZMid, + const Array2DReal &BruntVaisalaFreqSq) { + OMEGA_SCOPE(GradBuoyEdgeInterface, this->GradBuoyEdgeInterface); + + const auto &DcEdge = Mesh->DcEdge; + const auto &CellsOnEdge = Mesh->CellsOnEdge; + + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + const auto NVertLayers = VCoord->NVertLayers; + const auto NVertLayersP1 = VCoord->NVertLayersP1; + + parallelForOuter( + LaunchConfig({Mesh->NEdgesAll}, + TeamScratch(3 * NVertLayers + 3 * NVertLayersP1)), + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + int MinLyrEdgeBot = MinLayerEdgeBot(IEdge); + int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge); + + ScratchArray1DReal GradBuoyEdge(teamScratch(Team), NVertLayers); + ScratchArray1DReal GradGeomZMidEdge(teamScratch(Team), NVertLayers); + ScratchArray1DReal BVFSqEdge(teamScratch(Team), NVertLayersP1); + ScratchArray1DReal SpecVolEdge(teamScratch(Team), NVertLayers); + + // Horizontal interpolations and gradients + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop}, INNER_LAMBDA(int K) { + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + + const Real InvDcEdge = 1._Real / DcEdge(IEdge); + + // interpolate Brunt-Vaisala freq to edges + BVFSqEdge(K) = 0.5_Real * (BruntVaisalaFreqSq(JCell1, K) + + BruntVaisalaFreqSq(JCell0, K)); + + // interpolate SpecVol to Edges + SpecVolEdge(K) = + 0.5_Real * (SpecVol(JCell1, K) + SpecVol(JCell0, K)); + + // Compute grad of GeomZMid + GradGeomZMidEdge(K) = + InvDcEdge * (GeomZMid(JCell1, K) - GeomZMid(JCell0, K)); + + // Compute grad of buoyancy + GradBuoyEdge(K) = -Gravity * RhoSw * InvDcEdge * + (SpecVol(JCell1, K) - SpecVol(JCell0, K)); + }); + + // Interpolate Brunt-Vaisala freq to edges at the bottom interface + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + const int K = MaxLyrEdgeTop + 1; + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + BVFSqEdge(K) = 0.5_Real * (BruntVaisalaFreqSq(JCell1, K) + + BruntVaisalaFreqSq(JCell0, K)); + }); + + teamBarrier(Team); + + ScratchArray1DReal GradGeomZMidEdgeInterface(teamScratch(Team), + NVertLayersP1); + ScratchArray1DReal SpecVolEdgeInterface(teamScratch(Team), + NVertLayersP1); + + // Vertical interpolations + + // This interpolation can only be carried out on non-boundary edges + if (MaxLyrEdgeTop >= MinLyrEdgeBot) { + + parallelForInner( + Team, Range{MinLyrEdgeBot + 1, MaxLyrEdgeTop}, + INNER_LAMBDA(int K) { + const Real PseudoThickKm1 = + MeanPseudoThickEdge(IEdge, K - 1); + const Real PseudoThickK = MeanPseudoThickEdge(IEdge, K); + + const Real CoeffKm1 = + PseudoThickKm1 / (PseudoThickKm1 + PseudoThickK); + const Real CoeffK = + PseudoThickK / (PseudoThickKm1 + PseudoThickK); + + // interpolate GeomZMid gradient to interfaces + GradGeomZMidEdgeInterface(K) = + CoeffKm1 * GradGeomZMidEdge(K - 1) + + CoeffK * GradGeomZMidEdge(K); + + // interpolate SpecVol to interfaces + SpecVolEdgeInterface(K) = + CoeffKm1 * SpecVolEdge(K - 1) + CoeffK * SpecVolEdge(K); + + // interpolate GradBuoyEdge to interfaces + GradBuoyEdgeInterface(IEdge, K) = + CoeffKm1 * GradBuoyEdge(K - 1) + + CoeffK * GradBuoyEdge(K); + }); + + teamBarrier(Team); + + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + SpecVolEdgeInterface(MinLyrEdgeBot) = + SpecVolEdge(MinLyrEdgeBot); + SpecVolEdgeInterface(MaxLyrEdgeTop + 1) = + SpecVolEdge(MaxLyrEdgeTop); + + GradGeomZMidEdgeInterface(MinLyrEdgeBot) = + GradGeomZMidEdge(MinLyrEdgeBot); + GradGeomZMidEdgeInterface(MaxLyrEdgeTop + 1) = + GradGeomZMidEdge(MaxLyrEdgeTop); + + GradBuoyEdgeInterface(IEdge, MinLyrEdgeBot) = + GradBuoyEdge(MinLyrEdgeBot); + GradBuoyEdgeInterface(IEdge, MaxLyrEdgeTop + 1) = + GradBuoyEdge(MaxLyrEdgeTop); + }); + + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop + 1}, + INNER_LAMBDA(int K) { + GradBuoyEdgeInterface(IEdge, K) += + GradGeomZMidEdgeInterface(K) * RhoSw * + SpecVolEdgeInterface(K) * BVFSqEdge(K); + }); + } + }); +} + +void SubmesoEddies::computeEddyVelocity( + const Array2DReal &BruntVaisalaFreqSq, + const Array2DReal &MeanPseudoThickEdge) { + + OMEGA_SCOPE(GradBuoyEdgeInterface, this->GradBuoyEdgeInterface); + OMEGA_SCOPE(DenMixLayerIndex, this->DenMixLayerIndex); + OMEGA_SCOPE(DenMixLayerDepth, this->DenMixLayerDepth); + OMEGA_SCOPE(TimeScale, this->TimeScale); + OMEGA_SCOPE(LfMin, this->LfMin); + OMEGA_SCOPE(DsMax, this->DsMax); + OMEGA_SCOPE(Ce, this->Ce); + OMEGA_SCOPE(EddyVelocity, this->EddyVelocity); + + const auto &DcEdge = Mesh->DcEdge; + const auto &CellsOnEdge = Mesh->CellsOnEdge; + + const auto &GeomZInterface = VCoord->GeomZInterface; + + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + const auto NVertLayersP1 = VCoord->NVertLayersP1; + + // Replace with global constant when added + const Real Tiny = 1e-12_Real; + + parallelForOuter( + LaunchConfig({Mesh->NEdgesAll}, TeamScratch(NVertLayersP1)), + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int MinLyrEdgeBot = MinLayerEdgeBot(IEdge); + const int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge); + + if (MaxLyrEdgeTop >= MinLyrEdgeBot) { + + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + + const int MinLayerCell0 = MinLayerCell(JCell0); + const int MinLayerCell1 = MinLayerCell(JCell1); + + Real PseudoThickML; + Real GradBuoyML; + Real BVFreqML; + + const I4 IndexMLEdge = Kokkos::min(DenMixLayerIndex(JCell0), + DenMixLayerIndex(JCell1)); + + // compute mixed layer averages of buoyancy gradient and + // Brunt-Vaisala frequency + parallelReduceInner( + Team, Range{MinLyrEdgeBot, IndexMLEdge}, + INNER_LAMBDA(int K, Real &AccumThick, Real &AccumGradBuoy, + Real &AccumBVFreq) { + const Real PseudoThickKm1 = + ((K - 1) >= MinLyrEdgeBot + ? MeanPseudoThickEdge(IEdge, K - 1) + : 0); + const Real PseudoThickK = MeanPseudoThickEdge(IEdge, K); + const Real PseudoThickAvg = + 0.5_Real * (PseudoThickKm1 + PseudoThickK); + + const Real BVFSq0 = BruntVaisalaFreqSq(JCell0, K); + const Real BVFSq1 = BruntVaisalaFreqSq(JCell1, K); + const Real GradBuoy = GradBuoyEdgeInterface(IEdge, K); + + const Real BVFEdge = + Kokkos::sqrt(0.5_Real * (Kokkos::max(Tiny, BVFSq0) + + Kokkos::max(Tiny, BVFSq1))); + + AccumThick += PseudoThickAvg; + AccumGradBuoy += PseudoThickAvg * GradBuoy; + AccumBVFreq += PseudoThickAvg * BVFEdge; + }, + PseudoThickML, GradBuoyML, BVFreqML); + + GradBuoyML /= PseudoThickML; + BVFreqML /= PseudoThickML; + + // compute stream function + ScratchArray1DReal StreamFunction(teamScratch(Team), + NVertLayersP1); + parallelForInner( + Team, NVertLayersP1, + INNER_LAMBDA(int K) { StreamFunction(K) = 0; }); + + const Real MLDepthEdge = Kokkos::min(DenMixLayerDepth(JCell0), + DenMixLayerDepth(JCell1)); + + const Real TScale = TimeScale(IEdge); + const Real Ds = Kokkos::min(DcEdge(IEdge), DsMax); + + const Real Lf1 = + Kokkos::abs(GradBuoyML) * MLDepthEdge / (TScale * TScale); + const Real Lf2 = BVFreqML * MLDepthEdge / TScale; + const Real Lf = Kokkos::max(LfMin, Kokkos::max(Lf1, Lf2)); + + const Real Factor = + Ce * Ds / Lf * MLDepthEdge * MLDepthEdge * GradBuoyML / TScale; + + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop + 1}, + INNER_LAMBDA(int K) { + const Real ZEdge = + 0.5_Real * (GeomZInterface(JCell0, K) - + GeomZInterface(JCell0, MinLayerCell0) + + GeomZInterface(JCell1, K) - + GeomZInterface(JCell1, MinLayerCell1)); + + const Real Mu = shapeFunction(ZEdge, MLDepthEdge); + + StreamFunction(K) = Factor * Mu; + }); + + teamBarrier(Team); + + // compute eddy velocity + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop}, + INNER_LAMBDA(int K) { + const Real DZ = 0.5_Real * (GeomZInterface(JCell0, K) - + GeomZInterface(JCell0, K + 1) + + GeomZInterface(JCell1, K) - + GeomZInterface(JCell1, K + 1)); + EddyVelocity(IEdge, K) = + -(StreamFunction(K) - StreamFunction(K + 1)) / DZ; + }); + } + }); +} + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/SubmesoEddies.h b/components/omega/src/ocn/SubmesoEddies.h new file mode 100644 index 000000000000..870b022dc3cb --- /dev/null +++ b/components/omega/src/ocn/SubmesoEddies.h @@ -0,0 +1,124 @@ +#ifndef OMEGA_SUBMESOEDDIES_H +#define OMEGA_SUBMESOEDDIES_H + +#include "HorzMesh.h" +#include "VertCoord.h" + +namespace OMEGA { + +// Generic linear interpolation routine. This should be in OmegaMath.h or +// something like that once it exists. +KOKKOS_INLINE_FUNCTION Real linearInterp(Real x, Real y1, Real x1, Real y2, + Real x2) { + const Real A = (y1 - y2) / (x1 - x2); + const Real B = y1 - A * x1; + return (x1 == x2) ? y1 : A * x + B; +} + +// A class for the submesoscale eddy parametrization. +// It groups the variables related to this parametrization +// and provides methods to compute mixed layer depth, buoyancy gradient, +// and eddy velocity. +class SubmesoEddies { + public: + // Public methods + + // Initialize SubmesoEddies from config, default mesh, and default vertical + // coordinate + static void init(); + + // Get instance of SubmesoEddies + static SubmesoEddies *getInstance(); + + // Destroy instance (frees Kokkos views) + static void destroyInstance(); + + // Compute time scale array from time scale constant and Coriolis parameter + void computeTimeScale(); + + // Compute mixed layer index and depth based on the density difference + // criterion + void computeDenMixLayerDepth(const Array2DReal &SpecVol); + + // Compute buoyancy gradient + void computeBuoyGrad(const Array2DReal &SpecVol, + const Array2DReal &MeanPseudoThickEdge, + const Array2DReal &GeomZMid, + const Array2DReal &BruntVaisalaFreqSq); + + // Compute eddy velocity + void computeEddyVelocity(const Array2DReal &BruntVaisalaFreqSq, + const Array2DReal &MeanPseudoThickEdge); + + // Public member variables + + // Enable the parametrization + bool Enable; + + // Parametrization constants + + // Minimum width of submesoscale fronts + Real LfMin; + + // Efficiency coefficient + Real Ce; + + // Maximum edge length + Real DsMax; + + // Time scale constant + Real Tau; + + // Reference depth + Real ReferenceDepth = 10; + + // Density threshold for determining the mixed layer depth + Real DenThreshold = 0.03; + + // Mixed layer index + Array1DI4 DenMixLayerIndex; + + // Mixed layer depth + Array1DReal DenMixLayerDepth; + + // Buoyancy gradient + Array2DReal GradBuoyEdgeInterface; + + // Eddy Velocity + Array2DReal EddyVelocity; + + // Time scale array + Array1DReal TimeScale; + + private: + // Private methods + + // Constructor + SubmesoEddies(const HorzMesh *Mesh, const VertCoord *VCoord); + + // Destructor + ~SubmesoEddies() = default; + + // Delete copy and move constructors and assignment operators + SubmesoEddies(const SubmesoEddies &) = delete; + SubmesoEddies &operator=(const SubmesoEddies &) = delete; + SubmesoEddies(SubmesoEddies &&) = delete; + SubmesoEddies &operator=(SubmesoEddies &&) = delete; + + // Define fields and metadata + void defineFields(); + + // Private member variables + + // Instance pointer + static SubmesoEddies *Instance; + + // Pointer to horizontal mesh + const HorzMesh *Mesh; + + // Pointer to vertical coordinate + const VertCoord *VCoord; +}; + +} // namespace OMEGA +#endif diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 5bb131f32ab6..2003f69da845 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -566,8 +566,6 @@ void Tendencies::computePseudoThicknessTendenciesOnly( OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); - Pacer::start("Tend:computePseudoThicknessTendenciesOnly", 1); parallelForOuter( @@ -581,6 +579,8 @@ void Tendencies::computePseudoThicknessTendenciesOnly( }); // Compute pseudo-thickness flux divergence + const Array2DReal &NormalTransportVelocity = + AuxState->TransportAux.NormalTransportVelocity; const Array2DReal &ThickFluxEdge = AuxState->PseudoThicknessAux.FluxPseudoThickEdge; @@ -595,7 +595,7 @@ void Tendencies::computePseudoThicknessTendenciesOnly( parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { LocThicknessFluxDiv(LocPseudoThicknessTend, ICell, KChunk, - ThickFluxEdge, NormalVelEdge); + ThickFluxEdge, NormalTransportVelocity); }); }); Pacer::stop("Tend:thicknessFluxDiv", 2); @@ -843,7 +843,8 @@ void Tendencies::computeTracerTendenciesOnly( }); // compute tracer horizotal advection - Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); + const Array2DReal &NormalTransportVelocity = + AuxState->TransportAux.NormalTransportVelocity; const Array2DReal &FluxPseudoThickEdge = AuxState->PseudoThicknessAux.FluxPseudoThickEdge; if (LocTracerHorzAdv.Enabled) { @@ -857,7 +858,8 @@ void Tendencies::computeTracerTendenciesOnly( parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { LocTracerHorzAdv(L, IEdge, KChunk, TracerArray, - FluxPseudoThickEdge, NormalVelEdge); + FluxPseudoThickEdge, + NormalTransportVelocity); }); }); parallelForOuter( @@ -952,36 +954,15 @@ void Tendencies::computeTracerTendenciesOnly( void Tendencies::computePseudoThicknessTendencies( const OceanState *State, ///< [in] State variables const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + const Array3DReal &TracerArray, ///< [in] Tracer array int ThickTimeLevel, ///< [in] Time level int VelTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { - // only need PseudoThicknessAux on edge - Array2DReal PseudoThick = State->getPseudoThickness(ThickTimeLevel); - Array2DReal NormVel = State->getNormalVelocity(VelTimeLevel); - OMEGA_SCOPE(PseudoThicknessAux, AuxState->PseudoThicknessAux); - OMEGA_SCOPE(PseudoThickCell, PseudoThick); - OMEGA_SCOPE(NormalVelEdge, NormVel); - OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); - OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); - Pacer::start("Tend:computePseudoThicknessTendencies", 1); - Pacer::start("Tend:computePseudoThickAux", 2); - parallelForOuter( - "computePseudoThickAux", {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) { - PseudoThicknessAux.computeVarsOnEdge( - IEdge, KChunk, PseudoThickCell, NormalVelEdge); - }); - }); - Pacer::stop("Tend:computePseudoThickAux", 2); + AuxState->computePseudoThicknessAux(State, TracerArray, ThickTimeLevel, + VelTimeLevel); computePseudoThicknessTendenciesOnly(State, AuxState, ThickTimeLevel, VelTimeLevel, Time); @@ -1018,33 +999,10 @@ void Tendencies::computeTracerTendencies( int VelTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { - Array2DReal PseudoThickCell = State->getPseudoThickness(ThickTimeLevel); - Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); - OMEGA_SCOPE(TracerAux, AuxState->TracerAux); - 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); - const auto &MeanPseudoThickEdge = - AuxState->PseudoThicknessAux.MeanPseudoThickEdge; - Pacer::start("Tend:computeTracerAuxCell", 2); - 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, - MeanPseudoThickEdge, TracerArray); - }); - }); - Pacer::stop("Tend:computeTracerAuxCell", 2); + AuxState->computeTracerAux(State, TracerArray, ThickTimeLevel, VelTimeLevel); computeTracerTendenciesOnly(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, Time); diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index c277547c9e62..2d541658b75a 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -84,6 +84,7 @@ class Tendencies { // Methods to compute tendency groups void computePseudoThicknessTendencies(const OceanState *State, const AuxiliaryState *AuxState, + const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel, TimeInstant Time); void computeVelocityTendencies(const OceanState *State, diff --git a/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.cpp new file mode 100644 index 000000000000..07efca066fc6 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.cpp @@ -0,0 +1,57 @@ +#include "TransportAuxVars.h" +#include "DataTypes.h" +#include "Field.h" + +#include + +namespace OMEGA { + +TransportAuxVars::TransportAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, + const VertCoord *VCoord) + : NormalTransportVelocity("NormalTransportVelocity" + AuxStateSuffix, + Mesh->NEdgesSize, VCoord->NVertLayers) {} + +void TransportAuxVars::registerFields( + const std::string &AuxGroupName, // name of Auxiliary field group + const std::string &MeshName // name of horizontal mesh +) const { + + // Create/define fields + int NDims = 2; + std::vector DimNames(NDims); + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + DimNames[0] = "NEdges" + DimSuffix; + DimNames[1] = "NVertLayers"; + + auto NormalTransportVelocityField = + Field::create(NormalTransportVelocity.label(), // field name + "horizontal velocity used to transport pseudo-thickness " + "and tracers", // long Name or description + "m/s", // units + "", // CF standard Name + std::numeric_limits::lowest(), // min valid value + std::numeric_limits::max(), // max valid value + NDims, // number of dimensions + DimNames // dimension names + ); + + // Add fields to FieldGroup + FieldGroup::addFieldToGroup(NormalTransportVelocity.label(), AuxGroupName); + + // Attach data + NormalTransportVelocityField->attachData( + NormalTransportVelocity); +} + +void TransportAuxVars::unregisterFields() const { + Field::destroy(NormalTransportVelocity.label()); +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.h b/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.h new file mode 100644 index 000000000000..63b45077e276 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/TransportAuxVars.h @@ -0,0 +1,25 @@ +#ifndef OMEGA_AUX_TRANSPORT_H +#define OMEGA_AUX_TRANSPORT_H + +#include "DataTypes.h" +#include "HorzMesh.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +class TransportAuxVars { + public: + Array2DReal NormalTransportVelocity; + + TransportAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, + const VertCoord *VCoord); + + void registerFields(const std::string &AuxGroupName, + const std::string &MeshName) const; + void unregisterFields() const; +}; + +} // namespace OMEGA +#endif diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index e1bb57802b70..d5ddab5a9da7 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -73,8 +73,8 @@ void ForwardBackwardStepper::doStep( Pacer::stop("ForwardBackward:velHaloExch", 3); // R_h^{n} = RHS_h(u^{n+1}, h^{n}, t^{n}) - Tend->computePseudoThicknessTendencies(State, AuxState, ThickCurLevel, - VelNextLevel, SimTime); + Tend->computePseudoThicknessTendencies(State, AuxState, CurTracerArray, + ThickCurLevel, VelNextLevel, SimTime); // h^{n+1} = h^{n} + R_h^{n} updateThicknessByTend(State, ThickNextLevel, State, ThickCurLevel, TimeStep); diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 20fe917bdd48..3ca8a75f66e6 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -600,6 +600,7 @@ add_omega_test( ocn/SfcCouplingTest.cpp "-n;8" ) + ################## # Analysis tests ################## @@ -617,3 +618,14 @@ add_omega_test( analysis/AnalysisSystemTest.cpp "-n;8" ) + +########################## +# Submeso Eddies test +########################## + +add_omega_test( + SUBMESOEDDIES_TEST + testSubmesoEddies.exe + ocn/SubmesoEddiesTest.cpp + "-n;2" +) diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index f9bd3633bfae..51a1a0389414 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -63,23 +63,25 @@ int initState() { int NTracers = Tracers::getNumTracers(); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThickness(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.pseudoThickness(X, Y); + }, PseudoThickCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, - VCoord->MaxLayerCell); + VCoord->MaxLayerCell, nullptr); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.tracer(X, Y); }, TracerArray, Geom, Mesh, OnCell, VCoord->MinLayerCell, - VCoord->MaxLayerCell); + VCoord->MaxLayerCell, nullptr); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real Lon, Real Lat) { VecField[0] = Setup.velocityX(Lon, Lat); VecField[1] = Setup.velocityY(Lon, Lat); }, NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, - VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, ExchangeHalos::Yes, - CartProjection::No); + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, nullptr, + ExchangeHalos::Yes, CartProjection::No); return Err; } diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 6260519f38ed..ea08778acc4d 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -299,11 +299,13 @@ int initState(const Array2DReal &PseudoThickCell, TestSetup Setup; Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThickness(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.pseudoThickness(X, Y); + }, PseudoThickCell, Geom, Mesh, OnCell); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.velocityX(X, Y); VecField[1] = Setup.velocityY(X, Y); }, @@ -313,7 +315,9 @@ int initState(const Array2DReal &PseudoThickCell, const auto &FVertex = Mesh->FVertex; Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.planetaryVorticity(X, Y); }, + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { + return Setup.planetaryVorticity(X, Y); + }, FVertex, Geom, Mesh, OnVertex); return Err; @@ -332,13 +336,17 @@ int testKineticAuxVars(const Array2DReal &PseudoThicknessCell, Array2DReal ExactKineticEnergyCell("ExactKineticEnergyCell", Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.kineticEnergy(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.kineticEnergy(X, Y); + }, ExactKineticEnergyCell, Geom, Mesh, OnCell, ExchangeHalos::No); Array2DReal ExactVelocityDivCell("ExactVelocityDivCell", Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.divergence(X, Y); + }, ExactVelocityDivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical result @@ -385,7 +393,9 @@ int testPseudoThicknessAuxVars(const Array2DReal &PseudoThickCell, Array2DReal ExactThickEdge("ExactThickEdge", Mesh->NEdgesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThickness(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { + return Setup.pseudoThickness(X, Y); + }, ExactThickEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); // Compute numerical result @@ -437,13 +447,15 @@ int testVorticityAuxVars(const Array2DReal &PseudoThickCell, Array2DReal ExactRelVortVertex("ExactRelVortVertex", Mesh->NVerticesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.relativeVorticity(X, Y); }, + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { + return Setup.relativeVorticity(X, Y); + }, ExactRelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); Array2DReal ExactNormRelVortVertex("ExactNormRelVortVertex", Mesh->NVerticesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { return Setup.normalizedRelativeVorticity(X, Y); }, ExactNormRelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); @@ -451,7 +463,7 @@ int testVorticityAuxVars(const Array2DReal &PseudoThickCell, Array2DReal ExactNormPlanetVortVertex("ExactNormPlanetVortVertex", Mesh->NVerticesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { return Setup.normalizedPlanetaryVorticity(X, Y); }, ExactNormPlanetVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); @@ -496,7 +508,7 @@ int testVorticityAuxVars(const Array2DReal &PseudoThickCell, Array2DReal ExactNormRelVortEdge("ExactNormRelVortEdge", Mesh->NEdgesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return Setup.normalizedRelativeVorticity(X, Y); }, ExactNormRelVortEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); @@ -504,7 +516,7 @@ int testVorticityAuxVars(const Array2DReal &PseudoThickCell, Array2DReal ExactNormPlanetVortEdge("ExactNormPlanetVortEdge", Mesh->NEdgesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return Setup.normalizedPlanetaryVorticity(X, Y); }, ExactNormPlanetVortEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); @@ -554,20 +566,24 @@ int testVelocityDel2AuxVars(Real RTol) { Array2DReal ExactVelocityDivCell("ExactVelocityDivCell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.divergence(X, Y); + }, ExactVelocityDivCell, Geom, Mesh, OnCell); Array2DReal ExactRelVortVertex("ExactRelVortVertex", Mesh->NVerticesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.relativeVorticity(X, Y); }, + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { + return Setup.relativeVorticity(X, Y); + }, ExactRelVortVertex, Geom, Mesh, OnVertex); // Compute exact Del2 Array2DReal ExactDel2Edge("ExactDel2Edge", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.velocityDel2X(X, Y); VecField[1] = Setup.velocityDel2Y(X, Y); }, @@ -596,7 +612,9 @@ int testVelocityDel2AuxVars(Real RTol) { Array2DReal ExactDel2DivCell("ExactDel2DivCell", Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.velocityDel2Div(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.velocityDel2Div(X, Y); + }, ExactDel2DivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical Del2Div @@ -620,7 +638,9 @@ int testVelocityDel2AuxVars(Real RTol) { Array2DReal ExactDel2RelVortVertex("ExactDel2RelVortVertex", Mesh->NVerticesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.velocityDel2Curl(X, Y); }, + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { + return Setup.velocityDel2Curl(X, Y); + }, ExactDel2RelVortVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); // Compute numerical Del2RelVort @@ -663,13 +683,15 @@ int testTracerAuxVars(const Array2DReal &PseudoThickCell, Array3DReal TracersOnCell("TracersOnCell", NTracers, Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.tracer(X, Y); }, TracersOnCell, Geom, Mesh, OnCell); Array2DReal PseudoThickEdge("PseudoThickEdge", Mesh->NEdgesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThickness(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { + return Setup.pseudoThickness(X, Y); + }, PseudoThickEdge, Geom, Mesh, OnEdge); // Compute exact Del2TracerCell @@ -677,7 +699,9 @@ int testTracerAuxVars(const Array2DReal &PseudoThickCell, Array3DReal ExactDel2TrCell("ExactDel2TrCell", NTracers, Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.del2Tracer(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.del2Tracer(X, Y); + }, ExactDel2TrCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical Del2TracerCell diff --git a/components/omega/test/ocn/ForcingTest.cpp b/components/omega/test/ocn/ForcingTest.cpp index 2ac027553f05..e477862ee8af 100644 --- a/components/omega/test/ocn/ForcingTest.cpp +++ b/components/omega/test/ocn/ForcingTest.cpp @@ -84,7 +84,7 @@ int testSfcStressForcingVars(Real RTol) { Array1DReal ExactNormalStressEdge("ExactNormalStressEdge", Mesh->NEdgesOwned); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.sfcStressX(X, Y); VecField[1] = Setup.sfcStressY(X, Y); }, @@ -96,11 +96,15 @@ int testSfcStressForcingVars(Real RTol) { // Set inputs Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.sfcStressX(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.sfcStressX(X, Y); + }, SfcStressForcing.ZonalStressCell, Geom, Mesh, OnCell); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.sfcStressY(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.sfcStressY(X, Y); + }, SfcStressForcing.MeridStressCell, Geom, Mesh, OnCell); // Compute numerical result diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index f4af88fd369c..7595cee5e058 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -199,7 +199,7 @@ int testDivergence(Real RTol) { // Prepare operator input Array2DReal VecEdge("VecEdge", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, @@ -208,7 +208,9 @@ int testDivergence(Real RTol) { // Compute exact result Array2DReal ExactDivCell("ExactDivCell", Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.exactDivVec(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.exactDivVec(X, Y); + }, ExactDivCell, Geom, Mesh, OnCell, ExchangeHalos::No); // Compute numerical result @@ -243,7 +245,7 @@ int testGradient(Real RTol) { // Prepare operator input Array2DReal ScalarCell("ScalarCell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real Coord1, Real Coord2) { + KOKKOS_LAMBDA(int ICell, Real Coord1, Real Coord2) { return Setup.exactScalar(Coord1, Coord2); }, ScalarCell, Geom, Mesh, OnCell); @@ -251,7 +253,7 @@ int testGradient(Real RTol) { // Compute exact result Array2DReal ExactGradEdge("ExactGradEdge", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.exactGradScalarX(X, Y); VecField[1] = Setup.exactGradScalarY(X, Y); }, @@ -288,7 +290,7 @@ int testCurl(Real RTol) { // Prepare operator input Array2DReal VecEdge("VecEdge", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, @@ -298,7 +300,9 @@ int testCurl(Real RTol) { Array2DReal ExactCurlVertex("ExactCurlVertex", Mesh->NVerticesOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.exactCurlVec(X, Y); }, + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { + return Setup.exactCurlVec(X, Y); + }, ExactCurlVertex, Geom, Mesh, OnVertex, ExchangeHalos::No); // Compute numerical result @@ -335,7 +339,7 @@ int testTangentRecon(Real RTol) { // Prepare operator input Array2DReal VecEdge("VecEdge", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, @@ -345,7 +349,7 @@ int testTangentRecon(Real RTol) { Array2DReal ExactReconEdge("ExactReconEdge", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.exactVecX(X, Y); VecField[1] = Setup.exactVecY(X, Y); }, @@ -398,7 +402,7 @@ int testVectorRecon(Real RTol) { // layer, so we use rank-1 arrays here) Array1DReal VecEdge("VecEdge", Mesh->NEdgesSize); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real Lon, Real Lat) { VecField[0] = Setup.exactVecX(Lon, Lat); VecField[1] = Setup.exactVecY(Lon, Lat); }, @@ -407,7 +411,7 @@ int testVectorRecon(Real RTol) { // Compute exact magnitude of the vector field at cell centers Array1DReal ExactMagCell("ExactMagCell", Mesh->NCellsOwned); Err += setScalar( - KOKKOS_LAMBDA(Real Lon, Real Lat) { + KOKKOS_LAMBDA(int ICell, Real Lon, Real Lat) { return vecMagnitude(Setup.exactVecX(Lon, Lat), Setup.exactVecY(Lon, Lat)); }, @@ -451,7 +455,7 @@ int testInterpCellToEdge(Real RTol) { // Prepare operator input Array1DReal ScalarCell("ScalarCell", Mesh->NCellsSize); Err += setScalar( - KOKKOS_LAMBDA(Real Coord1, Real Coord2) { + KOKKOS_LAMBDA(int ICell, Real Coord1, Real Coord2) { return Setup.exactScalar(Coord1, Coord2); }, ScalarCell, Geom, Mesh, OnCell); @@ -459,7 +463,7 @@ int testInterpCellToEdge(Real RTol) { // Compute exact result Array1DReal ExactScalarEdge("ExactScalarEdge", Mesh->NEdgesOwned); Err += setScalar( - KOKKOS_LAMBDA(Real Coord1, Real Coord2) { + KOKKOS_LAMBDA(int IEdge, Real Coord1, Real Coord2) { return Setup.exactScalar(Coord1, Coord2); }, ExactScalarEdge, Geom, Mesh, OnEdge, ExchangeHalos::No); diff --git a/components/omega/test/ocn/OceanTestCommon.h b/components/omega/test/ocn/OceanTestCommon.h index 63f03d44a9db..ad67a89d4d88 100644 --- a/components/omega/test/ocn/OceanTestCommon.h +++ b/components/omega/test/ocn/OceanTestCommon.h @@ -8,6 +8,8 @@ #include "MachEnv.h" #include "OmegaKokkos.h" +#include + namespace OMEGA { // check if two real numbers are equal with a given relative tolerance @@ -86,10 +88,12 @@ template KOKKOS_FUNCTION int getVertBound(const T &VertBound, int 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, const VertMin &VMin, const VertMax &VMax, + VertArr ZCoord, // either array or nullptr to indicate 2D ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, SetBoundary SetBndOpt = SetBoundary::No) { @@ -98,8 +102,15 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, int NElementsOwned; int NElementsSize; Array1DReal XElement, YElement; + Array2DReal ZElement; Array1DReal LonElement, LatElement; + constexpr bool ZCoordNull = std::is_same_v; + + if constexpr (!ZCoordNull) { + ZElement = ZCoord; + } + switch (Element) { case OnCell: NElementsOwned = Mesh->NCellsOwned; @@ -139,41 +150,54 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, 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); - ScalarElement(IElement) = Fun(X, Y); + ScalarElement(IElement) = Fun(IElement, X, Y); } else { const Real Lon = LonElement(IElement); const Real Lat = LatElement(IElement); - ScalarElement(IElement) = Fun(Lon, Lat); + ScalarElement(IElement) = Fun(IElement, Lon, Lat); } }); } if constexpr (Array::rank == 2) { - const int NVertLayers = ScalarElement.extent_int(1); 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; + const int KMin = getVertBound(VMin, IElement); + const int KMax = getVertBound(VMax, IElement); parallelForInner( - Team, KRange, INNER_LAMBDA(int KOff) { - const int K = KMin + KOff; + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + Real X, Y; if (Geom == Geometry::Planar) { - const Real X = XElement(IElement); - const Real Y = YElement(IElement); - ScalarElement(IElement, K) = Fun(X, Y); + X = XElement(IElement); + Y = YElement(IElement); } else { - const Real Lon = LonElement(IElement); - const Real Lat = LatElement(IElement); - ScalarElement(IElement, K) = Fun(Lon, Lat); + X = LonElement(IElement); + Y = LatElement(IElement); } + + // Workaround for a CUDA issue with capturing variables + // inside `if constexpr` + const auto &LocZElement = ZElement; + const auto &LocFun = Fun; + constexpr auto LocZCoordNull = ZCoordNull; + + Real ScalarValue; + if constexpr (LocZCoordNull) { + ScalarValue = LocFun(IElement, X, Y); + } else { + const Real Z = LocZElement(IElement, K); + ScalarValue = LocFun(IElement, K, X, Y, Z); + } + + ScalarElement(IElement, K) = ScalarValue; }); }); } @@ -186,21 +210,34 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, 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; + const int KMin = getVertBound(VMin, IElement); + const int KMax = getVertBound(VMax, IElement); parallelForInner( - Team, KRange, INNER_LAMBDA(int KOff) { - const int K = KMin + KOff; + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + Real X, Y; if (Geom == Geometry::Planar) { - const Real X = XElement(IElement); - const Real Y = YElement(IElement); - ScalarElement(L, IElement, K) = Fun(X, Y); + X = XElement(IElement); + Y = YElement(IElement); } else { - const Real Lon = LonElement(IElement); - const Real Lat = LatElement(IElement); - ScalarElement(L, IElement, K) = Fun(Lon, Lat); + X = LonElement(IElement); + Y = LatElement(IElement); } + + // Workaround for a CUDA issue with capturing variables + // inside `if constexpr` + const auto &LocZElement = ZElement; + const auto &LocFun = Fun; + constexpr auto LocZCoordNull = ZCoordNull; + + Real ScalarValue; + if constexpr (LocZCoordNull) { + ScalarValue = LocFun(IElement, X, Y); + } else { + const Real Z = LocZElement(IElement, K); + ScalarValue = LocFun(IElement, K, X, Y, Z); + } + + ScalarElement(L, IElement, K) = ScalarValue; }); }); } @@ -222,17 +259,18 @@ int setScalar(const Functor &Fun, const Array &ScalarElement, Geometry Geom, const int VMin = 0; const int VMax = ScalarElement.extent_int(Array::rank - 1) - 1; return setScalar(Fun, ScalarElement, Geom, Mesh, Element, VMin, VMax, - ExchangeHalosOpt); + nullptr, 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, + const VertMin &VMin, const VertMax &VMax, VertArr ZCoord, ExchangeHalos ExchangeHalosOpt = ExchangeHalos::Yes, CartProjection CartProjectionOpt = CartProjection::Yes, SetBoundary SetBndOpt = SetBoundary::No) { @@ -260,14 +298,34 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, const int NEdgesSize = Mesh->NEdgesSize; const int NEdgesSet = Mesh->NEdgesOwned + static_cast(SetBndOpt); - auto ProjectVector = KOKKOS_LAMBDA(int IEdge) { + constexpr bool ZCoordNull = std::is_same_v; + + Array2DReal ZElement; + if constexpr (!ZCoordNull) { + ZElement = ZCoord; + } + + auto ProjectVector = KOKKOS_LAMBDA(int IEdge, int K) { Real VecFieldEdge; + + // Workaround for a CUDA issue with capturing variables inside `if + // constexpr` + const auto &LocZElement = ZElement; + const auto &LocFun = Fun; + constexpr auto LocZCoordNull = ZCoordNull; + if (Geom == Geometry::Planar) { const Real XE = XEdge(IEdge); const Real YE = YEdge(IEdge); Real VecField[2]; - Fun(VecField, XE, YE); + + if constexpr (LocZCoordNull) { + LocFun(VecField, IEdge, XE, YE); + } else { + const Real ZE = LocZElement(IEdge, K); + LocFun(VecField, IEdge, K, XE, YE, ZE); + } if (EdgeComp == EdgeComponent::Normal) { const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); @@ -287,7 +345,12 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, const Real LatE = LatEdge(IEdge); Real VecField[2]; - Fun(VecField, LonE, LatE); + if constexpr (LocZCoordNull) { + LocFun(VecField, IEdge, LonE, LatE); + } else { + const Real ZE = LocZElement(IEdge, K); + LocFun(VecField, IEdge, K, LonE, LatE, ZE); + } if (CartProjectionOpt == CartProjection::Yes) { Real VecFieldCart[3]; @@ -344,7 +407,7 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, if (SetBndOpt == SetBoundary::Yes) { IEdge = IEdge == 0 ? (NEdgesSize - 1) : (IEdge - 1); } - VectorFieldEdge(IEdge) = ProjectVector(IEdge); + VectorFieldEdge(IEdge) = ProjectVector(IEdge, 0); }); } @@ -354,13 +417,11 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, 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; + const int KMin = getVertBound(VMin, IEdge); + const int KMax = getVertBound(VMax, IEdge); parallelForInner( - Team, KRange, INNER_LAMBDA(int KOff) { - const int K = KMin + KOff; - VectorFieldEdge(IEdge, K) = ProjectVector(IEdge); + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + VectorFieldEdge(IEdge, K) = ProjectVector(IEdge, K); }); }); } @@ -385,7 +446,7 @@ int setVectorEdge(const Functor &Fun, const Array &VectorFieldEdge, 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); + nullptr, ExchangeHalosOpt, CartProjectionOpt); } template Real reduceArray(const Array1DReal &Arr, int Extent0) { diff --git a/components/omega/test/ocn/SubmesoEddiesTest.cpp b/components/omega/test/ocn/SubmesoEddiesTest.cpp new file mode 100644 index 000000000000..342263c68cf5 --- /dev/null +++ b/components/omega/test/ocn/SubmesoEddiesTest.cpp @@ -0,0 +1,750 @@ +#include "SubmesoEddies.h" +#include "FillValues.h" +#include "GlobalConstants.h" +#include "IO.h" +#include "IOStream.h" +#include "OceanDriver.h" +#include "OceanTestCommon.h" +#include "Pacer.h" +#include "TimeStepper.h" + +#include "mpi.h" +#include +#include + +using namespace OMEGA; + +constexpr Geometry Geom = Geometry::Spherical; + +KOKKOS_FUNCTION Real verticalSpacing(Real Lon, Real Lat, Real Scale) { + return Scale * 10 * + (0.5 + 0.3 * Kokkos::sin(Lon) * Kokkos::sin(Lon) * Kokkos::cos(Lat)); +} + +KOKKOS_FUNCTION Real surfaceHeight(Real Lon, Real Lat) { + return -20 * Kokkos::sin(Lon) * Kokkos::cos(Lat); +} + +KOKKOS_FUNCTION Real geomHeightMid(int K, Real Lon, Real Lat, Real Scale) { + const Real Z0 = surfaceHeight(Lon, Lat); + const Real DZ = verticalSpacing(Lon, Lat, Scale); + return Z0 - (K + 1) * DZ + DZ / 2; +} + +KOKKOS_FUNCTION Real geomHeightInterface(int K, Real Lon, Real Lat, + Real Scale) { + const Real Z0 = surfaceHeight(Lon, Lat); + const Real DZ = verticalSpacing(Lon, Lat, Scale); + return Z0 - K * DZ; +} + +constexpr Real H0 = 2e4; +constexpr Real Alpha = 0.4; + +KOKKOS_FUNCTION Real specVol(Real Lon, Real Lat, Real Z) { + return 1 / RhoSw * + (1 + Alpha * Kokkos::cos(Lon) * Kokkos::pow(Kokkos::cos(Lat), 4)) * + Kokkos::exp(-(Z * Z) / (H0 * H0)); +} + +KOKKOS_FUNCTION Real gradBuoyancyX(Real Lon, Real Lat, Real Z) { + return Alpha * Gravity / REarth * Kokkos::sin(Lon) * + Kokkos::pow(Kokkos::cos(Lat), 3) * Kokkos::exp(-(Z * Z) / (H0 * H0)); +} + +KOKKOS_FUNCTION Real gradBuoyancyY(Real Lon, Real Lat, Real Z) { + return Alpha * Gravity / REarth * 4 * Kokkos::cos(Lon) * + Kokkos::pow(Kokkos::cos(Lat), 3) * Kokkos::sin(Lat) * + Kokkos::exp(-(Z * Z) / (H0 * H0)); +} + +// Buoyancy gradient averaged between Z1 and Z2 +KOKKOS_FUNCTION Real meanGradBuoyancyX(Real Lon, Real Lat, Real Z1, Real Z2) { + const Real Tmp = Alpha * Gravity / REarth * Kokkos::sin(Lon) * + Kokkos::pow(Kokkos::cos(Lat), 3); + return Tmp * H0 * Kokkos::sqrt(Pi) / 2 * + (Kokkos::erf(Z2 / H0) - Kokkos::erf(Z1 / H0)) / (Z2 - Z1); +} + +KOKKOS_FUNCTION Real meanGradBuoyancyY(Real Lon, Real Lat, Real Z1, Real Z2) { + const Real Tmp = Alpha * Gravity / REarth * 4 * Kokkos::cos(Lon) * + Kokkos::pow(Kokkos::cos(Lat), 3) * Kokkos::sin(Lat); + return Tmp * H0 * Kokkos::sqrt(Pi) / 2 * + (Kokkos::erf(Z2 / H0) - Kokkos::erf(Z1 / H0)) / (Z2 - Z1); +} + +KOKKOS_FUNCTION Real bruntVaisalaFreqSq(Real Lon, Real Lat, Real Z) { + return Gravity * (-2 * Z / (H0 * H0)); +} + +// Brunt-Vaisala frequency averaged between Z1 and Z2 +KOKKOS_FUNCTION Real meanBVML(Real Lon, Real Lat, Real Z1, Real Z2) { + const Real Tmp = Gravity * (2 / (H0 * H0)); + return 2._Real / 3 * Kokkos::sqrt(Tmp) * + (Z2 * Kokkos::sqrt(-Z2) - Z1 * Kokkos::sqrt(-Z1)) / (Z2 - Z1); +} + +KOKKOS_FUNCTION Real shapeFunctionDeriv(Real Z, Real H) { + const Real Tmp = (2 * Z) / H + 1; + const Real Fac1 = (1 - Tmp * Tmp); + const Real Fac2 = (1 + 5 * Tmp * Tmp / 21); + + if (Fac1 * Fac2 >= 0) { + return -4 / H * Tmp * Fac2 + 20 / (21 * H) * Tmp * Fac1; + } else { + return 0; + } +} + +Error initSubmesoEddiesTest(MPI_Comm Comm, const std::string MeshFile, + int NVertLayers) { + Error Err; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + // Initialize the Logging system + initLogging(DefEnv); + + LOG_INFO("------ SubmesoEddies Unit Tests ------"); + + // Open config file + Config("Omega"); + Config::readAll("omega.yml"); + + // First step of time stepper initialization needed for IOstream + TimeStepper::init1(); + + // Get model clock + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); + + // Initialize the IO system + IO::init(DefComm); + + // Create the default decomposition (initializes the decomposition) + Decomp::init(MeshFile); + + // Initialize streams + IOStream::init(); + + // Initialize the default halo + Halo::init(); + + // Initialize the default mesh + HorzMesh::init(ModelClock); + + // Initialize the default vertical coordinate + if (NVertLayers == 0) { + VertCoord::init(); + } else { + VertCoord::init(false, NVertLayers); + } + + return Err; +} + +void finalizeSubmesoEddiesTest() { + + IOStream::finalize(); + Tracers::clear(); + VertAdv::clear(); + VertCoord::clear(); + TimeStepper::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); +} + +std::pair setupVerticalCoord() { + auto *Mesh = HorzMesh::getDefault(); + auto *VCoord = VertCoord::getDefault(); + + const auto &LonCellH = Mesh->LonCellH; + const auto &LatCellH = Mesh->LatCellH; + auto LonCell = createDeviceMirrorCopy(LonCellH); + auto LatCell = createDeviceMirrorCopy(LatCellH); + + const int NVertLayers = VCoord->NVertLayers; + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + auto &GeomZMid = VCoord->GeomZMid; + auto &GeomZInterface = VCoord->GeomZInterface; + + const Real Scale = 64.0 / NVertLayers; + + // Set vertical coordinates + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const Real Lon = LonCell(ICell); + const Real Lat = LatCell(ICell); + + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + GeomZMid(ICell, K) = geomHeightMid(K, Lon, Lat, Scale); + }); + + parallelForInner( + Team, Range{KMin, KMax + 1}, INNER_LAMBDA(int K) { + GeomZInterface(ICell, K) = + geomHeightInterface(K, Lon, Lat, Scale); + }); + }); + + // for this test we also need vertical coordinates on edges + Array2DReal GeomZMidEdge("GeomZMidEdge", Mesh->NEdgesSize, + VCoord->NVertLayers); + Array2DReal GeomZInterfaceEdge("GeomZInterfaceEdge", Mesh->NEdgesSize, + VCoord->NVertLayersP1); + + const auto &LonEdgeH = Mesh->LonEdgeH; + const auto &LatEdgeH = Mesh->LatEdgeH; + auto LonEdge = createDeviceMirrorCopy(LonEdgeH); + auto LatEdge = createDeviceMirrorCopy(LatEdgeH); + + const auto &MinLayerEdgeTop = VCoord->MinLayerEdgeTop; + const auto &MaxLayerEdgeBot = VCoord->MaxLayerEdgeBot; + + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const Real Lon = LonEdge(IEdge); + const Real Lat = LatEdge(IEdge); + + const int KMin = MinLayerEdgeTop(IEdge); + const int KMax = MaxLayerEdgeBot(IEdge); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + GeomZMidEdge(IEdge, K) = geomHeightMid(K, Lon, Lat, Scale); + }); + + parallelForInner( + Team, Range{KMin, KMax + 1}, INNER_LAMBDA(int K) { + GeomZInterfaceEdge(IEdge, K) = + geomHeightInterface(K, Lon, Lat, Scale); + }); + }); + + return {GeomZMidEdge, GeomZInterfaceEdge}; +} + +Array2DReal computePseudoThickOnEdges() { + auto *Mesh = HorzMesh::getDefault(); + const auto &CellsOnEdge = Mesh->CellsOnEdge; + + auto *VCoord = VertCoord::getDefault(); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + const auto &GeomZInterface = VCoord->GeomZInterface; + + Array2DReal PseudoThick("PseudoThick", Mesh->NCellsSize, + VCoord->NVertLayers); + + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + // The unit tests assume that pseudo thickness is equal to + // geometric thickness + PseudoThick(ICell, K) = + GeomZInterface(ICell, K) - GeomZInterface(ICell, K + 1); + }); + }); + + Array2DReal MeanPseudoThickEdge("MeanPseudoThickEdge", Mesh->NEdgesSize, + VCoord->NVertLayers); + + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + MeanPseudoThickEdge(IEdge, K) = + (PseudoThick(JCell0, K) + PseudoThick(JCell1, K)) / 2; + }); + }); + + return MeanPseudoThickEdge; +} + +Array2DReal computeExactGradInterface(const Array2DReal &GeomZInterfaceEdge) { + auto *Mesh = HorzMesh::getDefault(); + auto *VCoord = VertCoord::getDefault(); + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + + Array1DI4 MaxLayerEdgeTopP1("MaxLayerEdgeTopP1", Mesh->NEdgesSize); + + parallelFor( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { + MaxLayerEdgeTopP1(IEdge) = MaxLayerEdgeTop(IEdge) + 1; + }); + + Array2DReal ExactGradInterface("ExactGradInterface", Mesh->NEdgesSize, + VCoord->NVertLayersP1); + deepCopy(ExactGradInterface, FillValueReal); + + setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, int K, Real Lon, Real Lat, + Real Z) { + VecField[0] = gradBuoyancyX(Lon, Lat, Z); + VecField[1] = gradBuoyancyY(Lon, Lat, Z); + }, + ExactGradInterface, EdgeComponent::Normal, Geom, Mesh, MinLayerEdgeBot, + MaxLayerEdgeTopP1, GeomZInterfaceEdge, ExchangeHalos::No, + CartProjection::No); + + // Boundary edges should have fill value + parallelFor( + {Mesh->NEdgesOwned}, KOKKOS_LAMBDA(int IEdge) { + int MinLyrEdgeBot = MinLayerEdgeBot(IEdge); + int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge); + if (MaxLyrEdgeTop < MinLyrEdgeBot) { + ExactGradInterface(IEdge, MinLyrEdgeBot) = FillValueReal; + } + }); + + return ExactGradInterface; +} + +Array1DReal computeExactGradML(const Array2DReal &GeomZInterfaceEdge, + const Array1DReal &DenMixLayerDepth) { + auto *Mesh = HorzMesh::getDefault(); + const auto &CellsOnEdge = Mesh->CellsOnEdge; + + auto *VCoord = VertCoord::getDefault(); + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + + Array1DReal ExactGradML("ExactGradML", Mesh->NEdgesSize); + + setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real Lon, Real Lat) { + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + const Real MLDepthEdge = + Kokkos::min(DenMixLayerDepth(JCell0), DenMixLayerDepth(JCell1)); + + const Real Z0 = GeomZInterfaceEdge(IEdge, MinLayerEdgeBot(IEdge)); + + VecField[0] = meanGradBuoyancyX(Lon, Lat, Z0, Z0 - MLDepthEdge); + VecField[1] = meanGradBuoyancyY(Lon, Lat, Z0, Z0 - MLDepthEdge); + }, + ExactGradML, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No, + CartProjection::No); + + return ExactGradML; +} + +Array1DReal computeExactBVML(const Array2DReal &GeomZInterfaceEdge, + const Array1DReal &DenMixLayerDepth) { + auto *Mesh = HorzMesh::getDefault(); + const auto &CellsOnEdge = Mesh->CellsOnEdge; + + auto *VCoord = VertCoord::getDefault(); + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + + Array1DReal ExactBVML("ExactBVML", Mesh->NEdgesSize); + + setScalar( + KOKKOS_LAMBDA(int IEdge, Real Lon, Real Lat) { + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + const Real MLDepthEdge = + Kokkos::min(DenMixLayerDepth(JCell0), DenMixLayerDepth(JCell1)); + + const Real Z0 = GeomZInterfaceEdge(IEdge, MinLayerEdgeBot(IEdge)); + + return meanBVML(Lon, Lat, Z0, Z0 - MLDepthEdge); + }, + ExactBVML, Geom, Mesh, OnEdge, ExchangeHalos::No); + + return ExactBVML; +} + +// If SetupExact is true then this test is set up such that ML depth is located +// exactly on layer midpoints. Since no interpolation is necessary, +// very low error norms are expected. +Error testDenMixedLayerDepth(bool SetupExact) { + Error Err; + + auto *Mesh = HorzMesh::getDefault(); + + auto *VCoord = VertCoord::getDefault(); + const int NVertLayers = VCoord->NVertLayers; + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + auto &GeomZMid = VCoord->GeomZMid; + auto &GeomZInterface = VCoord->GeomZInterface; + + auto *SubEddies = SubmesoEddies::getInstance(); + + Array2DReal ReferenceSpecVol("ReferenceSpecVol", Mesh->NCellsAll, + VCoord->NVertLayers); + Array1DReal ExactMixLayerDepth("ExactMixLayerDepth", Mesh->NCellsAll); + + const Real DenThreshold = SubEddies->DenThreshold; + + parallelForOuter( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + + const Real MaxDepth = + GeomZInterface(ICell, KMin) - GeomZMid(ICell, KMax); + + const int KDen = Kokkos::max(10, ICell % NVertLayers); + + const bool KDenValid = KDen <= KMax; + + const Real DepthKDen = + KDenValid ? GeomZInterface(ICell, KMin) - GeomZMid(ICell, KDen) + : MaxDepth; + + const Real ExactDepth = SetupExact ? DepthKDen : 0.95 * DepthKDen; + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + const Real Depth = + GeomZInterface(ICell, KMin) - GeomZMid(ICell, K); + + Real Dens; + if (SetupExact) { + Dens = K < KDen ? RhoSw : RhoSw + (1 + 1e-6) * DenThreshold; + } else { + Dens = RhoSw + DenThreshold * (Depth / ExactDepth) * + (Depth / ExactDepth); + } + + ReferenceSpecVol(ICell, K) = 1._Real / Dens; + }); + + ExactMixLayerDepth(ICell) = KDenValid ? ExactDepth : MaxDepth; + }); + + SubEddies->computeDenMixLayerDepth(ReferenceSpecVol); + const auto &NumMixLayerDepth = SubEddies->DenMixLayerDepth; + + ErrorMeasures MixLayerDepthErrors; + computeErrors(MixLayerDepthErrors, NumMixLayerDepth, ExactMixLayerDepth, + Mesh, OnCell); + + if (SetupExact && + (MixLayerDepthErrors.L2 > 1e-7 || MixLayerDepthErrors.LInf > 1e-7)) { + Err += Error(ErrorCode::Fail, "denMixedLayerDepth Exact FAIL"); + } + + if (!SetupExact && + (MixLayerDepthErrors.L2 > 5e-2 || MixLayerDepthErrors.LInf > 5e-2)) { + Err += Error(ErrorCode::Fail, "denMixedLayerDepth NonExact FAIL"); + } + + return Err; +} + +Error testBuoyancyGrad(const Array2DReal &MeanPseudoThickEdge, + const Array2DReal &GeomZInterfaceEdge) { + Error Err; + + auto *Mesh = HorzMesh::getDefault(); + + auto *VCoord = VertCoord::getDefault(); + auto &GeomZMid = VCoord->GeomZMid; + auto &GeomZInterface = VCoord->GeomZInterface; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + // Compute specific volume + Array2DReal SpecVol("SpecVol", Mesh->NCellsSize, VCoord->NVertLayers); + setScalar( + KOKKOS_LAMBDA(int ICell, int K, Real Lon, Real Lat, Real Z) { + return specVol(Lon, Lat, Z); + }, + SpecVol, Geom, Mesh, OnCell, VCoord->MinLayerCell, VCoord->MaxLayerCell, + GeomZMid); + + // Compute squared Brunt-Vaisala frequency + Array1DI4 MaxLayerCellP1("MaxLayerCellP1", Mesh->NCellsSize); + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + MaxLayerCellP1(ICell) = MaxLayerCell(ICell) + 1; + }); + + Array2DReal BVFreqSq("BVFreqSq", Mesh->NCellsSize, VCoord->NVertLayersP1); + setScalar( + KOKKOS_LAMBDA(int ICell, int K, Real Lon, Real Lat, Real Z) { + return bruntVaisalaFreqSq(Lon, Lat, Z); + }, + BVFreqSq, Geom, Mesh, OnCell, VCoord->MinLayerCell, MaxLayerCellP1, + GeomZInterface); + + // Compute numerical buoyancy gradient + auto *SubEddies = SubmesoEddies::getInstance(); + + SubEddies->computeBuoyGrad(SpecVol, MeanPseudoThickEdge, GeomZMid, BVFreqSq); + + // Compute exact buoyancy gradient + Array2DReal ExactGradInterface = + computeExactGradInterface(GeomZInterfaceEdge); + + // Compute errors and check they are reasonable + ErrorMeasures InterfaceGradErrors; + computeErrors(InterfaceGradErrors, SubEddies->GradBuoyEdgeInterface, + ExactGradInterface, Mesh, OnEdge); + + const Real MaxL2Error = 4.1e-4; + if (InterfaceGradErrors.L2 > MaxL2Error) { + Err += Error(ErrorCode::Fail, "buoyancyGrad L2 FAIL, {:e} > {:e}", + InterfaceGradErrors.L2, MaxL2Error); + } + + const Real MaxLInfError = 5.5e-4; + if (InterfaceGradErrors.LInf > MaxLInfError) { + Err += Error(ErrorCode::Fail, "buoyancyGrad LInf FAIL, {:e} > {:e}", + InterfaceGradErrors.LInf, MaxLInfError); + } + + return Err; +} + +Error testEddyVelocity(const Array2DReal &GeomZInterfaceEdge, + const Array2DReal &MeanPseudoThickEdge) { + Error Err; + + auto *Mesh = HorzMesh::getDefault(); + const auto &CellsOnEdge = Mesh->CellsOnEdge; + const auto &CellsOnCell = Mesh->CellsOnCell; + const auto &NEdgesOnCell = Mesh->NEdgesOnCell; + const auto &DcEdge = Mesh->DcEdge; + const int NCellsAll = Mesh->NCellsAll; + + auto *VCoord = VertCoord::getDefault(); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + auto &GeomZMid = VCoord->GeomZMid; + auto &GeomZInterface = VCoord->GeomZInterface; + + auto *SubEddies = SubmesoEddies::getInstance(); + + const Real LfMin = 3500; + SubEddies->LfMin = LfMin; + const Real DsMax = SubEddies->DsMax; + + const auto &DenMixLayerDepth = SubEddies->DenMixLayerDepth; + const auto &DenMixLayerIndex = SubEddies->DenMixLayerIndex; + const auto &TimeScale = SubEddies->TimeScale; + const Real Ce = SubEddies->Ce; + + // Pick initial mixed layer indices arbitrarily + Array1DI4 DenMixLayerIndexTmp("DenMixLayerIndexTmp", Mesh->NCellsSize); + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + const int KMax = MaxLayerCell(ICell); + + int MLIndex = Kokkos::min(KMax, 15 + ICell % 10); + DenMixLayerIndexTmp(ICell) = MLIndex; + }); + + // Make sure there are no large differences in mixed layer depth between + // neighbouring cells + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + const int KMin = MinLayerCell(ICell); + + int MLIndex = DenMixLayerIndexTmp(ICell); + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + const int JCell = CellsOnCell(ICell, J); + if (JCell != NCellsAll) { + MLIndex = Kokkos::min(MLIndex, DenMixLayerIndexTmp(JCell)); + } + } + DenMixLayerIndex(ICell) = MLIndex; + DenMixLayerDepth(ICell) = + GeomZInterface(ICell, KMin) - GeomZMid(ICell, MLIndex); + }); + + // Compute exact buoyancy gradient + SubEddies->GradBuoyEdgeInterface = + computeExactGradInterface(GeomZInterfaceEdge); + + Array2DReal BVFreqSq("BVFreqSq", Mesh->NCellsSize, VCoord->NVertLayersP1); + + // Compute exact Brunt-Vaisala frequency + Array1DI4 MaxLayerCellP1("MaxLayerCellP1", Mesh->NCellsSize); + parallelFor( + {Mesh->NCellsAll}, KOKKOS_LAMBDA(int ICell) { + MaxLayerCellP1(ICell) = MaxLayerCell(ICell) + 1; + }); + + setScalar( + KOKKOS_LAMBDA(int ICell, int K, Real Lon, Real Lat, Real Z) { + return bruntVaisalaFreqSq(Lon, Lat, Z); + }, + BVFreqSq, Geom, Mesh, OnCell, VCoord->MinLayerCell, MaxLayerCellP1, + GeomZInterface); + + // Compute numerical eddy velocity + SubEddies->computeEddyVelocity(BVFreqSq, MeanPseudoThickEdge); + + const auto &EddyVelocity = SubEddies->EddyVelocity; + + // Compute exact mixed layer average of buoyancy gradient + Array1DReal MeanBuoyGrad = + computeExactGradML(GeomZInterfaceEdge, DenMixLayerDepth); + + // Compute exact mixed layer average of Brunt-Vaisala frequency + Array1DReal MeanBV = computeExactBVML(GeomZInterfaceEdge, DenMixLayerDepth); + + // Compute exact eddy velocity + Array2DReal ExactEddyVelocity("ExactEddyVelocity", Mesh->NEdgesSize, + VCoord->NVertLayers); + deepCopy(ExactEddyVelocity, FillValueReal); + + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int MinLyrEdgeBot = MinLayerEdgeBot(IEdge); + const int MaxLyrEdgeTop = MaxLayerEdgeTop(IEdge); + + const int JCell0 = CellsOnEdge(IEdge, 0); + const int JCell1 = CellsOnEdge(IEdge, 1); + + const Real TScale = TimeScale(IEdge); + + const Real MLDepthEdge = + Kokkos::min(DenMixLayerDepth(JCell0), DenMixLayerDepth(JCell1)); + + const Real Ds = Kokkos::min(DcEdge(IEdge), DsMax); + + const Real GradBuoyML = MeanBuoyGrad(IEdge); + const Real BVFreqML = MeanBV(IEdge); + + const Real Lf1 = + Kokkos::abs(GradBuoyML) * MLDepthEdge / (TScale * TScale); + const Real Lf2 = BVFreqML * MLDepthEdge / TScale; + + const Real Lf = Kokkos::max(LfMin, Kokkos::max(Lf1, Lf2)); + + const Real Factor = + Ce * Ds / Lf * MLDepthEdge * MLDepthEdge * GradBuoyML / TScale; + + if (MaxLyrEdgeTop >= MinLyrEdgeBot) { + + parallelForInner( + Team, Range{MinLyrEdgeBot, MaxLyrEdgeTop}, + INNER_LAMBDA(int K) { + const Real ZEdge = + 0.5_Real * + (GeomZMid(JCell0, K) - + GeomZInterface(JCell0, MinLayerCell(JCell0)) + + GeomZMid(JCell1, K) - + GeomZInterface(JCell1, MinLayerCell(JCell1))); + + ExactEddyVelocity(IEdge, K) = + -Factor * shapeFunctionDeriv(ZEdge, MLDepthEdge); + }); + } + }); + + // Compute errors and check that they are reasonable + ErrorMeasures EddyVelocityErrors; + computeErrors(EddyVelocityErrors, EddyVelocity, ExactEddyVelocity, Mesh, + OnEdge); + + const Real MaxL2Error = 2.4e-1; + if (EddyVelocityErrors.L2 > MaxL2Error) { + Err += Error(ErrorCode::Fail, "eddyVelocity L2 FAIL, {:e} > {:e}", + EddyVelocityErrors.L2, MaxL2Error); + } + + const Real MaxLInfError = 5.3e-1; + if (EddyVelocityErrors.LInf > MaxLInfError) { + Err += Error(ErrorCode::Fail, "eddyVelocity LInf FAIL, {:e} > {:e}", + EddyVelocityErrors.LInf, MaxLInfError); + } + + return Err; +} + +Error testSubmesoEddies(MPI_Comm Comm, std::string MeshFile, int NVertLayers) { + Error Err; + + // Initialize Omega modules needed for this test + Err += initSubmesoEddiesTest(Comm, MeshFile, NVertLayers); + + // Setup vertical coordinates + auto [GeomZMidEdge, GeomZInterfaceEdge] = setupVerticalCoord(); + + // Initialize submesoscale eddy parametrization + SubmesoEddies::init(); + + // Test retrieval + if (!SubmesoEddies::getInstance()) { + ABORT_ERROR("SubmesoEddiesTest: SubmesoEddies retrieval FAIL"); + } + + // Test mixed layer depth computation + for (bool SetupExact : {true, false}) { + Err += testDenMixedLayerDepth(SetupExact); + } + + // Compute pseudo-thickness on edges for the subsequent tests + auto MeanPseudoThickEdge = computePseudoThickOnEdges(); + + // Test buoyancy gradient computation + Err += testBuoyancyGrad(MeanPseudoThickEdge, GeomZInterfaceEdge); + + // Test eddy velocity computation + Err += testEddyVelocity(GeomZInterfaceEdge, MeanPseudoThickEdge); + + // Destroy submesoscale eddy parametrization + SubmesoEddies::destroyInstance(); + + finalizeSubmesoEddiesTest(); + + return Err; +} + +int main(int argc, char *argv[]) { + Error Err; + + const MPI_Comm Comm = MPI_COMM_WORLD; + + MPI_Init(&argc, &argv); + Pacer::initialize(Comm); + Pacer::setPrefix("Omega:"); + + try { + Kokkos::initialize(argc, argv); + { Err += testSubmesoEddies(Comm, "OmegaMesh.nc", 0); } + Kokkos::finalize(); + } catch (const std::exception &Ex) { + Err += Error(ErrorCode::Fail, Ex.what() + std::string(": FAIL")); + } catch (...) { + Err += Error(ErrorCode::Fail, "Unknown: FAIL"); + } + + CHECK_ERROR_ABORT(Err, "Submeso Eddies Unit Tests FAIL"); + + Pacer::finalize(); + MPI_Finalize(); + + return 0; + +} // end of main +//===-----------------------------------------------------------------------===/ diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index a747f87f6d63..07f43193bfdb 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -80,23 +80,25 @@ int initState() { deepCopy(TracersCell, NAN); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThickness(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.pseudoThickness(X, Y); + }, PseudoThickCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, - VCoord->MaxLayerCell, ExchangeHalos::Yes, SetBoundary::Yes); + VCoord->MaxLayerCell, nullptr, ExchangeHalos::Yes, SetBoundary::Yes); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.tracer(X, Y); }, TracersCell, Geom, Mesh, OnCell, VCoord->MinLayerCell, - VCoord->MaxLayerCell, ExchangeHalos::Yes, SetBoundary::Yes); + VCoord->MaxLayerCell, nullptr, ExchangeHalos::Yes, SetBoundary::Yes); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real Lon, Real Lat) { VecField[0] = Setup.velocityX(Lon, Lat); VecField[1] = Setup.velocityY(Lon, Lat); }, NormalVelEdge, EdgeComponent::Normal, Geom, Mesh, - VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, ExchangeHalos::Yes, - CartProjection::No, SetBoundary::Yes); + VCoord->MinLayerEdgeTop, VCoord->MaxLayerEdgeBot, nullptr, + ExchangeHalos::Yes, CartProjection::No, SetBoundary::Yes); return Err; } diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index e7f5038be0ce..a22c83311b06 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -347,7 +347,7 @@ int setupBottomDragTestFields(const int NVertLayers, const Real Coeff, // Note: this computes bottom drag at every layer. Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.bottomDragX(X, Y, Coeff); VecField[1] = Setup.bottomDragY(X, Y, Coeff); }, @@ -365,20 +365,20 @@ int setupBottomDragTestFields(const int NVertLayers, const Real Coeff, }); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, NormalVelEdge, EdgeComponent::Normal, Geom, Mesh); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarA(X, Y) * Setup.scalarA(X, Y) / 2; }, KECell, Geom, Mesh, OnCell); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarB(X, Y); }, PseudoThickEdge, Geom, Mesh, OnEdge); return Err; @@ -397,7 +397,9 @@ int testThickFluxDiv(int NVertLayers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.divergence(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return -Setup.divergence(X, Y); + }, ExactThickFluxDiv, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input array @@ -408,7 +410,7 @@ int testThickFluxDiv(int NVertLayers, Real RTol) { deepCopy(OnesEdge, 1); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, @@ -453,7 +455,7 @@ int testPotVortHAdv(int NVertLayers, Real RTol) { NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = (Setup.normRelVort(X, Y) + Setup.normPlanetVort(X, Y)) * Setup.pseudoThick(X, Y) * Setup.vectorX(X, Y); VecField[1] = (Setup.normRelVort(X, Y) + Setup.normPlanetVort(X, Y)) * @@ -467,26 +469,32 @@ int testPotVortHAdv(int NVertLayers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normRelVort(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { + return Setup.normRelVort(X, Y); + }, NormRelVortEdge, Geom, Mesh, OnEdge); Array2DReal NormPlanetVortEdge("NormPlanetVortEdge", Mesh->NEdgesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normPlanetVort(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { + return Setup.normPlanetVort(X, Y); + }, NormPlanetVortEdge, Geom, Mesh, OnEdge); Array2DReal PseudoThickEdge("PseudoThickEdge", Mesh->NEdgesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.pseudoThick(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { + return Setup.pseudoThick(X, Y); + }, PseudoThickEdge, Geom, Mesh, OnEdge); Array2DReal NormVelEdge("NormVelEdge", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, @@ -530,7 +538,7 @@ int testKEGrad(int NVertLayers, Real RTol) { Array2DReal ExactKEGrad("ExactKEGrad", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = -Setup.gradX(X, Y); VecField[1] = -Setup.gradY(X, Y); }, @@ -540,8 +548,8 @@ int testKEGrad(int NVertLayers, Real RTol) { Array2DReal KECell("KECell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, KECell, - Geom, Mesh, OnCell); + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalar(X, Y); }, + KECell, Geom, Mesh, OnCell); // Compute numerical result Array2DReal NumKEGrad("NumKEGrad", Mesh->NEdgesOwned, NVertLayers); @@ -579,7 +587,7 @@ int testSSHGrad(int NVertLayers, Real RTol) { Array2DReal ExactSSHGrad("ExactSSHGrad", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = -Gravity * Setup.gradX(X, Y); VecField[1] = -Gravity * Setup.gradY(X, Y); }, @@ -588,8 +596,8 @@ int testSSHGrad(int NVertLayers, Real RTol) { Array1DReal SSHCell("SSHCell", Mesh->NCellsSize); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, SSHCell, - Geom, Mesh, OnCell); + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalar(X, Y); }, + SSHCell, Geom, Mesh, OnCell); // Compute numerical result Array2DReal NumSSHGrad("NumSSHGrad", Mesh->NEdgesOwned, NVertLayers); @@ -635,7 +643,7 @@ int testVelDiff(int NVertLayers, Real RTol) { Array2DReal ExactVelDiff("ExactVelDiff", Mesh->NEdgesOwned, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = ViscDel2 * Setup.laplaceVecX(X, Y); VecField[1] = ViscDel2 * Setup.laplaceVecY(X, Y); }, @@ -645,14 +653,16 @@ int testVelDiff(int NVertLayers, Real RTol) { Array2DReal DivCell("DivCell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.divergence(X, Y); + }, DivCell, Geom, Mesh, OnCell); Array2DReal RVortVertex("RVortVertex", Mesh->NVerticesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.curl(X, Y); }, RVortVertex, - Geom, Mesh, OnVertex); + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { return Setup.curl(X, Y); }, + RVortVertex, Geom, Mesh, OnVertex); // Compute numerical result Array2DReal NumVelDiff("NumVelDiff", Mesh->NEdgesOwned, NVertLayers); @@ -705,7 +715,7 @@ int testVelHyperDiff(int NVertLayers, Real RTol) { NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = -ViscDel4 * Setup.laplaceVecX(X, Y); VecField[1] = -ViscDel4 * Setup.laplaceVecY(X, Y); }, @@ -715,14 +725,16 @@ int testVelHyperDiff(int NVertLayers, Real RTol) { Array2DReal DivCell("DivCell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.divergence(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.divergence(X, Y); + }, DivCell, Geom, Mesh, OnCell); Array2DReal RVortVertex("RVortVertex", Mesh->NVerticesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.curl(X, Y); }, RVortVertex, - Geom, Mesh, OnVertex); + KOKKOS_LAMBDA(int IVertex, Real X, Real Y) { return Setup.curl(X, Y); }, + RVortVertex, Geom, Mesh, OnVertex); // Compute numerical result Array2DReal NumVelHyperDiff("NumVelHyperDiff", Mesh->NEdgesOwned, @@ -763,7 +775,7 @@ int testSfcStressForcing(int NVertLayers) { // Note: this computes surface stress forcing at every layer Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.sfcStressForcingX(X, Y); VecField[1] = Setup.sfcStressForcingY(X, Y); }, @@ -779,7 +791,7 @@ int testSfcStressForcing(int NVertLayers) { Array1DReal NormalStressEdge("NormalStressEdge", Mesh->NEdgesSize); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, @@ -789,7 +801,7 @@ int testSfcStressForcing(int NVertLayers) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return Setup.scalarB(X, Y); }, PseudoThickEdge, Geom, Mesh, OnEdge); // Compute numerical result @@ -910,7 +922,7 @@ int testImplicitBottomDrag(int NVertLayers, Real RTol) { deepCopy(NumImplicitBottomDrag, 0.0_Real); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarB(X, Y); }, PseudoThickCell, Geom, Mesh, OnCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -986,7 +998,7 @@ int testBottomDragInactiveEdges(int NVertLayers) { Array2DReal NormalVelEdge("NormalVelEdge", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, @@ -995,7 +1007,7 @@ int testBottomDragInactiveEdges(int NVertLayers) { Array2DReal KECell("KECell", Mesh->NCellsSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarA(X, Y) * Setup.scalarA(X, Y) / 2; }, KECell, Geom, Mesh, OnCell); @@ -1004,7 +1016,7 @@ int testBottomDragInactiveEdges(int NVertLayers) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return Setup.scalarB(X, Y); }, PseudoThickEdge, Geom, Mesh, OnEdge); // Save the layer indices so they can be restored for the later tests @@ -1084,14 +1096,16 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracerFluxDiv(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.tracerFluxDiv(X, Y); + }, ExactTrFluxDiv, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays Array2DReal NormalVelocity("NormalVelocity", Mesh->NEdgesSize, NVertLayers); Err += setVectorEdge( - KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + KOKKOS_LAMBDA(Real(&VecField)[2], int IEdge, Real X, Real Y) { VecField[0] = Setup.vectorX(X, Y); VecField[1] = Setup.vectorY(X, Y); }, @@ -1101,12 +1115,14 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { Array2DReal ThickEdge("ThickEdh", Mesh->NEdgesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.pseudoThick(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return -Setup.pseudoThick(X, Y); + }, TrCell, Geom, Mesh, OnCell); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return 1; }, ThickEdge, Geom, Mesh, - OnEdge); + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return 1; }, ThickEdge, Geom, + Mesh, OnEdge); // Compute numerical result Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, @@ -1154,7 +1170,9 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracerDiff(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.tracerDiff(X, Y); + }, ExactTracerDiff, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays @@ -1162,14 +1180,14 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarA(X, Y); }, TracerCell, Geom, Mesh, OnCell); Array2DReal PseudoThickEdge("PseudoThickEdge", Mesh->NEdgesSize, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int IEdge, Real X, Real Y) { return Setup.scalarB(X, Y); }, PseudoThickEdge, Geom, Mesh, OnEdge); // Compute numerical result @@ -1212,7 +1230,9 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { Mesh->NCellsOwned, NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.tracerHyperDiff(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return -Setup.tracerHyperDiff(X, Y); + }, ExactTracerHyperDiff, Geom, Mesh, OnCell, ExchangeHalos::No); // Set input arrays @@ -1220,7 +1240,7 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { NVertLayers); Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarC(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarC(X, Y); }, TrDel2Cell, Geom, Mesh, OnCell); // Compute numerical result @@ -1296,7 +1316,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { // Set Input Field values. Use a combination of scalarB and vectorX to // ensure the full surface tracer restoring logic is exercised. Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { return Setup.scalarB(X, Y) + Setup.vectorX(X, Y); }, InputField, Geom, Mesh, OnCell); @@ -1304,7 +1324,9 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { // Set TracersOnCell values (use scalarB for simplicity, but could be // any field). Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + KOKKOS_LAMBDA(int ICell, Real X, Real Y) { + return Setup.scalarB(X, Y); + }, TracersOnCell, Geom, Mesh, OnCell); parallelFor( {NTracers, Mesh->NCellsSize, NVertLayers},