diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 08164bdaa373..cc7e62cfae65 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -32,6 +32,9 @@ Omega: FluxThicknessType: Center HorzTracerFluxOrder: 2 VerticalTracerFluxLimiterEnable: true + HorzTracerFluxLimiterEnable: false + HorzTracerFluxLimiterBudgetsEnable: false + HorzTracerFluxLimiterMonotonicityCheckEnable: false VerticalTracerFluxOrder: 3 SfcStress: InterpType: Isotropic diff --git a/components/omega/doc/design/Tendencies.md b/components/omega/doc/design/Tendencies.md index cfc5e00ca364..89d2b6bf86c1 100644 --- a/components/omega/doc/design/Tendencies.md +++ b/components/omega/doc/design/Tendencies.md @@ -52,11 +52,22 @@ class Tendencies{ VelocityHyperDiffOnEdge VelocityHyperDiff; SfcStressForcingOnEdge SfcStressForcing; BottomDragOnEdge BottomDrag; + TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; - TracerHorzAdvOnCell TracerHorzAdv; - TracerHighOrderHorzAdvOnCell TracerHighOrderHorzAdv; + SurfaceTracerRestoringOnCell SurfaceTracerRestoring; + private: + const HorzMesh *Mesh; ///< Pointer to horizontal mesh + VertCoord *VCoord; ///< Pointer to vertical coordinate + VertAdv *VAdv; ///< Pointer to vertical advection + CustomTendencyType CustomThicknessTend; + CustomTendencyType CustomVelocityTend; + Eos *EqState; ///< Pointer to equation of state + PressureGrad *PGrad; ///< Pointer to pressure gradient + VertMix *VMix; ///< Pointer to vertical mixing + I4 NTracers; ///< Number of tracers + TimeInterval TimeStep; ///< Time step static Tendencies *DefaultTendencies; static std::map> AllTendencies; }; @@ -106,7 +117,7 @@ void Tendencies::computeAllTendencies(const OceanState *State, const AuxilarySta ``` The layer thickness tendencies will be computed with a method: ```c++ -void Tendencies::computeThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time); +void Tendencies::computePseudoThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time); ``` The normal velocity tendencies will be computed with a method: ```c++ diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 5fa72197132f..aa99ac309ea4 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -38,7 +38,6 @@ implemented: - `SfcStressForcingOnEdge` - `BottomDragOnEdge` - `TracerHorzAdvOnCell` -- `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` - `SurfaceTracerRestoringOnCell` diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index e6ee85a5045c..f6a20ec5c3a7 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -15,7 +15,6 @@ tendency terms are currently implemented: | VelocityDiffusionOnEdge | Laplacian horizontal mixing, defined on edges | VelocityHyperDiffOnEdge | biharmonic horizontal mixing, defined on edges | TracerHorzAdvOnCell | horizontal advection of thickness-weighted tracers -| TracerHighOrderHorzAdvOnCell | second order horizontal advection of thickness-weighted tracers | TracerDiffOnCell | horizontal diffusion of thickness-weighted tracers | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers | SfcStressForcingOnEdge | forcing by surface stress (e.g. wind), defined on edges @@ -45,9 +44,11 @@ the currently available tendency terms: | | ViscDel4 | horizontal biharmonic mixing coefficient for normal velocity | | DivFactor | scale factor for the divergence term | TracerHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term -| | HorzTracerFluxOrder | 1 for standard linear advection -| TracerHighOrderHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term -| | HorzTracerFluxOrder | 2 for second order advection algorithm +| | HorzTracerFluxOrder | 2 for standard linear advection +| | HorzTracerFluxOrder | 3 for second order advection algorithm +| | HorzTracerFluxLimiterEnable | enable/disable monotonic flux corrected transport (FCT) +| | HorzTracerFluxLimiterBudgetsEnable | enable/disable budgets if FCT is enabled +| | HorzTracerFluxLimiterMonotonicityCheckEnable | enable/disable check for mon-mononic values if FCT is enabled | TracerDiffOnCell | TracerDiffTendencyEnable | enable/disable term | | EddyDiff2 | horizontal diffusion coefficient | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term @@ -92,13 +93,13 @@ $$ $$ Where $u$ is $\mathbf{V}\cdot\mathbf{n}$ where $\mathbf{n}$ is the unit normal along the edge being evaluated and $F$ is the evaluation of $\psi$ for the elements on each side of the edge being evaluated. -This is used when the TracerHorzAdvOnCell user option is active and HorzTracerFluxOrder is 1. +This is used when the TracerHorzAdvOnCell user option is active and HorzTracerFluxOrder is 2. To make the comparison to higher order methods explicit, this first order method is equivalent to defining $\psi$ as a linear function, $\psi = c_0 + c_x x + c_y y$, between the two elements sharing the edge $i$ and taking the directed derivative along the edge. For higher order flux calculations, higher order derivatives of $\psi$ are needed and for the user option of -TracerHorzAdvOnCell along with HorzTracerFluxOrder set to 2, a quadratic approximation of $\psi$ is used, +TracerHorzAdvOnCell along with HorzTracerFluxOrder set to 3, a quadratic approximation of $\psi$ is used, $$ \psi = c_0 + c_x x + c_y y + c_{xx} x^2 + c_{xy} xy + c_{yy} y^2, $$ @@ -138,6 +139,60 @@ is about order 1.7 as shown in {numref}`tracer-higher-order-convergence`: Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate ``` +### Monotonic Flux Limiting Transport for Second Order Horizontal Advection + +It is well known that higher order (order 2 and above) scheme for numerically integrating fluxes suffer +from dispersive "ripples" in the results particularly near steep gradients. Lower order schemes +produce no ripples but suffer from excessive numerical diffusion. Flux-corrected transport (FCT) is a +technique which embodies the best of both schemes. + +The implementation of a monotonic flux-corrected transport (FCT) in Omega follows a standard algorithm as given in + +Zalesak, S. T. (1979). Fully multidimensional flux-corrected transport algorithms for fluids. +Journal of Computational Physics, 31(3), 335–362. DOI: 10.1016/0021-9991(79)90051-2 + +The procedure is as follows given the notation above: +1. Compute $F^L_{i+1/2}(u\psi)$, the transportive flux by some low order scheme guarenteed to give +monotonic (ripple-free) results. +2. Compute $F^H_{i+1/2}(u\psi)$, the transportative flux by some high order scheme. +3. Define the "antififfusive flux": +$$ +A_{i+1/2}(u\psi) = F^H_{i+1/2}(u\psi) - F^L_{i+1/2}(u\psi) +$$ +4. Compute the updated low order ("transported and diffused") "td" solution: +$$ + w^{td}_i = w^n_i - \frac{1}{\Delta x}[F_{i+1/2}(u\psi) - F_{i-1/2}(u\psi)] +$$ +5. Limit he $A_{i+1/2}(u\psi)$ in a manner such that $w^{n+1}$ as computed below is free of extrema +not found in $w^{td}$ = $w^n$; +$$ +A^C_{i+1/2}(u\psi) = C_{i+1/2}(u\psi) A_{i+1/2}(u\psi), \quad \quad 0 \leq C_{i+1/2} \leq 1 +$$ +6. Finally apply hte limited antidiffusive fluxes: +$$ + w^{n+1}_i = w^{td}_i - \frac{1}{\Delta x}[A^C_{i+1/2}(u\psi) - A^C_{i-1/2}(u\psi)] +$$ + + +```{figure} images/higher_order_tracer_convergence_on_sphere_FCT.jpeg +:name: tracer-higher-order-convergence_FCT +:align: center +:width: 600 px +Tracer higer order convergence example of a cosine bell advected on a sphere with FCT showing an order 2.23 convergence rate +``` + + +When monotonic flux limiting is used then there are two diagnostic options that can be enabled, +HorzTracerFluxLimiterBudgetsEnable and HorzTracerFluxLimiterMonotonicityCheckEnable. The +HorzTracerFluxLimiterMonotonicityCheckEnable flag enables a post-transport check that the +resulting values are indeed monotone as required by the algorithm and any discrepancies are +printed along with information on where they occur. + +The HorzTracerFluxLimiterBudgetsEnable option enables the output of two diagnostic fields to the +output mesh file. One is the "FCTActiveTracerHorizontalAdvectionEdgeFlux" variable which is the +edge flux across every side of every element. The other is the "FCTActiveTracerHorizontalAdvectionTendency" +variable which is the cell-centered value of the advection tendency. + ## See Also Additional information on forcing (currently wind forcing and surface tracer diff --git a/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere_FCT.jpeg b/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere_FCT.jpeg new file mode 100644 index 000000000000..4f5557164610 Binary files /dev/null and b/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere_FCT.jpeg differ diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index a85e24d37dc6..c5a20c9222de 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -7,8 +7,8 @@ // tendency terms are enabled. // //===----------------------------------------------------------------------===// +#include -#include "Tendencies.h" #include "CustomTendencyTerms.h" #include "Eos.h" #include "Error.h" @@ -17,6 +17,7 @@ #include "OceanState.h" #include "PGrad.h" #include "Pacer.h" +#include "Tendencies.h" #include "TimeStepper.h" #include "Tracers.h" #include "VertAdv.h" @@ -226,6 +227,22 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config if (Order == 4) { this->TracerHorzAdv.Coef3rdOrder = 0; } + Err += AdvectConfig.get("HorzTracerFluxLimiterEnable", TracerHorzAdv.FCT); + CHECK_ERROR_ABORT( + Err, + "Tendencies: HorzTracerFluxLimiterEnable not found in AdvectConfig"); + if (TracerHorzAdv.FCT) { + Err += AdvectConfig.get("HorzTracerFluxLimiterBudgetsEnable", + TracerHorzAdv.ComputeBudgets); + CHECK_ERROR_ABORT(Err, + "Tendencies: HorzTracerFluxLimiterBudgetsEnable not " + "found in AdvectConfig"); + Err += AdvectConfig.get("HorzTracerFluxLimiterMonotonicityCheckEnable", + TracerHorzAdv.MonotonicityCheck); + CHECK_ERROR_ABORT( + Err, "Tendencies: HorzTracerFluxLimiterMonotonicityCheckEnable " + "not found in AdvectConfig"); + } } Err += TendConfig.get("TracerDiffTendencyEnable", this->TracerDiffusion.Enabled); @@ -419,7 +436,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), - TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), + TracerHorzAdv(Mesh, VCoord, VAdv), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad), VMix(VMix) { @@ -559,7 +576,7 @@ void Tendencies::computeVelocityTendenciesOnly( AuxState->PseudoThicknessAux.FluxPseudoThickEdge; const Array2DReal &NormRVortEdge = AuxState->VorticityAux.NormRelVortEdge; const Array2DReal &NormFEdge = AuxState->VorticityAux.NormPlanetVortEdge; - Array2DReal NormVelEdge = State->getNormalVelocity(VelTimeLevel); + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); if (LocPotentialVortHAdv.Enabled) { Pacer::start("Tend:PotentialVortHAdv", 2); parallelForOuter( @@ -572,7 +589,7 @@ void Tendencies::computeVelocityTendenciesOnly( Team, KRange, INNER_LAMBDA(int KChunk) { LocPotentialVortHAdv(LocNormalVelocityTend, IEdge, KChunk, NormRVortEdge, NormFEdge, - FluxPseudoThickEdge, NormVelEdge); + FluxPseudoThickEdge, NormalVelEdge); }); }); Pacer::stop("Tend:PotentialVortHAdv", 2); @@ -653,7 +670,7 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::start("Tend:computeVelocityVAdvTend", 2); // Compute velocity tendency from vertical advection - VAdv->computeVelocityVAdvTend(NormalVelocityTend, NormVelEdge, + VAdv->computeVelocityVAdvTend(NormalVelocityTend, NormalVelEdge, FluxPseudoThickEdge); Pacer::stop("Tend:computeVelocityVAdvTend", 2); @@ -684,7 +701,7 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::start("Tend:bottomDrag", 2); parallelFor( {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { - LocBottomDrag(LocNormalVelocityTend, IEdge, NormVelEdge, KECell, + LocBottomDrag(LocNormalVelocityTend, IEdge, NormalVelEdge, KECell, MeanPseudoThickEdge); }); Pacer::stop("Tend:bottomDrag", 2); @@ -729,52 +746,211 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); - OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocMinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(LocMaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); + Array2DReal PseudoThickCell = State->getPseudoThickness(ThickTimeLevel); + Array2DReal FluxPseudoThickEdge = + AuxState->PseudoThicknessAux.FluxPseudoThickEdge; + + OMEGA_SCOPE(LocPseudoThickCell, PseudoThickCell); + OMEGA_SCOPE(LocNormalVelEdge, NormalVelEdge); + OMEGA_SCOPE(LocFluxPseudoThickEdge, FluxPseudoThickEdge); Pacer::start("Tend:computeTracerTendenciesOnly", 1); parallelForOuter( {NTracers, Mesh->NCellsAll}, KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); parallelForInner( Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { LocTracerTend(L, ICell, K) = 0; }); }); - // compute tracer horizotal advection - Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); - const Array2DReal &FluxPseudoThickEdge = - AuxState->PseudoThicknessAux.FluxPseudoThickEdge; + R8 Dt = 0; + TimeStep.get(Dt, TimeUnits::Seconds); if (LocTracerHorzAdv.Enabled) { Pacer::start("Tend:tracerHorzAdv", 2); - parallelForOuter( - {NTracers, Mesh->NEdgesAll}, - KOKKOS_LAMBDA(int L, 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) { - LocTracerHorzAdv(L, IEdge, KChunk, TracerArray, - FluxPseudoThickEdge, NormalVelEdge); - }); - }); - parallelForOuter( - {NTracers, Mesh->NCellsAll}, - KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocTracerHorzAdv(LocTracerTend, L, ICell, KChunk); - }); - }); + if (LocTracerHorzAdv.FCT) { + parallelForOuter( + {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTProvisionaLayerThicknesses( + ICell, KChunk, Dt, LocFluxPseudoThickEdge, + LocPseudoThickCell, LocNormalVelEdge); + }); + }); + Kokkos::fence(); + for (int L = 0; L < NTracers; ++L) { + parallelForOuter( + {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTTracerCurFill(L, ICell, KChunk, + TracerArray); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTTracerMinMax( + ICell, KChunk, LocMinLayerCell, LocMaxLayerCell); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTHighAndLowOrderFlux( + IEdge, KChunk, LocMinLayerCell, LocMaxLayerCell, + LocFluxPseudoThickEdge, LocNormalVelEdge); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NCellsAll + 1}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTInitFluxInOut(ICell, KChunk); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NCellsHaloH(0)}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTFluxInOut(ICell, KChunk, Dt, + LocPseudoThickCell); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NEdgesHaloH(0)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTRescaleHighOrderFlux(IEdge, + KChunk); + }); + }); + Kokkos::fence(); + parallelForOuter( + {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTAccumulateHighOrderFlux( + L, ICell, KChunk, Dt, LocTracerTend, + LocPseudoThickCell); + }); + }); + Kokkos::fence(); + if (LocTracerHorzAdv.ComputeBudgets) { + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTComputeBudgetAdvectionEdgeFlux( + L, IEdge, KChunk); + }); + }); + parallelForOuter( + {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTComputeBudgetAdvectionTendency( + L, ICell, KChunk); + }); + }); + } + if (LocTracerHorzAdv.MonotonicityCheck) { + parallelForOuter( + {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv.FCTMonotonicityCheck( + L, ICell, KChunk, TracerArray); + }); + }); + } + } + } else { + parallelForOuter( + {NTracers, Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int L, int IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv(L, IEdge, KChunk, TracerArray, + LocFluxPseudoThickEdge, + LocNormalVelEdge); + }); + }); + parallelForOuter( + {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocTracerHorzAdv(LocTracerTend, L, ICell, KChunk); + }); + }); + } Pacer::stop("Tend:tracerHorzAdv", 2); } @@ -786,8 +962,8 @@ void Tendencies::computeTracerTendenciesOnly( parallelForOuter( {NTracers, Mesh->NCellsAll}, KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( @@ -806,8 +982,8 @@ void Tendencies::computeTracerTendenciesOnly( parallelForOuter( {NTracers, Mesh->NCellsAll}, KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( @@ -842,7 +1018,7 @@ void Tendencies::computeTracerTendenciesOnly( parallelFor( {NTracersToRestore, Mesh->NCellsAll}, KOKKOS_LAMBDA(int R, int ICell) { - const int KMin = MinLayerCell(ICell); + const int KMin = LocMinLayerCell(ICell); const int L = TracerIdsToRestore(R); LocSurfaceTracerRestoring(LocTracerTend, L, ICell, KMin, TracersMonthlySurfClimo, TracerArray); @@ -861,11 +1037,11 @@ void Tendencies::computePseudoThicknessTendencies( TimeInstant Time ///< [in] Time ) { // only need PseudoThicknessAux on edge - Array2DReal PseudoThick = State->getPseudoThickness(ThickTimeLevel); - Array2DReal NormVel = State->getNormalVelocity(VelTimeLevel); + Array2DReal PseudoThick = State->getPseudoThickness(ThickTimeLevel); + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); OMEGA_SCOPE(PseudoThicknessAux, AuxState->PseudoThicknessAux); OMEGA_SCOPE(PseudoThickCell, PseudoThick); - OMEGA_SCOPE(NormalVelEdge, NormVel); + OMEGA_SCOPE(LocNormalVelEdge, NormalVelEdge); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -882,7 +1058,7 @@ void Tendencies::computePseudoThicknessTendencies( parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { PseudoThicknessAux.computeVarsOnEdge( - IEdge, KChunk, PseudoThickCell, NormalVelEdge); + IEdge, KChunk, PseudoThickCell, LocNormalVelEdge); }); }); Pacer::stop("Tend:computePseudoThickAux", 2); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 5b142888142c..abefbe7be16f 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -7,13 +7,14 @@ // defined in the corresponding header file. // //===----------------------------------------------------------------------===// +#include +#include -#include "TendencyTerms.h" -#include "AuxiliaryState.h" #include "DataTypes.h" +#include "Error.h" #include "HorzMesh.h" #include "HorzOperators.h" -#include "OceanState.h" +#include "TendencyTerms.h" #include "Tracers.h" namespace OMEGA { @@ -72,23 +73,37 @@ BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh, MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, - const VertCoord *VCoord) + const VertCoord *VCoord, + const VertAdv *VAdv) : HorzontalMesh(Mesh), VerticalCoord(VCoord), + NVertLayers(VCoord->NVertLayers), NAdvCellsForEdge("NumberOfCellsContribToAdvectionAtEdge", Mesh->NEdgesAll), AdvCellsForEdge("IndexOfCellsContributingToAdvection", Mesh->NEdgesAll, Mesh->MaxEdges2 + 2), AdvMaskHighOrder("MaskForHighOrderAdvectionTerms", Mesh->NEdgesAll, VCoord->NVertLayers), + CellsOnCell(Mesh->CellsOnCell), AdvCoefs("CommonAdvectionCoefficients", Mesh->MaxEdges2 + 2, Mesh->NEdgesAll), AdvCoefs3rd("CommonAdvectionCoeffsForHighOrder", Mesh->MaxEdges2 + 2, Mesh->NEdgesAll), HighOrderFlxHorz("HigherOrderHorizontalFlux", Tracers::getNumTracers(), Mesh->NEdgesAll, VCoord->NVertLayers), - NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), - CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), - DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {} + TracerCur(), NEdgesOnCell(Mesh->NEdgesOnCell), + EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MinLayerEdgeTop), + EdgeSignOnCell(Mesh->EdgeSignOnCell), DvEdge(Mesh->DvEdge), + AreaCell(Mesh->AreaCell), + VerticalPseudoVelocity(VAdv->VerticalPseudoVelocity), HProvInv(), + HNewInv(), HProv(), TracerMax(), TracerMin(), HighOrderFlx(), + LowOrderFlx(), MinLayerCell(VCoord->MinLayerCell), + MaxLayerCell(VCoord->MaxLayerCell), WorkTend(), FlxIn(), FlxOut(), + ActiveTracerHorizontalAdvectionEdgeFlux(), + ActiveTracerHorizontalAdvectionTendency() { + deepCopy(HighOrderFlxHorz, 0); +} TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) @@ -142,6 +157,70 @@ void TracerHorzAdvOnCell::init() { parallelFor( {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); }); Kokkos::fence(); + if (FCT) { + HProvInv = Array2DReal("FCTProvesionalLayerThickness", Mesh->NEdgesAll, + NVertLayers); + HNewInv = Array2DReal("FCTProvesionalNewInverse", NEdgesAll, NVertLayers); + HProv = Array2DReal("FCTProvesionalThickness", NCellsAll, NVertLayers); + TracerCur = Array2DReal("TracerCur", NCellsAll + 1, NVertLayers), + TracerMax = Array2DReal("FCTTracerMax", NCellsAll, NVertLayers); + TracerMin = Array2DReal("FCTTracerMin", NCellsAll, NVertLayers); + HighOrderFlx = + Array2DReal("FCTHighOrderFlx", std::max(NEdgesAll, NCellsAll) + 1, + NVertLayers + 1); + LowOrderFlx = + Array2DReal("FCTLowOrderFlx", std::max(NEdgesAll, NCellsAll) + 1, + NVertLayers + 1); + WorkTend = Array2DReal("WorkTend", NCellsAll + 1, NVertLayers); + FlxIn = Array2DReal("FlxIn", NCellsAll + 1, NVertLayers); + FlxOut = Array2DReal("FlxOut", NCellsAll + 1, NVertLayers); + if (ComputeBudgets) { + const int NTracers = Tracers::getNumTracers(); + const int NEdges = Mesh->NEdgesHaloH(1); + ActiveTracerHorizontalAdvectionEdgeFlux = + Array3DReal("FCTActiveTracerHorizontalAdvectionEdgeFlux", NTracers, + NEdges, NVertLayers); + ActiveTracerHorizontalAdvectionTendency = + Array3DReal("FCTActiveTracerHorizontalAdvectionTendency", NTracers, + NCellsAll, NVertLayers); + deepCopy(ActiveTracerHorizontalAdvectionEdgeFlux, 0.0); + deepCopy(ActiveTracerHorizontalAdvectionTendency, 0.0); + const int NDims = 1; + const std::string GroupName = "AuxiliaryState"; + std::vector FluxDimNames(NDims, "NEdges"); + auto BudgetAdvectionEdgeFlux = Field::create( + ActiveTracerHorizontalAdvectionEdgeFlux.label(), // field name + "Tracer FCT Horizontal Advection Edge Flux Budget", // long name or + // description + "", // units + "", // CF standard Name + 0, // min valid value + std::numeric_limits::max(), // max valid value + NDims, // number of dimensions + FluxDimNames // dimension names + ); + BudgetAdvectionEdgeFlux->attachData( + ActiveTracerHorizontalAdvectionEdgeFlux); + FieldGroup::addFieldToGroup( + ActiveTracerHorizontalAdvectionEdgeFlux.label(), GroupName); + std::vector TendDimNames(NDims, "NCells"); + auto BudgetAdvectionCellTend = Field::create( + ActiveTracerHorizontalAdvectionTendency.label(), // field name + "Tracer FCT Horizontal Advection Cell Flux " + "Tendency", // long name or description + "", // units + "", // CF standard Name + 0, // min valid value + std::numeric_limits::max(), // max valid value + NDims, // number of dimensions + TendDimNames // dimension names + ); + BudgetAdvectionCellTend->attachData( + ActiveTracerHorizontalAdvectionTendency); + FieldGroup::addFieldToGroup( + ActiveTracerHorizontalAdvectionTendency.label(), GroupName); + } + } } } // end namespace OMEGA diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 389d6bb1aafb..756c360a2f11 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -12,12 +12,14 @@ #include "AuxiliaryState.h" #include "GlobalConstants.h" +#include "Halo.h" #include "HorzMesh.h" #include "MachEnv.h" #include "OceanState.h" #include "VertCoord.h" #include // for std::copysign +#include namespace OMEGA { @@ -163,7 +165,6 @@ class KEGradOnEdge { /// Gradient of sea surface height defined on edges multipled by gravitational /// acceleration, for momentum equation -/// NOTE: This term is only appropriate for shallow water (Omega v0) simulations class SSHGradOnEdge { public: bool Enabled = false; @@ -367,11 +368,16 @@ class BottomDragOnEdge { // Tracer horizontal advection term class TracerHorzAdvOnCell { public: - bool Enabled = false; - bool ForceLowOrder = false; + bool Enabled = true; + bool ForceLowOrder = false; + bool FCT = false; + bool ComputeBudgets = false; + bool MonotonicityCheck = false; // coefficient for blending high-order terms Real Coef3rdOrder = 0.25; - TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord, + const VertAdv *VAdv); + virtual ~TracerHorzAdvOnCell() {} void init(); KOKKOS_FUNCTION void operator()(const I4 L, const I4 IEdge, const I4 KChunk, const Array3DReal &TracerCell, @@ -424,29 +430,336 @@ class TracerHorzAdvOnCell { for (int I = 0; I < NEdgesOnCell(ICell); ++I) { const I4 IEdge = EdgesOnCell(ICell, I); - for (int K = KStart; K < KEnd; ++K) { + for (I4 K = KStart; K < KEnd; ++K) { Tend(L, ICell, K) += EdgeSignOnCell(ICell, I) * HighOrderFlxHorz(L, IEdge, K) * InvAreaCell; } } } - private: + KOKKOS_FUNCTION void + FCTProvisionaLayerThicknesses(const I4 ICell, const I4 KChunk, const Real Dt, + const Array2DReal &FluxPseudoThickEdge, + const Array2DReal &LayerThickness, + const Array2DReal &NormVelEdge) const { + + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + + const Real InvAreaCell = Dt / AreaCell(ICell); + for (I4 K = KStart; K < KEnd; ++K) { + HProv(ICell, K) = LayerThickness(ICell, K); + } + for (I4 I = 0; I < NEdgesOnCell(ICell); ++I) { + const I4 IEdge = EdgesOnCell(ICell, I); + const Real SignedFactor = + InvAreaCell * DvEdge(IEdge) * EdgeSignOnCell(ICell, I); + // Provisional layer thickness is after horizontal + // thickness flux only + for (I4 K = KStart; K < KEnd; ++K) { + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + HProv(ICell, K) += SignedFactor * NormalThicknessFlux; + } + } + // New layer thickness is after horizontal and vertical + // thickness flux + for (I4 K = KStart; K < KEnd; ++K) { + HProvInv(ICell, K) = 1.0_Real / HProv(ICell, K); + HNewInv(ICell, K) = + 1.0_Real / + (HProv(ICell, K) - Dt * VerticalPseudoVelocity(ICell, K) + + Dt * VerticalPseudoVelocity(ICell, K + 1)); + } + } + + KOKKOS_FUNCTION void FCTTracerCurFill(const I4 L, const I4 ICell, + const I4 KChunk, + const Array3DReal &TracerArray) const { + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) + TracerCur(ICell, K) = TracerArray(L, ICell, K); + } + KOKKOS_FUNCTION void FCTTracerMinMax(const I4 ICell, const I4 KChunk, + const Array1DI4 &MinLayerCell, + const Array1DI4 &MaxLayerCell) const { + + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) { + TracerMin(ICell, K) = TracerCur(ICell, K); + TracerMax(ICell, K) = TracerCur(ICell, K); + } + for (I4 I = 0; I < NEdgesOnCell(ICell); ++I) { + const I4 ICell2 = CellsOnCell(ICell, I); + const I4 KMin = MinLayerCell(ICell2); + const I4 KMax = MaxLayerCell(ICell2); + const I4 KRange = vertRangeChunked(KMin, KMax); + for (I4 KChunk = 0; KChunk < KRange; ++KChunk) { + const I4 KStart1 = Kokkos::max(KStart, KChunk * VecLength); + const I4 KEnd1 = Kokkos::min(KEnd, KStart1 + VecLength); + for (I4 K = KStart1; K < KEnd1; ++K) { + TracerMax(ICell, K) = + Kokkos::max(TracerMax(ICell, K), TracerCur(ICell2, K)); + TracerMin(ICell, K) = + Kokkos::min(TracerMin(ICell, K), TracerCur(ICell2, K)); + } + } + } + } + + KOKKOS_FUNCTION void FCTHighAndLowOrderFlux( + const I4 IEdge, const I4 KChunk, const Array1DI4 &MinLayerCell, + const Array1DI4 &MaxLayerCell, const Array2DReal &FluxPseudoThickEdge, + const Array2DReal &NormVelEdge) const { + const Real Coef3rdOrder = 0.25; + + const I4 ICell1 = CellsOnEdge(IEdge, 0); + const I4 ICell2 = CellsOnEdge(IEdge, 1); + + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + + // compute some common intermediate factors + for (I4 K = KStart; K < KEnd; ++K) + HighOrderFlx(IEdge, K) = 0.0_Real; + + // Compute 3rd or 4th fluxes where requested. + for (int I = 0; I < NAdvCellsForEdge(IEdge); ++I) { + const I4 ICell = AdvCellsForEdge(IEdge, I); + const Real Coef1 = AdvCoefs(I, IEdge); + const Real Coef3 = AdvCoefs3rd(I, IEdge) * Coef3rdOrder; + const I4 KStartCell = MinLayerCell(ICell); + const I4 KLenCell = MaxLayerCell(ICell); + const I4 KEndCell = KStartCell + KLenCell; + for (I4 K = KStartCell; K <= KEndCell; ++K) { + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + const Real TracerWgt = + NormalThicknessFlux * + (Coef1 + Coef3 * std::copysign(1.0_Real, NormalThicknessFlux)); + HighOrderFlx(IEdge, K) += + TracerWgt * TracerCur(ICell, K) * AdvMaskHighOrder(IEdge, K); + } + } + // Compute 2nd order fluxes where needed. + // Also compute low order upwind horizontal flux (monotonic) + // Remove low order flux from the high order flux + // Store left over high order flux in highOrderFlx array + for (I4 K = MinLayerEdgeBot(IEdge); K <= MaxLayerEdgeTop(IEdge); ++K) { + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + const Real TracerWeight = (1.0_Real - AdvMaskHighOrder(IEdge, K)) * + (DvEdge(IEdge) * 0.5_Real) * + NormalThicknessFlux; + + LowOrderFlx(IEdge, K) = + DvEdge(IEdge) * (Kokkos::max(0.0_Real, NormalThicknessFlux) * + TracerCur(ICell1, K) + + Kokkos::min(0.0_Real, NormalThicknessFlux) * + TracerCur(ICell2, K)); + + HighOrderFlx(IEdge, K) += + TracerWeight * (TracerCur(ICell1, K) + TracerCur(ICell2, K)); + HighOrderFlx(IEdge, K) -= LowOrderFlx(IEdge, K); + } + } + + KOKKOS_FUNCTION void FCTInitFluxInOut(const I4 ICell, + const I4 KChunk) const { + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) + WorkTend(ICell, K) = 0; + for (I4 K = KStart; K < KEnd; ++K) + FlxIn(ICell, K) = 0; + for (I4 K = KStart; K < KEnd; ++K) + FlxOut(ICell, K) = 0; + } + + KOKKOS_FUNCTION void FCTFluxInOut(const I4 ICell, const I4 KChunk, + const Real Dt, + const Array2DReal &LayerThickness) const { + const Real InvAreaCell = 1._Real / AreaCell(ICell); + const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); + const I4 KEndCell = KStartCell + KLenCell - 1; + for (I4 K = KStartCell; K <= KEndCell; ++K) { + // Finish computing the low order horizontal fluxes + // Upwind fluxes are accumulated in workTend + for (I4 I = 0; I < NEdgesOnCell(ICell); ++I) { + const I4 IEdge = EdgesOnCell(ICell, I); + const Real SignedFactor = EdgeSignOnCell(ICell, I) * InvAreaCell; + for (I4 K = MinLayerEdgeBot(IEdge); K <= MaxLayerEdgeTop(IEdge); + ++K) { + // Here workTend is the advection tendency due to the + // upwind (low order) fluxes. + WorkTend(ICell, K) += SignedFactor * LowOrderFlx(IEdge, K); + + // Accumulate remaining high order fluxes + FlxOut(ICell, K) += + Kokkos::min(0.0_Real, SignedFactor * HighOrderFlx(IEdge, K)); + FlxIn(ICell, K) += + Kokkos::max(0.0_Real, SignedFactor * HighOrderFlx(IEdge, K)); + } + } + // Build the factors for the FCT + // Computed using the bounds that were computed previously, + // and the bounds on the newly updated value + // Factors are placed in the flxIn and flxOut arrays + for (I4 K = MinLayerCell(ICell); K <= MaxLayerCell(ICell); ++K) { + // Here workTend is the upwind tendency + const Real TracerUpwindNew = + (TracerCur(ICell, K) * LayerThickness(ICell, K) + + Dt * WorkTend(ICell, K)) * + HProvInv(ICell, K); + const Real TracerMinNew = + TracerUpwindNew + Dt * FlxOut(ICell, K) * HProvInv(ICell, K); + const Real TracerMaxNew = + TracerUpwindNew + Dt * FlxIn(ICell, K) * HProvInv(ICell, K); + const Real ScaleFactorIn = (TracerMax(ICell, K) - TracerUpwindNew) / + (TracerMaxNew - TracerUpwindNew + Eps); + FlxIn(ICell, K) = + Kokkos::min(1.0_Real, Kokkos::max(0.0_Real, ScaleFactorIn)); + const Real ScaleFactorOut = + (TracerUpwindNew - TracerMin(ICell, K)) / + (TracerUpwindNew - TracerMinNew + Eps); + FlxOut(ICell, K) = + Kokkos::min(1.0_Real, Kokkos::max(0.0_Real, ScaleFactorOut)); + } + } + } + + KOKKOS_FUNCTION void FCTRescaleHighOrderFlux(const I4 IEdge, + const I4 KChunk) const { + const I4 ICell1 = CellsOnEdge(IEdge, 0); + const I4 ICell2 = CellsOnEdge(IEdge, 1); + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) { + HighOrderFlx(IEdge, K) = + Kokkos::max(0.0_Real, HighOrderFlx(IEdge, K)) * + Kokkos::min(FlxOut(ICell1, K), FlxIn(ICell2, K)) + + Kokkos::min(0.0_Real, HighOrderFlx(IEdge, K)) * + Kokkos::min(FlxIn(ICell1, K), FlxOut(ICell2, K)); + } + } + + KOKKOS_FUNCTION void + FCTAccumulateHighOrderFlux(const I4 L, const I4 ICell, const I4 KChunk, + const Real Dt, const Array3DReal &TracerArray, + const Array2DReal &LayerThickness) const { + + // Accumulate the scaled high order vertical tendencies + // and the upwind tendencies + const Real InvAreaCell1 = 1.0_Real / AreaCell(ICell); + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + + // Accumulate the scaled high order horizontal tendencies + for (I4 I = 0; I < NEdgesOnCell(ICell); ++I) { + const I4 IEdge = EdgesOnCell(ICell, I); + const Real SignedFactor = InvAreaCell1 * EdgeSignOnCell(ICell, I); + for (I4 K = MinLayerEdgeBot(IEdge); K <= MaxLayerEdgeTop(IEdge); ++K) { + // WorkTend on RHS is upwind tendency + // WorkTend on LHS is total horiz advect tendency + WorkTend(ICell, K) += SignedFactor * HighOrderFlx(IEdge, K); + } + } + for (I4 K = KStart; K < KEnd; ++K) { + // workTend on RHS is total horiz advection tendency + // TracerCur on LHS is provisional tracer after + // horizontal fluxes only. + TracerCur(ICell, K) = (TracerCur(ICell, K) * LayerThickness(ICell, K) + + Dt * WorkTend(ICell, K)) * + HProvInv(ICell, K); + TracerArray(L, ICell, K) += WorkTend(ICell, K); + } + } + + KOKKOS_FUNCTION void + FCTComputeBudgetAdvectionEdgeFlux(const I4 L, const I4 IEdge, + const I4 KChunk) const { + // Save u*h*T flux on edge for analysis. This variable will be + // divided by h at the end of the time step. + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) { + ActiveTracerHorizontalAdvectionEdgeFlux(L, IEdge, K) = + (LowOrderFlx(IEdge, K) + HighOrderFlx(IEdge, K)) / DvEdge(IEdge); + } + } + + KOKKOS_FUNCTION void + FCTComputeBudgetAdvectionTendency(const I4 L, const I4 ICell, + const I4 KChunk) const { + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) { + ActiveTracerHorizontalAdvectionTendency(L, ICell, K) = + WorkTend(ICell, K); + } + } + + KOKKOS_FUNCTION void + FCTMonotonicityCheck(const I4 L, const I4 ICell, const I4 KChunk, + const Array3DReal &TracerArray) const { + // Check tracer values against local min,max to detect + // non-monotone values and write warning if found + + const I4 KStart = KChunk * VecLength; + const I4 KEnd = KStart + VecLength; + for (I4 K = KStart; K < KEnd; ++K) { + if (TracerArray(L, ICell, K) < TracerMin(ICell, K) - Eps) { + printf("Horizontal minimum out of bounds on tracer: %i %lg %lg\n", + L, TracerMin(ICell, K), TracerArray(L, ICell, K)); + } + if (TracerArray(L, ICell, K) > TracerMax(ICell, K) + Eps) { + printf("Horizontal maximum out of bounds on tracer: %i %lg %lg\n", + L, TracerMax(ICell, K), TracerArray(L, ICell, K)); + } + } + } + + protected: + const Real Eps = 1.e-10_Real; const HorzMesh *HorzontalMesh; const VertCoord *VerticalCoord; + const I4 NVertLayers; + Array1DI4 NAdvCellsForEdge; Array2DI4 AdvCellsForEdge; Array2DI4 AdvMaskHighOrder; + Array2DI4 CellsOnCell; Array2DReal AdvCoefs; Array2DReal AdvCoefs3rd; Array3DReal HighOrderFlxHorz; + Array2DReal TracerCur; Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; Array2DI4 CellsOnEdge; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; Array2DReal EdgeSignOnCell; Array1DReal DvEdge; Array1DReal AreaCell; + Array2DReal VerticalPseudoVelocity; + Array2DReal HProvInv; + Array2DReal HNewInv; + Array2DReal HProv; + Array2DReal TracerMax; + Array2DReal TracerMin; + Array2DReal HighOrderFlx; + Array2DReal LowOrderFlx; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array2DReal WorkTend; + Array2DReal FlxIn; + Array2DReal FlxOut; + Array3DReal ActiveTracerHorizontalAdvectionEdgeFlux; + Array3DReal ActiveTracerHorizontalAdvectionTendency; }; // Tracer horizontal diffusion term diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index fb2b178d634f..b12d8bc43dce 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -56,6 +56,7 @@ struct TestSetupPlane { 0.00134324628763667899}; ErrorMeasures ExpectedTrHAdvErrors = {0.0029211089892916243, 0.0024583038518548855}; + ErrorMeasures ExpectedFCTTracerHorzAdvErrors = {0.00, 0.00}; ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, 0.00290978146207349032}; ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, @@ -196,6 +197,7 @@ struct TestSetupSphere { 0.270546491554748}; ErrorMeasures ExpectedTrHAdvErrors = {0.013259410329645643, 0.004094907022292395}; + ErrorMeasures ExpectedFCTTracerHorzAdvErrors = {0.00, 0.00}; ErrorMeasures ExpectedTrDel2Errors = {0.04865718541236144, 0.005105510870642706}; ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, @@ -848,7 +850,6 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -884,7 +885,11 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { // Compute numerical result Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, NVertLayers); - TracerHorzAdvOnCell TrHorzAdvOnC(Mesh, VCoord); + + VertAdv::init(); + const auto VAdv = VertAdv::getDefault(); + + TracerHorzAdvOnCell TrHorzAdvOnC(Mesh, VCoord, VAdv); TrHorzAdvOnC.ForceLowOrder = true; TrHorzAdvOnC.init(); @@ -894,11 +899,13 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { TrHorzAdvOnC(L, IEdge, KLayer, TrCell, ThickEdge, NormalVelocity); }); + Kokkos::fence(); parallelFor( {NTracers, Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int L, int ICell, int KLayer) { TrHorzAdvOnC(NumTrFluxDiv, L, ICell, KLayer); }); + Kokkos::fence(); ErrorMeasures TrHAdvErrors; Err += @@ -911,6 +918,132 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { LOG_INFO("TendencyTermsTest: TracerHorzAdv PASS"); } + VertAdv::clear(); + return Err; +} // end testTracerHorzAdvOnCell + +int testFCTTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { + + class TracerHorzAdvOnCellTest : public TracerHorzAdvOnCell { + public: + TracerHorzAdvOnCellTest(const HorzMesh *Mesh, const VertCoord *VCoord, + const VertAdv *VAdv) + : TracerHorzAdvOnCell(Mesh, VCoord, VAdv){}; + Array2DReal GetHProvInv() { return HProvInv; }; + Array2DReal GetHProv() { return HProv; }; + Array2DReal GetHNewInv() { return HNewInv; }; + Array2DReal GetTracerCur() { return TracerCur; }; + }; + + I4 Err = 0; + TestSetup Setup; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + // Compute exact result + Array3DReal ExactTrFluxDiv("ExactTrFluxDiv", NTracers, Mesh->NCellsOwned, + NVertLayers); + + Err += setScalar( + KOKKOS_LAMBDA(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) { + VecField[0] = Setup.vectorX(X, Y); + VecField[1] = Setup.vectorY(X, Y); + }, + NormalVelocity, EdgeComponent::Normal, Geom, Mesh); + + Array3DReal TrCell("TrCell", NTracers, Mesh->NCellsSize, NVertLayers); + Array2DReal FluxPseudoThickEdge("FluxPseudoThickEdge", Mesh->NEdgesSize, + NVertLayers); + Array2DReal LayerThickness("LayerThickness", Mesh->NCellsSize, NVertLayers); + Array2DReal NormVelEdge("NormVelEdge", Mesh->NEdgesSize, NVertLayers); + Array2DReal HProvInvExact("HProvInvExact", Mesh->NCellsSize, NVertLayers); + Array2DReal HProvExact("HProvExact", Mesh->NCellsSize, NVertLayers); + Array2DReal HNewInvExact("HProvExact", Mesh->NCellsSize, NVertLayers); + Array2DReal TracerSubView("TracerSubView", Mesh->NCellsSize, NVertLayers); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.pseudoThick(X, Y); }, + TrCell, Geom, Mesh, OnCell); + + deepCopy(FluxPseudoThickEdge, 1._Real); + deepCopy(LayerThickness, 1._Real); + deepCopy(HProvExact, 1._Real); + deepCopy(HProvInvExact, 1._Real); + deepCopy(HNewInvExact, 1._Real); + deepCopy(NormVelEdge, 1._Real); + + // Compute numerical result + Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, + NVertLayers); + + VertAdv::init(); + const auto VAdv = VertAdv::getDefault(); + deepCopy(VAdv->VerticalPseudoVelocity, 10._Real); + + TracerHorzAdvOnCellTest TrHorzAdvOnC(Mesh, VCoord, VAdv); + TrHorzAdvOnC.ForceLowOrder = false; + TrHorzAdvOnC.FCT = true; + TrHorzAdvOnC.init(); + + const Real Dt = 1; + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int KLayer) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, KLayer, Dt, FluxPseudoThickEdge, LayerThickness, + NormVelEdge); + }); + Kokkos::fence(); + + const auto HProvInv = TrHorzAdvOnC.GetHProvInv(); + const auto HProv = TrHorzAdvOnC.GetHProv(); + const auto HNewInv = TrHorzAdvOnC.GetHNewInv(); + // for (int i=0; iNCellsOwned; ++i) + // for (int j=0; jNCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(int ICell, int KLayer) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, KLayer, TrCell); + }); + const auto TracerCur = TrHorzAdvOnC.GetTracerCur(); + const auto Tracers = Kokkos::subview(TrCell, L, Kokkos::ALL, Kokkos::ALL); + deepCopy(TracerSubView, Tracers); + Err += computeErrors(FCTTracerHorzAdvErrors, TracerCur, TracerSubView, + Mesh, OnCell); + Err += + checkErrors("TendencyTermsTest", "TracerCur", FCTTracerHorzAdvErrors, + Setup.ExpectedFCTTracerHorzAdvErrors, RTol); + } + + if (Err == 0) { + LOG_INFO("TendencyTermsTest: FCTTracerHorzAdv PASS"); + } + + VertAdv::clear(); return Err; } // end testTracerHorzAdvOnCell @@ -1217,6 +1350,8 @@ int tendencyTermsTest(const std::string &MeshFile = DefaultMeshFile) { Err += testTracerHorzAdvOnCell(NVertLayers, NTracers, RTol); + Err += testFCTTracerHorzAdvOnCell(NVertLayers, NTracers, RTol); + Err += testTracerDiffOnCell(NVertLayers, NTracers, RTol); Err += testTracerHyperDiffOnCell(NVertLayers, NTracers, RTol);