diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9f589f9466cf..b19536e27664 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 fdca63beac38..6d5252a6fbfa 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 ExplicitBottomDrag; + 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 3db2b4c10098..4aedbd264bbd 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). Requires HorzTracerFluxOrder to be greater than 2. +| | 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 @@ -94,13 +95,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, $$ @@ -140,6 +141,87 @@ 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. + + +### Monotonic Flux Limiting Transport Example of Limiting Ripple Effect + +The effect of monotonic higher order FCT horizontal transport verses just higher order horizontal transport +can be seen in the Polaris test of the slotted cylinder. The cylinder is advected for one revolution under +a constant velocity field. The result is advecting the solution back to the original configuration as +an exact solution. The difference between the initial condition and the advected solution is the error +in the advection scheme. This difference is shown in the following two plots for non-monotonic and monotonic +advection. + +```{figure} images/higher_order_slotted_cylinder_non_monotonic_diff.jpg +:name: tracer-error-in-higher-order-convergence_non-monotonic +:align: center +:width: 600 px +Error in a tracer higher order convergence example of a slotted cylinder advected on a sphere with the standard third order advection +algorithm with no FCT flux correction. Notice the ringing due to the sharp edges of the tracer distribution. +``` + +```{figure} images/higher_order_slotted_cylinder_FCT_diff.jpg +:name: tracer-error-in-higher-order-convergence_FCT +:align: center +:width: 600 px +Error Tracer higher order convergence example of a slotted cylinder advected on a sphere with FCT monotone +flux correction. Notice the reduction in ringing around the sharp edges of the tracer distribution. +``` + + ## See Also Additional information on forcing (currently wind forcing and surface tracer diff --git a/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_FCT_diff.jpg b/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_FCT_diff.jpg new file mode 100644 index 000000000000..80da864e47bc Binary files /dev/null and b/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_FCT_diff.jpg differ diff --git a/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_non_monotonic_diff.jpg b/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_non_monotonic_diff.jpg new file mode 100644 index 000000000000..34e5553315f8 Binary files /dev/null and b/components/omega/doc/userGuide/images/higher_order_slotted_cylinder_non_monotonic_diff.jpg differ 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 5bb131f32ab6..f71753b35fee 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" @@ -240,6 +241,27 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config if (Order == 4) { this->TracerHorzAdv.Coef3rdOrder = 0; } + Err += AdvectConfig.get("HorzTracerFluxLimiterEnable", TracerHorzAdv.FCT); + OMEGA_REQUIRE( + !TracerHorzAdv.FCT || 2 < Order, + "HorzTracerFluxOrder: Since HorzTracerFluxLimiterEnable is true, " + "HorzTracerFluxOrder must be greater than 2. Found Order={}", + Order); + 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); @@ -476,7 +498,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), SfcStressForcing(Mesh, VCoord), ExplicitBottomDrag(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) { @@ -655,7 +677,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( @@ -668,7 +690,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); @@ -749,7 +771,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); @@ -780,7 +802,7 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::start("Tend:explicitBottomDrag", 2); parallelFor( {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { - LocExplicitBottomDrag(LocNormalVelocityTend, IEdge, NormVelEdge, + LocExplicitBottomDrag(LocNormalVelocityTend, IEdge, NormalVelEdge, KECell, MeanPseudoThickEdge); }); Pacer::stop("Tend:explicitBottomDrag", 2); @@ -820,57 +842,136 @@ void Tendencies::computeTracerTendenciesOnly( int VelTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { + OMEGA_SCOPE(LocTracerArray, TracerArray); OMEGA_SCOPE(LocTracerTend, TracerTend); OMEGA_SCOPE(LocTracerHorzAdv, TracerHorzAdv); 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}, + "Tend:LocTracerTend.init", {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( + "Tend:FCTProvisionaLayerThicknesses", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTProvisionaLayerThicknesses( + Team, ICell, Dt, LocFluxPseudoThickEdge, LocPseudoThickCell, + LocNormalVelEdge); + }); + for (int L = 0; L < NTracers; ++L) { + parallelForOuter( + "Tend:FCTTracerCurFill", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTTracerCurFill(Team, L, ICell, + LocTracerArray); + }); + parallelForOuter( + "Tend:FCTTracerMinMax", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTTracerMinMax(Team, ICell); + }); + parallelForOuter( + "Tend:FCTHighAndLowOrderFlux", {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + LocTracerHorzAdv.FCTHighAndLowOrderFlux( + Team, IEdge, LocFluxPseudoThickEdge, LocNormalVelEdge); + }); + parallelForOuter( + "Tend:FCTFluxInOut", {Mesh->NCellsHaloH(0)}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTFluxInOut(Team, ICell, Dt, + LocPseudoThickCell); + }); + parallelForOuter( + "Tend:FCTRescaleHighOrderFlux", {Mesh->NEdgesHaloH(0)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + LocTracerHorzAdv.FCTRescaleHighOrderFlux(Team, IEdge); + }); + parallelForOuter( + "Tend:FCTAccumulateHighOrderFlux", {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const Array2DReal Tend = Kokkos::subview( + LocTracerTend, L, Kokkos::ALL, Kokkos::ALL); + LocTracerHorzAdv.FCTAccumulateHighOrderFlux( + Team, ICell, Dt, Tend, LocPseudoThickCell); + }); + if (LocTracerHorzAdv.ComputeBudgets) { + parallelForOuter( + "Tend:FCTComputeBudgetAdvectionEdgeFlux", + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + LocTracerHorzAdv.FCTComputeBudgetAdvectionEdgeFlux( + Team, L, IEdge); + }); + parallelForOuter( + "Tend:FCTComputeBudgetAdvectionTendency", + {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTComputeBudgetAdvectionTendency( + Team, L, ICell); + }); + } + if (LocTracerHorzAdv.MonotonicityCheck) { + const I4 NCellsOwned = Mesh->NCellsOwned; + parallelForOuter( + "Tend:FCTMonotonicityCheck", {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + LocTracerHorzAdv.FCTMonotonicityCheck(Team, ICell); + }); + } + } + } else { + parallelForOuter( + "Tend:TracerHorzAdv.init", {NTracers, Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int L, int IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocTracerHorzAdv(L, IEdge, K, LocTracerArray, + LocFluxPseudoThickEdge, + LocNormalVelEdge); + }); + }); + parallelForOuter( + "Tend:TracerHorzAdv.exec", {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocTracerHorzAdv(LocTracerTend, L, ICell, K); + }); + }); + } Pacer::stop("Tend:tracerHorzAdv", 2); } @@ -880,16 +981,16 @@ void Tendencies::computeTracerTendenciesOnly( if (LocTracerDiffusion.Enabled) { Pacer::start("Tend:tracerDiffusion", 2); parallelForOuter( - {NTracers, Mesh->NCellsAll}, + "Tend::TracerDiffusion", {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( Team, KRange, INNER_LAMBDA(int KChunk) { LocTracerDiffusion(LocTracerTend, L, ICell, KChunk, - TracerArray, MeanPseudoThickEdge); + LocTracerArray, MeanPseudoThickEdge); }); }); Pacer::stop("Tend:tracerDiffusion", 2); @@ -902,8 +1003,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( @@ -923,7 +1024,7 @@ void Tendencies::computeTracerTendenciesOnly( } else if (VAdv->VertAdvChoice == VertAdvOption::FCT) { ThicknessForVAdv = AuxState->PseudoThicknessAux.ProvPseudoThickness; } - VAdv->computeTracerVAdvTend(LocTracerTend, TracerArray, ThicknessForVAdv, + VAdv->computeTracerVAdvTend(LocTracerTend, LocTracerArray, ThicknessForVAdv, TimeStep); Pacer::stop("Tend:computeTracerVAdvTend", 2); @@ -938,10 +1039,10 @@ 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); + TracersMonthlySurfClimo, LocTracerArray); }); Pacer::stop("Tend:surfaceTracerRestoring", 2); } @@ -957,11 +1058,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); @@ -978,7 +1079,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..3cf52d9690cc 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->MaxLayerEdgeTop), + EdgeSignOnCell(Mesh->EdgeSignOnCell), DvEdge(Mesh->DvEdge), + AreaCell(Mesh->AreaCell), + TotalVerticalPseudoVelocity(VAdv->TotalVerticalPseudoVelocity), + 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) @@ -119,6 +134,8 @@ void TracerHorzAdvOnCell::init() { const auto MaxEdges2 = Mesh->MaxEdges2; const auto NEdgesAll = Mesh->NEdgesAll; const auto NCellsAll = Mesh->NCellsAll; + const auto NCellsSize = Mesh->NCellsSize; + const auto NEdgesSize = Mesh->NEdgesSize; // Allocate Kokkos arrays in member data if (ForceLowOrder) { @@ -134,14 +151,82 @@ void TracerHorzAdvOnCell::init() { {NCellsAll}, KOKKOS_LAMBDA(int ICell) { secondDerivativeOnCell(DerivTwo, ICell); }); // Compute masks and coefficients - Kokkos::fence(); MasksAndCoefficients masksAndCoefficients( Mesh, VCoord, DerivTwo, NAdvCellsForEdge, AdvCellsForEdge, AdvMaskHighOrder, AdvCoefs, AdvCoefs3rd); - Kokkos::fence(); parallelFor( {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); }); - Kokkos::fence(); + if (FCT) { + const int NVertsFCT = NVertLayers + 1; + HProvInv = Array2DReal("FCTProvesionalLayerThickness", NCellsSize, NVertsFCT); + HNewInv = Array2DReal("FCTProvesionalNewInverse", NCellsSize, NVertsFCT); + HProv = Array2DReal("FCTProvesionalThickness", NCellsSize, NVertsFCT); + TracerCur = Array2DReal("TracerCur", NCellsSize, NVertsFCT), + TracerMax = Array2DReal("FCTTracerMax", NCellsSize, NVertsFCT); + TracerMin = Array2DReal("FCTTracerMin", NCellsSize, NVertsFCT); + HighOrderFlx = Array2DReal("FCTHighOrderFlx", NEdgesSize, NVertsFCT); + LowOrderFlx = Array2DReal("FCTLowOrderFlx", NEdgesSize, NVertsFCT); + WorkTend = Array2DReal("WorkTend", NCellsSize, NVertsFCT); + FlxIn = Array2DReal("FlxIn", NCellsSize, NVertsFCT); + FlxOut = Array2DReal("FlxOut", NCellsSize, NVertsFCT); + deepCopy(HProvInv, 0.0); + deepCopy(HNewInv, 0.0); + deepCopy(HProv, 0.0); + deepCopy(TracerCur, 0.0); + deepCopy(TracerMax, 0.0); + deepCopy(TracerMin, 0.0); + deepCopy(HighOrderFlx, 0.0); + deepCopy(LowOrderFlx, 0.0); + deepCopy(WorkTend, 0.0); + deepCopy(FlxIn, 0.0); + deepCopy(FlxOut, 0.0); + 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 548832b290bf..4723e84e9d04 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; @@ -375,86 +376,409 @@ 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); + KOKKOS_FUNCTION ~TracerHorzAdvOnCell() {} void init(); - KOKKOS_FUNCTION void operator()(const I4 L, const I4 IEdge, const I4 KChunk, + KOKKOS_FUNCTION void operator()(const I4 L, const I4 IEdge, const I4 K, const Array3DReal &TracerCell, const Array2DReal &FluxPseudoThickEdge, const Array2DReal &NormVelEdge) const { - const I4 KStart = KChunk * VecLength; - const I4 KEnd = KStart + VecLength; - for (int K = KStart; K < KEnd; ++K) - HighOrderFlxHorz(L, IEdge, K) = 0; + HighOrderFlxHorz(L, IEdge, K) = 0; // Stay at low order at boundaries - for (int K = KStart; K < KEnd; ++K) { - const I4 JCell0 = CellsOnEdge(IEdge, 0); - const I4 JCell1 = CellsOnEdge(IEdge, 1); - const Real NormalThicknessFlux = - FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); - const Real TracerWgt = DvEdge(IEdge) * 0.5_Real * NormalThicknessFlux; - HighOrderFlxHorz(L, IEdge, K) += - TracerWgt * (1._Real - AdvMaskHighOrder(IEdge, K)) * - (TracerCell(L, JCell1, K) + TracerCell(L, JCell0, K)); - } + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + const Real TracerWgt = DvEdge(IEdge) * 0.5_Real * NormalThicknessFlux; + HighOrderFlxHorz(L, IEdge, K) += + TracerWgt * (1._Real - AdvMaskHighOrder(IEdge, K)) * + (TracerCell(L, JCell1, K) + TracerCell(L, JCell0, K)); // High order (3rd or 4th) fluxes elsewhere when requested // - If HorzTracerFluxOrder = 2, NAdvCellsForEdge = 0 and // this loop is skipped. for (int I = 0; I < NAdvCellsForEdge(IEdge); ++I) { const I4 ICell = AdvCellsForEdge(IEdge, I); - for (int K = KStart; K < KEnd; ++K) { - const Real NormalThicknessFlux = - FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); - const Real TracerWgt = - (AdvCoefs(I, IEdge) + - Coef3rdOrder * std::copysign(1._Real, NormalThicknessFlux) * - AdvCoefs3rd(I, IEdge)) * - NormalThicknessFlux; - HighOrderFlxHorz(L, IEdge, K) += TracerWgt * - TracerCell(L, ICell, K) * - AdvMaskHighOrder(IEdge, K); - } + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + const Real TracerWgt = + (AdvCoefs(I, IEdge) + + Coef3rdOrder * std::copysign(1._Real, NormalThicknessFlux) * + AdvCoefs3rd(I, IEdge)) * + NormalThicknessFlux; + HighOrderFlxHorz(L, IEdge, K) += + TracerWgt * TracerCell(L, ICell, K) * AdvMaskHighOrder(IEdge, K); } } KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, const I4 L, - const I4 ICell, const I4 KChunk) const { - const I4 KStart = KChunk * VecLength; - const I4 KEnd = KStart + VecLength; + const I4 ICell, const I4 K) const { const Real InvAreaCell = 1._Real / AreaCell(ICell); - for (int K = KStart; K < KEnd; ++K) - Tend(L, ICell, K) = 0; + Tend(L, ICell, K) = 0; for (int I = 0; I < NEdgesOnCell(ICell); ++I) { const I4 IEdge = EdgesOnCell(ICell, I); - for (int K = KStart; K < KEnd; ++K) { - Tend(L, ICell, K) += EdgeSignOnCell(ICell, I) * - HighOrderFlxHorz(L, IEdge, K) * InvAreaCell; - } + Tend(L, ICell, K) += EdgeSignOnCell(ICell, I) * + HighOrderFlxHorz(L, IEdge, K) * InvAreaCell; } } - private: + KOKKOS_FUNCTION void + FCTProvisionaLayerThicknesses(const TeamMember &Team, const I4 ICell, + const Real Dt, + const Array2DReal &FluxPseudoThickEdge, + const Array2DReal &LayerThickness, + const Array2DReal &NormVelEdge) const { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + FCTProvisionaLayerThicknesses(ICell, K, Dt, FluxPseudoThickEdge, + LayerThickness, NormVelEdge); + }); + } + + KOKKOS_FUNCTION void + FCTProvisionaLayerThicknesses(const I4 ICell, const I4 K, const Real Dt, + const Array2DReal &FluxPseudoThickEdge, + const Array2DReal &LayerThickness, + const Array2DReal &NormVelEdge) const { + + const Real InvAreaCell = Dt / AreaCell(ICell); + 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 + const Real NormalThicknessFlux = + FluxPseudoThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + HProv(ICell, K) += SignedFactor * NormalThicknessFlux; + } + // New layer thickness is after horizontal and vertical + // thickness flux + HProvInv(ICell, K) = 1.0_Real / HProv(ICell, K); + HNewInv(ICell, K) = + 1.0_Real / + (HProv(ICell, K) - Dt * TotalVerticalPseudoVelocity(ICell, K) + + Dt * TotalVerticalPseudoVelocity(ICell, K + 1)); + } + + KOKKOS_FUNCTION void FCTTracerCurFill(const TeamMember &Team, const I4 L, + const I4 ICell, + const Array3DReal &TracerArray) const { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), + INNER_LAMBDA(int K) { FCTTracerCurFill(L, ICell, K, TracerArray); }); + } + + KOKKOS_FUNCTION void FCTTracerCurFill(const I4 L, const I4 ICell, const I4 K, + const Array3DReal &TracerArray) const { + TracerCur(ICell, K) = TracerArray(L, ICell, K); + TracerMin(ICell, K) = TracerArray(L, ICell, K); + TracerMax(ICell, K) = TracerArray(L, ICell, K); + } + KOKKOS_FUNCTION void FCTTracerMinMax(const TeamMember &Team, + const I4 ICell) const { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 NEdges = NEdgesOnCell(ICell); + for (I4 I = 0; I < NEdges; ++I) { + const I4 ICell2 = CellsOnCell(ICell, I); + const I4 KMin1 = Kokkos::max(KMin, MinLayerCell(ICell2)); + const I4 KMax1 = Kokkos::min(KMax, MaxLayerCell(ICell2)); + parallelForInner( + Team, Range(KMin1, KMax1), INNER_LAMBDA(int 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 TeamMember &Team, const I4 IEdge, + const Array2DReal &FluxPseudoThickEdge, + const Array2DReal &NormVelEdge) const { + + const I4 ICell1 = CellsOnEdge(IEdge, 0); + const I4 ICell2 = CellsOnEdge(IEdge, 1); + + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + HighOrderFlx(IEdge, K) = 0; + LowOrderFlx(IEdge, K) = 0; + }); + // 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 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int 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 + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int 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 FCTFluxInOut(const TeamMember &Team, const I4 ICell, + const Real Dt, + const Array2DReal &LayerThickness) const { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const Real InvAreaCell = 1._Real / AreaCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + WorkTend(ICell, K) = 0; + FlxIn(ICell, K) = 0; + FlxOut(ICell, K) = 0; + }); + // 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; + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int 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 + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int 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 TeamMember &Team, + const I4 IEdge) const { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range(KMin, KMax), + INNER_LAMBDA(int K) { FCTRescaleHighOrderFlux(IEdge, K); }); + } + KOKKOS_FUNCTION void FCTRescaleHighOrderFlux(const I4 IEdge, + const I4 K) const { + const I4 ICell1 = CellsOnEdge(IEdge, 0); + const I4 ICell2 = CellsOnEdge(IEdge, 1); + 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 TeamMember &Team, const I4 ICell, + const Real Dt, const Array2DReal &TracerArray, + const Array2DReal &LayerThickness) const { + + // Accumulate the scaled high order vertical tendencies + // and the upwind tendencies + const Real InvAreaCell1 = 1.0_Real / AreaCell(ICell); + // 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); + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + // WorkTend on RHS is upwind tendency + // WorkTend on LHS is total horiz advect tendency + WorkTend(ICell, K) += SignedFactor * HighOrderFlx(IEdge, K); + }); + } + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int 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(ICell, K) += WorkTend(ICell, K); + }); + } + + KOKKOS_FUNCTION void + FCTComputeBudgetAdvectionEdgeFlux(const TeamMember &Team, const I4 L, + const I4 IEdge) const { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + FCTComputeBudgetAdvectionEdgeFlux(L, IEdge, K); + }); + } + KOKKOS_FUNCTION void FCTComputeBudgetAdvectionEdgeFlux(const I4 L, + const I4 IEdge, + const I4 K) const { + // Save u*h*T flux on edge for analysis. This variable will be + // divided by h at the end of the time step. + ActiveTracerHorizontalAdvectionEdgeFlux(L, IEdge, K) = + (LowOrderFlx(IEdge, K) + HighOrderFlx(IEdge, K)) / DvEdge(IEdge); + } + + KOKKOS_FUNCTION void + FCTComputeBudgetAdvectionTendency(const TeamMember &Team, const I4 L, + const I4 ICell) const { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), INNER_LAMBDA(int K) { + FCTComputeBudgetAdvectionTendency(L, ICell, K); + }); + } + KOKKOS_FUNCTION void FCTComputeBudgetAdvectionTendency(const I4 L, + const I4 ICell, + const I4 K) const { + ActiveTracerHorizontalAdvectionTendency(L, ICell, K) = WorkTend(ICell, K); + } + + KOKKOS_FUNCTION void FCTMonotonicityCheck(const TeamMember &Team, + const I4 ICell) const { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + parallelForInner( + Team, Range(KMin, KMax), + INNER_LAMBDA(int K) { FCTMonotonicityCheck(ICell, K); }); + } + KOKKOS_FUNCTION void FCTMonotonicityCheck(const I4 ICell, const I4 K) const { + // Check tracer values against local min,max to detect + // non-monotone values and write warning if found + if (TracerCur(ICell, K) < TracerMin(ICell, K) - Eps) { + Kokkos::printf( + "Horizontal minimum out of bounds on cell %d, level %d " + "new tracer value %lg is smaller than previous tracer minimum " + "%lg\n", + ICell, K, TracerCur(ICell, K), TracerMin(ICell, K)); + } + if (TracerCur(ICell, K) > TracerMax(ICell, K) + Eps) { + Kokkos::printf( + "Horizontal maximum out of bounds on cell %d, level %d " + "new tracer value %lg is larger than previous tracer maximum " + "%lg\n", + ICell, K, TracerCur(ICell, K), TracerMax(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; + Array1DI4 NEdgesOnCell; 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 TotalVerticalPseudoVelocity; + 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 e7f5038be0ce..4ea8654875fc 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -38,6 +38,7 @@ #include #include +#include #include using namespace OMEGA; @@ -57,6 +58,16 @@ struct TestSetupPlane { 0.00134324628763667899}; ErrorMeasures ExpectedTrHAdvErrors = {0.0029211089892916243, 0.0024583038518548855}; + ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; + ErrorMeasures ExpectedFCTHProv = {0.00, 0.00}; + ErrorMeasures ExpectedFCTHInv = {0.00, 0.00}; + ErrorMeasures ExpectedFCTHNew = {0.00, 0.00}; + ErrorMeasures ExpectedFCT_High = {1.0408340855860843e-17, + 1.2836244918860014e-17}; + ErrorMeasures ExpectedFCT_Low = {1.0095554103513882e-15, + 1.0095554103514163e-15}; + ErrorMeasures ExpectedFCTToNonFCT = {9.524542988823467e-16, + 5.336197203674528e-16}; ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, 0.00290978146207349032}; ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, @@ -66,6 +77,10 @@ struct TestSetupPlane { ErrorMeasures ExpectedBottomDragErrors = {0.033848740052302935, 0.01000133508329411}; + KOKKOS_FUNCTION Real velEdgeX(Real X, Real Y) const { return 2 / Lx; } + + KOKKOS_FUNCTION Real velEdgeY(Real X, Real Y) const { return 3 / Ly; } + KOKKOS_FUNCTION Real vectorX(Real X, Real Y) const { return std::sin(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } @@ -187,16 +202,27 @@ struct TestSetupSphere { // TODO: get this from the mesh Real Radius = REarth; - ErrorMeasures ExpectedDivErrors = {0.013652414501664885, - 0.0036904315983599676}; - ErrorMeasures ExpectedPVErrors = {0.0219217796608757037, - 0.0122537418367830303}; - ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, - 0.0015218320661105702}; - ErrorMeasures ExpectedLaplaceErrors = {0.28193638497826856, - 0.270546491554748}; - ErrorMeasures ExpectedTrHAdvErrors = {0.013259410329645643, - 0.004094907022292395}; + ErrorMeasures ExpectedDivErrors = {0.013652414501664885, + 0.0036904315983599676}; + ErrorMeasures ExpectedPVErrors = {0.0219217796608757037, + 0.0122537418367830303}; + ErrorMeasures ExpectedGradErrors = {0.0019094381714837498, + 0.0015218320661105702}; + ErrorMeasures ExpectedLaplaceErrors = {0.28193638497826856, + 0.270546491554748}; + ErrorMeasures ExpectedTrHAdvErrors = {0.013259410329645643, + 0.004094907022292395}; + ErrorMeasures ExpectedFCTErrors = {0.00, 0.00}; + ErrorMeasures ExpectedFCTHProv = {3.0542657508680904e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCTHInv = {3.0541724683419424e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCTHNew = {3.0541724683419424e-05, + 1.0779233406323438e-06}; + ErrorMeasures ExpectedFCT_High = {0.00146484375, 16478.526025524854}; + ErrorMeasures ExpectedFCT_Low = {0.2958192441638326, 0.09980000200239414}; + ErrorMeasures ExpectedFCTToNonFCT = {1.5120472343676013e-15, + 6.215666716489436e-16}; ErrorMeasures ExpectedTrDel2Errors = {0.04865718541236144, 0.005105510870642706}; ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, @@ -206,6 +232,10 @@ struct TestSetupSphere { ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, 0.0014897009917655022}; + KOKKOS_FUNCTION Real velEdgeX(Real X, Real Y) const { return 2; } + + KOKKOS_FUNCTION Real velEdgeY(Real X, Real Y) const { return 3; } + KOKKOS_FUNCTION Real vectorX(Real Lon, Real Lat) const { return -Radius * std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); } @@ -1075,7 +1105,6 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -1111,7 +1140,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(); @@ -1138,6 +1171,408 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { LOG_INFO("TendencyTermsTest: TracerHorzAdv PASS"); } + VertAdv::clear(); + return Err; +} // end testTracerHorzAdvOnCell + +class TracerHorzAdvOnCellTest : public TracerHorzAdvOnCell { + public: + TracerHorzAdvOnCellTest(const HorzMesh *Mesh, const VertCoord *VCoord, + const VertAdv *VAdv) + : TracerHorzAdvOnCell(Mesh, VCoord, VAdv) {}; + Array2DReal GetHProvInv() const { return HProvInv; }; + Array2DReal GetHProv() const { return HProv; }; + Array2DReal GetHNewInv() const { return HNewInv; }; + Array2DReal GetTracerCur() const { return TracerCur; }; + Array2DReal GetTracerMin() const { return TracerMin; }; + Array2DReal GetTracerMax() const { return TracerMax; }; + Array3DReal GetHighOrderFlxHorz() const { return HighOrderFlxHorz; }; + Array2DReal GetLowOrderFlx() const { return LowOrderFlx; }; + Array2DReal GetHighOrderFlx() const { return HighOrderFlx; }; + Array2DReal GetFlxIn() const { return FlxIn; }; + Array2DReal GetFlxOut() const { return FlxOut; }; +}; + +int testFCTTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) { + + I4 Err = 0; + TestSetup Setup; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + const int NCellsFCT = Mesh->NCellsSize; + const int NEdgesFCT = Mesh->NEdgesSize; + const int NVertsFCT = NVertLayers + 1; + + // Compute exact result + Array3DReal ExactTrFluxDiv("ExactTrFluxDiv", NTracers, Mesh->NCellsOwned, + NVertsFCT); + + 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", NEdgesFCT, NVertsFCT); + + Array3DReal TrCell("TrCell", NTracers, NCellsFCT, NVertsFCT); + Array2DReal Tend("Tend", NCellsFCT, NVertsFCT); + Array3DReal TendNoFCT("TendNoFCT", NTracers, NCellsFCT, NVertsFCT); + Array2DReal FluxPseudoThickEdge("FluxPseudoThickEdge", NEdgesFCT, NVertsFCT); + Array2DReal LayerThickness("LayerThickness", NCellsFCT, NVertsFCT); + Array2DReal NormVelEdge("NormVelEdge", NEdgesFCT, NVertsFCT); + Array2DReal HProvInvExact("HProvInvExact", NCellsFCT, NVertsFCT); + Array2DReal HProvExact("HProvExact", NCellsFCT, NVertsFCT); + Array2DReal HNewInvExact("HProvExact", NCellsFCT, NVertsFCT); + Array2DReal TracerSubView("TracerSubView", NCellsFCT, NVertsFCT); + Array2DReal CellSubView("CellSubView", NCellsFCT, NVertsFCT); + Array2DReal FluxSubView("FluxSubView", NEdgesFCT, NVertsFCT); + Array2DReal ReferenceSolution("ReferenceSolution", NEdgesFCT, NVertsFCT); + deepCopy(TrCell, 1._Real); + deepCopy(Tend, 0._Real); + deepCopy(TendNoFCT, 0._Real); + deepCopy(Kokkos::subview(TrCell, Kokkos::ALL, Kokkos::ALL, NVertLayers), + 0._Real); + deepCopy(FluxPseudoThickEdge, 1._Real); + deepCopy(LayerThickness, 1._Real); + deepCopy(HProvExact, 1._Real); + deepCopy(Kokkos::subview(HProvExact, Kokkos::ALL, NVertLayers), 0._Real); + deepCopy(HProvInvExact, 1._Real); + deepCopy(Kokkos::subview(HProvInvExact, Kokkos::ALL, NVertLayers), 0._Real); + deepCopy(HNewInvExact, 1._Real); + deepCopy(Kokkos::subview(HNewInvExact, Kokkos::ALL, NVertLayers), 0._Real); + deepCopy(NormVelEdge, 1._Real); + deepCopy(ReferenceSolution, 0._Real); + deepCopy(CellSubView, 0._Real); + + // Compute numerical result + Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned, + NVertsFCT); + + VertAdv::init(); + const auto VAdv = VertAdv::getDefault(); + deepCopy(VAdv->VerticalPseudoVelocity, 10._Real); + deepCopy(VAdv->TotalVerticalPseudoVelocity, 10._Real); + + TracerHorzAdvOnCellTest TrHorzAdvOnC(Mesh, VCoord, VAdv); + TrHorzAdvOnC.ForceLowOrder = false; + TrHorzAdvOnC.FCT = true; + TrHorzAdvOnC.init(); + + const Real Dt = 1; + ErrorMeasures FCTErrors; + + { + setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.velEdgeX(X, Y); + VecField[1] = Setup.velEdgeY(X, Y); + }, + NormalVelocity, EdgeComponent::Normal, Geom, Mesh); + + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, K, Dt, FluxPseudoThickEdge, LayerThickness, + NormalVelocity); + }); + Kokkos::fence(); + + const auto HProvInv = TrHorzAdvOnC.GetHProvInv(); + const auto HProv = TrHorzAdvOnC.GetHProv(); + const auto HNewInv = TrHorzAdvOnC.GetHNewInv(); + + const Real ATol = 1.0e-10; + Err = 0; + Err += computeErrors(FCTErrors, HProv, HProvExact, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", "FCTHProv", FCTErrors, + Setup.ExpectedFCTHProv, RTol, ATol); + + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTHProv PASS"); + + Err += computeErrors(FCTErrors, HProvInv, HProvInvExact, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", "FCTHProvInv", FCTErrors, + Setup.ExpectedFCTHInv, RTol, ATol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTHProvInv PASS"); + Err = 0; + Err += computeErrors(FCTErrors, HNewInv, HNewInvExact, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", "FCTHNewInv", FCTErrors, + Setup.ExpectedFCTHNew, RTol, ATol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTHNewInv PASS"); + Err = 0; + } + for (int L = 0; L < NTracers; ++L) { + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrCell(L, ICell, K) = L + ICell + K; + }); + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, K, TrCell); + }); + const Array2DReal TracerCur = TrHorzAdvOnC.GetTracerCur(); + const Array2DReal Tracers = + Kokkos::subview(TrCell, L, Kokkos::ALL, Kokkos::ALL); + deepCopy(TracerSubView, Tracers); + Err += computeErrors(FCTErrors, TracerCur, TracerSubView, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", "FCTTracerCurFill", FCTErrors, + Setup.ExpectedFCTErrors, RTol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTTracerCurFill_" + std::to_string(L) + + " PASS"); + Err = 0; + } + { + setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.velEdgeX(X, Y); + VecField[1] = Setup.velEdgeY(X, Y); + }, + NormalVelocity, EdgeComponent::Normal, Geom, Mesh); + deepCopy(TrCell, 1._Real); + deepCopy(Kokkos::subview(TrCell, Kokkos::ALL, Kokkos::ALL, NVertLayers), + 0._Real); + const int L = 0; + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, K, TrCell); + }); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, K, Dt, FluxPseudoThickEdge, LayerThickness, + NormalVelocity); + }); + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(const int IEdge, const TeamMember &Team) { + TrHorzAdvOnC.FCTHighAndLowOrderFlux( + Team, IEdge, FluxPseudoThickEdge, NormalVelocity); + }); + Kokkos::fence(); + + Array2DReal LowOrderFlx = TrHorzAdvOnC.GetLowOrderFlx(); + + parallelFor( + {Mesh->NEdgesHaloH(1), NVertLayers}, + KOKKOS_LAMBDA(const int IEdge, const int K) { + const double DvEdgeFromMeshFile = + Geom == Geometry::Planar ? 0.012028130608117 : 277371; + FluxSubView(IEdge, K) = + DvEdgeFromMeshFile * NormalVelocity(IEdge, K); + }); + + const Real ATol = 1.0e-10; + Err += computeErrors(FCTErrors, LowOrderFlx, FluxSubView, Mesh, OnEdge); + Err += checkErrors("TendencyTermsTest", "FCTHighAndLowOrderFlux_Low", + FCTErrors, Setup.ExpectedFCT_Low, RTol, ATol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTHighAndLowOrderFlux_Low PASS"); + Err = 0; + } + if (Geom == Geometry::Planar) { + 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); + + deepCopy(TrCell, 1._Real); + deepCopy(Kokkos::subview(TrCell, Kokkos::ALL, Kokkos::ALL, NVertLayers), + 0._Real); + const int L = 0; + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, K, TrCell); + }); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, K, Dt, FluxPseudoThickEdge, LayerThickness, + NormalVelocity); + }); + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(const int IEdge, const TeamMember &Team) { + TrHorzAdvOnC.FCTHighAndLowOrderFlux( + Team, IEdge, FluxPseudoThickEdge, NormalVelocity); + }); + Kokkos::fence(); + + parallelFor( + {Mesh->NEdgesHaloH(1), NVertLayers}, + KOKKOS_LAMBDA(const int IEdge, const int K) { + FluxSubView(IEdge, K) = 0; + }); + + const Real ATol = 1.0e-10; + Array2DReal HighOrderFlx = TrHorzAdvOnC.GetHighOrderFlx(); + Err += computeErrors(FCTErrors, HighOrderFlx, FluxSubView, Mesh, OnEdge); + Err += checkErrors("TendencyTermsTest", "FCTHighAndLowOrderFlux_High", + FCTErrors, Setup.ExpectedFCT_High, RTol, ATol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTHighAndLowOrderFlux_High PASS"); + Err = 0; + + Kokkos::fence(); + } + if (Geom == Geometry::Planar) { + setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { + VecField[0] = Setup.velEdgeX(X, Y); + VecField[1] = Setup.velEdgeY(X, Y); + }, + NormalVelocity, EdgeComponent::Normal, Geom, Mesh); + + deepCopy(TrCell, 1._Real); + deepCopy(Kokkos::subview(TrCell, Kokkos::ALL, Kokkos::ALL, NVertLayers), + 0._Real); + const int L = 0; + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, K, TrCell); + }); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, K, Dt, FluxPseudoThickEdge, LayerThickness, + NormalVelocity); + }); + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(const int IEdge, const TeamMember &Team) { + TrHorzAdvOnC.FCTHighAndLowOrderFlux( + Team, IEdge, FluxPseudoThickEdge, NormalVelocity); + }); + parallelForOuter( + {Mesh->NCellsHaloH(0)}, + KOKKOS_LAMBDA(const int ICell, const TeamMember &Team) { + TrHorzAdvOnC.FCTFluxInOut(Team, ICell, Dt, LayerThickness); + }); + + Kokkos::fence(); + Err = 0; + Array2DReal FlxOut = TrHorzAdvOnC.GetFlxOut(); + Err += computeErrors(FCTErrors, FlxOut, CellSubView, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", "FCTFluxOut", FCTErrors, + Setup.ExpectedFCTErrors, RTol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTFluxOut PASS"); + Err = 0; + Kokkos::fence(); + } + { + parallelFor( + {Mesh->NEdgesHaloH(0), NVertLayers}, + KOKKOS_LAMBDA(const int IEdge, const int K) { + TrHorzAdvOnC.FCTRescaleHighOrderFlux(IEdge, K); + }); + Kokkos::fence(); + const Array2DReal HighOrderFlx = TrHorzAdvOnC.GetHighOrderFlx(); + Err += computeErrors(FCTErrors, HighOrderFlx, ReferenceSolution, Mesh, + OnEdge); + Err += checkErrors("TendencyTermsTest", "FCTRescaleHighOrderFlux", + FCTErrors, Setup.ExpectedFCTErrors, RTol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTRescaleHighOrderFlux PASS"); + Err = 0; + Kokkos::fence(); + } + { + 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); + + deepCopy(TrCell, 1._Real); + deepCopy(Kokkos::subview(TrCell, Kokkos::ALL, Kokkos::ALL, NVertLayers), + 0._Real); + const int L = 0; + parallelFor( + {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTTracerCurFill(L, ICell, K, TrCell); + }); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(const int ICell, const int K) { + TrHorzAdvOnC.FCTProvisionaLayerThicknesses( + ICell, K, Dt, FluxPseudoThickEdge, LayerThickness, + NormalVelocity); + }); + parallelForOuter( + {Mesh->NEdgesHaloH(1)}, + KOKKOS_LAMBDA(const int IEdge, const TeamMember &Team) { + TrHorzAdvOnC.FCTHighAndLowOrderFlux( + Team, IEdge, FluxPseudoThickEdge, NormalVelocity); + }); + parallelForOuter( + {Mesh->NCellsHaloH(0)}, + KOKKOS_LAMBDA(const int ICell, const TeamMember &Team) { + TrHorzAdvOnC.FCTFluxInOut(Team, ICell, Dt, LayerThickness); + }); + parallelForOuter( + {Mesh->NCellsOwned}, + KOKKOS_LAMBDA(const int ICell, const TeamMember &Team) { + TrHorzAdvOnC.FCTAccumulateHighOrderFlux(Team, ICell, Dt, Tend, + LayerThickness); + }); + Kokkos::fence(); + + const auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + const auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + parallelForOuter( + "Tend:TracerHorzAdv.init", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + TrHorzAdvOnC(L, IEdge, K, TrCell, FluxPseudoThickEdge, + NormalVelocity); + }); + }); + parallelForOuter( + "Tend:TracerHorzAdv.exec", {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) { TrHorzAdvOnC(TendNoFCT, L, ICell, K); }); + }); + deepCopy(TracerSubView, + Kokkos::subview(TendNoFCT, L, Kokkos::ALL, Kokkos::ALL)); + + const Real ATol = 1.0e-10; + Err += computeErrors(FCTErrors, Tend, TracerSubView, Mesh, OnCell); + Err += checkErrors("TendencyTermsTest", + "FCTAccumulateHighOrderFlux_" + std::to_string(L), + FCTErrors, Setup.ExpectedFCTToNonFCT, RTol, ATol); + if (Err == 0) + LOG_INFO("TendencyTermsTest: FCTAccumulateHighOrderFlux_" + + std::to_string(L) + " PASS"); + Err = 0; + Kokkos::fence(); + } + + VertAdv::clear(); return Err; } // end testTracerHorzAdvOnCell @@ -1448,6 +1883,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);