diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index d592acd1359a..f7e408e38360 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -44,6 +44,7 @@ Omega: BottomDragTendencyEnable: false BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true + TracerHorzAdvTendencyOrder: 1 TracerDiffTendencyEnable: true EddyDiff2: 10.0 TracerHyperDiffTendencyEnable: true diff --git a/components/omega/doc/design/Tendencies.md b/components/omega/doc/design/Tendencies.md index 9bceff5a0b38..3f392698c12b 100644 --- a/components/omega/doc/design/Tendencies.md +++ b/components/omega/doc/design/Tendencies.md @@ -39,14 +39,23 @@ retrieval. ```c++ class Tendencies{ public: + // Arrays for accumulating tendencies Array2DReal LayerThicknessTend; Array2DReal NormalVelocityTend; + Array3DReal TracerTend; + // Instances of tendency terms ThicknessFluxDivOnCell ThicknessFluxDiv; PotentialVortHAdvOnEdge PotientialVortHAdv; KEGradOnEdge KEGrad; SSHGradOnEdge SSHGrad; VelocityDiffusionOnEdge VelocityDiffusion; VelocityHyperDiffOnEdge VelocityHyperDiff; + WindForcingOnEdge WindForcing; + BottomDragOnEdge BottomDrag; + TracerDiffOnCell TracerDiffusion; + TracerHyperDiffOnCell TracerHyperDiff; + TracerHorzAdvOnCell TracerHorzAdv; + TracerHighOrderHorzAdvOnCell TracerHighOrderHorzAdv; private: static Tendencies *DefaultTendencies; static std::map> AllTendencies; @@ -79,29 +88,29 @@ Tendencies* Tendencies::create(const std::string &Name, const HorzMesh *Mesh, in #### 4.2.2 Initialization The init method will create the default tendencies and return an error code: ```c++ -int Tendencies::init(); +static void Tendencies::init(); ``` #### 4.2.3 Retrieval There will be methods for getting the default and non-default tendencies instances: ```c++ -Tendencies *Tendencies::getDefault(); -Tendencies *Tendencies::get(const std::string &Name); +static Tendencies *Tendencies::getDefault(); +static Tendencies *Tendencies::get(const std::string &Name); ``` #### 4.2.4 Computation The 'computeAll' method will compute and accumulate all layer thickness and normal velocity tendencies using the ocean state and auxiliary state from a given time level: ```c++ -void Tendencies::computeAllTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel); +void Tendencies::computeAllTendencies(const OceanState *State, const AuxilaryState *AuxState, const Array3DReal &TracerArray, int TimeLevel, int VelTimeLevel, TimeInstant Time); ``` The layer thickness tendencies will be computed with a method: ```c++ -void Tendencies::computeThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel); +void Tendencies::computeThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time); ``` The normal velocity tendencies will be computed with a method: ```c++ -void Tendencies::computeVelocityTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel); +void Tendencies::computeVelocityTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time); ``` @@ -111,8 +120,8 @@ when they are no longer in scope. The erase method will remove a named tendencies instance, whereas the clear method will remove all of them. Both will call the destructor in the process. ```c++ -void Tendencies::erase(const std::string &Name); -void Tendencies::clear(); +static void Tendencies::erase(const std::string &Name); +static void Tendencies::clear(); ``` ## 5 Verification and Testing diff --git a/components/omega/doc/design/Tendency.md b/components/omega/doc/design/Tendency.md index a4f077638ce9..da797e057cf4 100644 --- a/components/omega/doc/design/Tendency.md +++ b/components/omega/doc/design/Tendency.md @@ -39,24 +39,25 @@ The class contains private member variables for any constant data that are defin Each tendency term will have a `bool` variable which can be set to enable/disable the computation of the tendency. ```c++ -class ThicknessFluxDivergenceOnCell { +class ThicknessFluxDivOnCell { public: bool enabled = false; - ThicknessFluxDivergenceOnCell(const HorzMesh *Mesh, Config *Options); + ThicknessFluxDivOnCell(const HorzMesh *Mesh, Config *Options); KOKKOS_FUNCTION void operator()(Array2DReal &Tend int ICell, int KChunk, - const Array2DReal &ThicknessFlux); + const Array2DReal &ThicknessFlux, + const Array2DReal &NormalVelEdge); private: Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; - Array1DR8 DvEdge; - Array1DR8 AreaCell; - Array2DR8 EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal AreaCell; + Array2DReal EdgeSignOnCell; }; ``` @@ -66,13 +67,10 @@ class ThicknessFluxDivergenceOnCell { Tendency functor constructors are responsible for initializing the private member variables: ```c++ -ThicknessFluxDivergenceOnCell::ThicknessFluxDivergenceOnCell(HorzMesh const *Mesh, Config *Options) +ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(HorzMesh const *Mesh) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), - EdgeSignOnCell(Mesh->EdgeSignOnCell) { - - enabled = Options->get("ThicknessFluxTendencyEnable") -} + EdgeSignOnCell(Mesh->EdgeSignOnCell) { } ``` #### 4.2.2 operator @@ -80,17 +78,27 @@ The operator method implements the tendency computation for a chunk of vertical The inner loop over a chunk of vertical layers enables CPU vectorization. ```c++ -KOKKOS_FUNCTION void ThicknessFluxDivergenceOnCell::operator()(Array2DReal &Tend - int ICell, - int KChunk, - const Array2DReal &ThicknessFlux) const { - const Real InvAreaCell = 1. / AreaCell(iCell); - for (int J = 0; J < NEdgesOnCell(ICell); ++J) { - const int JEdge = EdgesOnCell(ICell, J); - for (int K = KChunk * VecLength; K < (KChunk + 1) * VecLength; ++K) { - Tend(ICell, K) -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * ThicknessFlux(JEdge, K) * InvAreaCell; - } - } +KOKKOS_FUNCTION void ThicknessFluxDivOnCell::operator()(Array2DReal &Tend + int ICell, + int KChunk, + const Array2DReal &ThicknessFlux, + const Array2DReal &NormalVelEdge) const { + const I4 KStart = KChunk * VecLength; + const Real InvAreaCell = 1._Real / AreaCell(ICell); + Real DivTmp[VecLength] = {0}; + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + const I4 JEdge = EdgesOnCell(ICell, J); + for (int KVec = 0; KVec < VecLength; ++KVec) { + const I4 K = KStart + KVec; + DivTmp[KVec] -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * + ThicknessFlux(JEdge, K) * NormalVelEdge(JEdge, K) * + InvAreaCell; + } + } + for (int KVec = 0; KVec < VecLength; ++KVec) { + const I4 K = KStart + KVec; + Tend(ICell, K) -= DivTmp[KVec]; + } } ``` diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 536efbffb60c..e665fb03ed76 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -35,8 +35,9 @@ implemented: - `SSHGradOnEdge` - `VelocityDiffusionOnEdge` - `VelocityHyperDiffOnEdge` +- `WindForcingOnEdge` +- `BottomDragOnEdge` - `TracerHorzAdvOnCell` +- `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` -- `WindForcingOnEdge` -- `BottomDragOnEdge` diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 69ad8437d210..369d5673fd2d 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -15,6 +15,7 @@ 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 | WindForcingOnEdge | forcing by wind stress, defined on edges @@ -43,6 +44,9 @@ 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 +| | TracerHorzAdvTendencyOrder | 1 for standard linear advection +| TracerHighOrderHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term +| | TracerHorzAdvTendencyOrder | 2 for second order advection algorithm | TracerDiffOnCell | TracerDiffTendencyEnable | enable/disable term | | EddyDiff2 | horizontal diffusion coefficient | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term @@ -50,3 +54,85 @@ the currently available tendency terms: | WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term | | BottomDragCoeff | bottom drag coefficient + + +## Second Order Horizontal Advection Algorithm + +The horizontal advection is done independently within each ocean layer +so a two dimensional scheme is required to advect mixing ratios. +The second order horizontal advection scheme is described in the paper + +William C. Skamarock and Almut Gassmann, +"Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration" +Monthly Weather Review, Vol 139: Issue 9, pp 2962–2975, 2011. doi.org/10.1175/MWR-D-10-05056.1 + +and only a summary of a few equations from that paper are reproduced here. + +The conservative form of the scalar transport equation +within a fluid is given in equation (1) of Skamarock and Gassmann, +$$ + \frac{\partial(\rho\psi)}{\partial t} = -\nabla\cdot\mathbf{V}\rho\psi +$$ +where $\rho$ is the fluid density, $\psi$ is the mixing ratio, and $\mathbf{V}$ is +the fluid velocity. This is simplified by assuming that $\rho$ is constant and cancels out. +The advection scheme for a finite volume method involves approximating this equation along an element +edge at a certain time step and multiplying by the time step length to get an approximate scalar transport +from element to element during the time step. +The left hand side of this equation is a time derivative and for Omega this is handled by standard Runge—Kutta methods and +not discussed here but can be found in the Skamarock and Gassmann paper. +The right hand side of this equation is the spatial derivative that needs to be evaluated by the Omega horizontal advection scheme. +Omega ocean uses unstructured Voronoi meshes and a +finite volume scheme with $\psi$ defined on cell centers. A first order approximation +to the spatial derivative of the left hand termacross a single edge of an element would be the difference of $\psi$ at +element centroids divided by the distance from cell centroid to cell centroid. This is expressed in equation (5) of the paper. +For edge $i$ between elements $i+1/2$ and $i-1/2$, +$$ + \frac{\partial(u\psi_i)}{\partial x} = \frac{1}{\Delta x}[F_{i+1/2}(u\psi) - F_{i-1/2}(u\psi)] + O(\Delta x^2). +$$ +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 TracerHorzAdvTendencyOrder is 1. +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 TracerHorzAdvTendencyOrder set to 2, 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, +$$ +defined for each edge in the mesh. The six coefficients in this quadratic equation are calculated from +a least squares fit of $\psi$ in a neighborhood of elements around an edge. The +least squares fit is defined in equation (12) of the paper, +$$ + \mathbf{f = (P^T P)^{-1}P^T s = Bs} +$$ +where the vector $\mathbf{f} = [c_0,c_x,c_y,c_{xx},c_{xy},c_{yy}]$ and $\mathbf{P}$ is an $m\times 6$ matrix with each row, +$(1,x_i,y_i,x_i^2,x_iy_i,y_i^2)$, based on the cell center coordinates where there is one row for each of the $m$ elements in the neighborhood. +This results in a $6\times m$ matrix $\mathbf{B}$. The values of $\mathbf{s}$ are the mixing ratios, +$\mathbf{s}=[\psi_0, \psi_1,...,\psi_m]$, of the elements in the neighborhood at the time step being evaluated. +Note that $\mathbf{B}$ depends on the mesh geometry only, so only needs to +be computed once for each edge during initialization. Since only the second order terms are used to determine the +second order derivatives needed for the higher order flux corrections to the lower order linear flux given above +and the derivatives are taken along a line connecting +cell centers, the amount of data that needs to be stored is minimized and $\mathbf{B}$ reduces to a vector. + +The second order derivatives of $\psi$ are then combined as shown in equation (11) of Skamarock and Gassmann to get +a second order horizontal transport algorithm. + +The computation of these second order fluxes is complicated by having to do these computations on a sphere. +Figure 2 from the Skamarock and Gassmann paper describs the slight modifications needed to compute the entries +of the $\mathbf{P}$ matrix for a spherical geometry where the distances need to be measured as arc lengths along +the surface. + +In practice, the convergence of this higher order algorithm on a sphere is slightly below +second order. For the advection of a $\psi$ in the shape of a cosine bell +being advected around a sphere and compared with the analytic solution, the spatial convergence for grid refinement +is about order 1.7 as shown in {numref}`tracer-higher-order-convergence`: + +```{figure} images/higher_order_tracer_convergence_on_sphere.jpeg +:name: tracer-higher-order-convergence +:align: center +:width: 600 px +Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate +``` diff --git a/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere.jpeg b/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere.jpeg new file mode 100644 index 000000000000..fdcd0424a4c2 Binary files /dev/null and b/components/omega/doc/userGuide/images/higher_order_tracer_convergence_on_sphere.jpeg differ diff --git a/components/omega/src/infra/Config.cpp b/components/omega/src/infra/Config.cpp index 6e7077a60062..493ddc6eb744 100644 --- a/components/omega/src/infra/Config.cpp +++ b/components/omega/src/infra/Config.cpp @@ -33,17 +33,7 @@ int Config::NumReadGroups = 0; MPI_Comm Config::ConfigComm; Config Config::ConfigAll; -//------------------------------------------------------------------------------ -// Constructor that creates configuration with a given name and an emtpy -// YAML node that will be filled later with either a readAll or get call. -Config::Config(const std::string &InName // [in] name of config, node -) { - - // If this is the first Config created, set variables for reading the - // input configuration file. In particular, to avoid too many MPI tasks - // reading the same input stream, we divide the tasks into groups and - // only one group at a time loads the full configuration from a stream - +void Config::Initialize() { if (NotInitialized) { // Determine MPI variables MachEnv *DefEnv = MachEnv::getDefault(); @@ -58,6 +48,18 @@ Config::Config(const std::string &InName // [in] name of config, node NotInitialized = false; // now initialized for future calls } +} +//------------------------------------------------------------------------------ +// Constructor that creates configuration with a given name and an emtpy +// YAML node that will be filled later with either a readAll or get call. +Config::Config(const std::string &InName // [in] name of config, node +) { + + // If this is the first Config created, set variables for reading the + // input configuration file. In particular, to avoid too many MPI tasks + // reading the same input stream, we divide the tasks into groups and + // only one group at a time loads the full configuration from a stream + Initialize(); // Set the name - this is also used as the name of the root YAML node // (a map node) in this configuration. diff --git a/components/omega/src/infra/Config.h b/components/omega/src/infra/Config.h index 54cd24dc878e..7dcb3b2f510e 100644 --- a/components/omega/src/infra/Config.h +++ b/components/omega/src/infra/Config.h @@ -59,6 +59,8 @@ class Config { public: // Methods + static void Initialize(); + /// Constructor that creates a Configuration with a given name and an /// empty YAML::Node. The node will be filled later with either a read /// or get operation. diff --git a/components/omega/src/ocn/CustomTendencyTerms.h b/components/omega/src/ocn/CustomTendencyTerms.h index 7ae226d2a27a..0635e0483ed9 100644 --- a/components/omega/src/ocn/CustomTendencyTerms.h +++ b/components/omega/src/ocn/CustomTendencyTerms.h @@ -22,8 +22,8 @@ // //===----------------------------------------------------------------------===// -#include "HorzMesh.h" -#include "TendencyTerms.h" +#include "AuxiliaryState.h" +#include "OceanState.h" #include "TimeMgr.h" namespace OMEGA { diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 0c626b49905e..c3fd1b414f08 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -17,7 +17,6 @@ #include "Error.h" #include "IO.h" #include "Logging.h" -#include "MachEnv.h" #include "OmegaKokkos.h" namespace OMEGA { @@ -46,19 +45,21 @@ void HorzMesh::init() { HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh Decomp *MeshDecomp //< [in] Decomp for the new mesh -) { - + ) + : CellID(MeshDecomp->CellID) ///< global cell ID for each local cell +{ MeshName = Name; // Retrieve mesh files name from Decomp MeshFileName = MeshDecomp->MeshFileName; // Retrieve mesh cell/edge/vertex totals from Decomp - NCellsHalo = MeshDecomp->NCellsHalo; - NCellsHaloH = MeshDecomp->NCellsHaloH; - NCellsOwned = MeshDecomp->NCellsOwned; - NCellsAll = MeshDecomp->NCellsAll; - NCellsSize = MeshDecomp->NCellsSize; + NCellsHalo = MeshDecomp->NCellsHalo; + NCellsHaloH = MeshDecomp->NCellsHaloH; + NCellsOwned = MeshDecomp->NCellsOwned; + NCellsAll = MeshDecomp->NCellsAll; + NCellsSize = MeshDecomp->NCellsSize; + NCellsGlobal = MeshDecomp->NCellsGlobal; NEdgesHalo = MeshDecomp->NEdgesHalo; NEdgesHaloH = MeshDecomp->NEdgesHaloH; @@ -67,6 +68,7 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh NEdgesSize = MeshDecomp->NEdgesSize; MaxCellsOnEdge = MeshDecomp->MaxCellsOnEdge; MaxEdges = MeshDecomp->MaxEdges; + NEdgesGlobal = MeshDecomp->NEdgesGlobal; NVerticesHalo = MeshDecomp->NVerticesHalo; NVerticesHaloH = MeshDecomp->NVerticesHaloH; diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index 0b244c0e17a9..ff824479cf3d 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -98,6 +98,7 @@ class HorzMesh { I4 NCellsOwned; ///< Number of cells owned by this task I4 NCellsAll; ///< Total number of local cells (owned + all halo) I4 NCellsSize; ///< Array size (incl padding, bndy cell) for cell arrays + I4 NCellsGlobal; Array1DI4 NEdgesHalo; ///< num cells owned+halo for halo layer HostArray1DI4 NEdgesHaloH; ///< num cells owned+halo for halo layer @@ -106,7 +107,8 @@ class HorzMesh { I4 NEdgesSize; ///< Array length (incl padding, bndy) for edge dim I4 MaxCellsOnEdge; ///< Max number of cells sharing an edge I4 MaxEdges; ///< Max number of edges around a cell - I4 MaxEdges2; ///< Max number of edges around a cell x2 + I4 NEdgesGlobal; + I4 MaxEdges2; ///< Max number of edges around a cell x2 Array1DI4 NVerticesHalo; ///< num cells owned+halo for halo layer HostArray1DI4 NVerticesHaloH; ///< num cells owned+halo for halo layer @@ -147,6 +149,7 @@ class HorzMesh { Array2DI4 EdgesOnVertex; ///< Indx of edges sharing vertex as endpoint HostArray2DI4 EdgesOnVertexH; ///< Indx of edges sharing vertex as endpoint + Array1DI4 CellID; ///< global cell ID for each local cell // Coordinates Array1DReal XCell; ///< X Coordinates of cell centers (m) diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index b0f13dbe8d59..c6ce926f4bca 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -27,4 +27,47 @@ InterpCellToEdge::InterpCellToEdge(const HorzMesh *Mesh) KiteAreasOnVertex(Mesh->KiteAreasOnVertex), VertexDegree(Mesh->VertexDegree) {} +SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh) + : OnSphere(true), NCellsAll(Mesh->NCellsAll), MaxEdges(1 + Mesh->MaxEdges), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + CellsOnCell(Mesh->CellsOnCell), CellsOnEdge(Mesh->CellsOnEdge), + VerticesOnEdge(Mesh->VerticesOnEdge), + + XCell(Mesh->XCell), YCell(Mesh->YCell), + ZCell(createDeviceMirrorCopy(Mesh->ZCellH)), DvEdge(Mesh->DvEdge), + DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), + AreaCell(Mesh->AreaCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), + XVertex(createDeviceMirrorCopy(Mesh->XVertexH)), + YVertex(createDeviceMirrorCopy(Mesh->YVertexH)), + ZVertex(createDeviceMirrorCopy(Mesh->ZVertexH)), + + ThetaAbs("SphereAngle", NCellsAll), XPCell("XP", NCellsAll, MaxEdges), + YPCell("YP", NCellsAll, MaxEdges), + Angle2DCell("Angle2D", NCellsAll, MaxEdges), + BCell("WorkSpaceForLeastSquares", NCellsAll, 6, MaxEdges), + CellListCell("CellList", NCellsAll, MaxEdges) { + if (MaxMaxEdges <= Mesh->MaxEdges) + LOG_CRITICAL( + "SecondDerivativeOnCell::SecondDerivativeOnCell Max Edges exceeded:" + "Max Allowed: {} Found in Mesh:{}", + MaxMaxEdges - 1, Mesh->MaxEdges); +} + +MasksAndCoefficients::MasksAndCoefficients( + HorzMesh const *Mesh, const Array3DReal DerivTwo, + Array1DI4 NAdvCellsForEdge, Array2DI4 AdvCellsForEdge, + Array1DI4 AdvMaskHighOrder, Array2DReal AdvCoefs, Array2DReal AdvCoefs3rd) + : NCellsGlobal(Mesh->NCellsGlobal), NCellsAll(Mesh->NCellsAll), + NAdvCellsMax(Mesh->MaxEdges2), // PatchCellLists("PatchCellLists", + // Mesh->NEdgesOwned, Mesh->NEdgesAll+1), + NAdvCellsForEdge(NAdvCellsForEdge), AdvCellsForEdge(AdvCellsForEdge), + NEdgesOnEdge(Mesh->NEdgesOnEdge), NEdgesOnCell(Mesh->NEdgesOnCell), + CellIndx("CellIndx", Mesh->NEdgesOwned, Mesh->MaxEdges2 + 2), + CellIndxSorted("CellIndxSorted", Mesh->NEdgesOwned, 2, + Mesh->MaxEdges2 + 2), + CellID(Mesh->CellID), AdvMaskHighOrder(AdvMaskHighOrder), + EdgesOnEdge(Mesh->EdgesOnEdge), CellsOnCell(Mesh->CellsOnCell), + CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge), + DvEdge(Mesh->DvEdge), AdvCoefs(AdvCoefs), AdvCoefs3rd(AdvCoefs3rd), + DerivTwo(DerivTwo) {} } // namespace OMEGA diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index 6e187a2976e4..1ceec927f4d4 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -3,6 +3,7 @@ #include "DataTypes.h" #include "HorzMesh.h" +#include "HorzUtil.h" namespace OMEGA { @@ -186,5 +187,461 @@ class InterpCellToEdge { I4 VertexDegree; }; +class SecondDerivativeOnCell { + /* + \brief deriv two computation + \author Doug Jacobsen, Bill Skamarock + \date 03/09/12 + \details + This routine precomputes the second derivative values for tracer + advection. It computes cell coefficients for the polynomial fit + as described in: + Skamarock, W. C., & Gassmann, A. (2011). + Conservative Transport Schemes for Spherical Geodesic Meshs: + High-Order Flux Operators for ODE-Based Time Integration. + Monthly Weather Review, 139(9), 2962-2975. + doi:10.1175/MWR-D-10-05056.1 + This is performed during model initialization. + */ + public: + SecondDerivativeOnCell(HorzMesh const *Mesh); + KOKKOS_FUNCTION ~SecondDerivativeOnCell() {} + KOKKOS_FUNCTION void operator()(const Array3DReal &DerivTwo, + const int ICell) const { + const int NEdges = NEdgesOnCell(ICell); + if (MaxMaxEdges < NEdges) + printf("Error: Number of edges on cell:%d exceeds maximum " + "expected:%d for cell:%d", + NEdges, MaxMaxEdges, ICell); + + // check to see if we are reaching outside the halo + auto CellList = Kokkos::subview(CellListCell, ICell, Kokkos::ALL); + bool doCell = true; + CellList[0] = ICell; + for (int J = 1; J <= NEdges; ++J) + CellList[J] = CellsOnCell(ICell, J - 1); + + for (int I = 0; I <= NEdges; ++I) + if (NCellsAll <= CellList[I]) + doCell = false; + if (!doCell) + return; + + auto XP = Kokkos::subview(XPCell, ICell, Kokkos::ALL); + auto YP = Kokkos::subview(YPCell, ICell, Kokkos::ALL); + auto Angle2D = Kokkos::subview(Angle2DCell, ICell, Kokkos::ALL); + auto B = Kokkos::subview(BCell, ICell, Kokkos::ALL, Kokkos::ALL); + if (OnSphere) { + const Array1DI4 edgesOnCell = + Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); + DetermineSphericalPatchGeometry( + NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, XVertex, + YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs[ICell]); + } else { // On an x-y plane + const Array1DI4 edgesOnCell = + Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); + DeterminePlanerPatchGeometry(ICell, NEdges, edgesOnCell, CellsOnEdge, + AngleEdge, DcEdge, XP, YP, Angle2D); + } + LeastSquaresFit(XP, YP, NEdges, B); + + // fill second derivative stencil for rk advection + for (int I = 0; I < NEdges; ++I) { + const I4 IEdge = EdgesOnCell(ICell, I); + const I4 Ind = (ICell == CellsOnEdge(IEdge, 0)) ? 0 : 1; + const Real Theta = Angle2D[I]; + const Real x = Kokkos::cos(Theta); + const Real y = Kokkos::sin(Theta); + const Real xx = x * x; + const Real xy = x * y; + const Real yy = y * y; + for (int J = 0; J <= NEdges; ++J) + DerivTwo(J, Ind, IEdge) = 2._Real * xx * B(3, J) + + 2._Real * xy * B(4, J) + + 2._Real * yy * B(5, J); + } + } + + private: + // MaxMaxEdges is used to dimention arrays that include ICell and the + // neighbor cells, so it is technically one more than MaxEdges. + static const I4 MaxMaxEdges = 10; + static constexpr R8 Pii = 3.141592653589793_Real; + + const bool OnSphere; + const I4 NCellsAll; + const I4 MaxEdges; + Array1DI4 NEdgesOnCell; + Array2DI4 EdgesOnCell; + Array2DI4 CellsOnCell; + Array2DI4 CellsOnEdge; + Array2DI4 VerticesOnEdge; + + Array1DReal XCell; + Array1DReal YCell; + Array1DReal ZCell; + Array1DReal DvEdge; + Array1DReal DcEdge; + Array1DReal AngleEdge; + Array1DReal AreaCell; + Array2DReal EdgeSignOnCell; + Array1DReal XVertex; + Array1DReal YVertex; + Array1DReal ZVertex; + + Array1DReal ThetaAbs; + Array2DReal XPCell; + Array2DReal YPCell; + Array2DReal Angle2DCell; + Array3DReal BCell; + Array2DI4 CellListCell; + + protected: + KOKKOS_INLINE_FUNCTION static void LeastSquaresFit(const Array1DReal XP, + const Array1DReal YP, + const int NEdges, + Array2DReal B) { + constexpr int NA = 6; + Real P[MaxMaxEdges][NA] = {}; + + // (polynomial_order == 2) is the only order supported + // The first row is for the origin (0,0) since the data + // is relative to the ICell centroid. + P[0][0] = 1._Real; + for (int I = 1; I <= NEdges; ++I) { + P[I][0] = 1._Real; + P[I][1] = XP[I - 1]; + P[I][2] = YP[I - 1]; + P[I][3] = XP[I - 1] * XP[I - 1]; + P[I][4] = XP[I - 1] * YP[I - 1]; + P[I][5] = YP[I - 1] * YP[I - 1]; + } + poly_fit_2(P, B, 1 + NEdges); + } + + KOKKOS_INLINE_FUNCTION static void DeterminePlanerPatchGeometry( + const int ICell, const int NEdges, const Array1DI4 EdgesOnCell, + const Array2DI4 CellsOnEdge, const Array1DReal AngleEdge, + const Array1DReal DcEdge, Array1DReal XP, Array1DReal YP, + Array1DReal Angle2D) { + for (int I = 0; I < NEdges; ++I) { + const auto IEdge = EdgesOnCell[I]; + Angle2D[I] = AngleEdge(IEdge); + if (ICell != CellsOnEdge(IEdge, 0)) + Angle2D[I] -= Pii; + XP[I] = DcEdge(IEdge) * Kokkos::cos(Angle2D[I]); + YP[I] = DcEdge(IEdge) * Kokkos::sin(Angle2D[I]); + } + } + KOKKOS_INLINE_FUNCTION static void DetermineSphericalPatchGeometry( + const int NEdges, const Array1DI4 EdgesOnCell, + const Array2DI4 VerticesOnEdge, const Array1DReal XCell, + const Array1DReal YCell, const Array1DReal ZCell, + const Array1DReal XVertex, const Array1DReal YVertex, + const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, + Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { + const Real length_scale = 1._Real; + const Real sphereRadius = distance(XCell(0), YCell(0), ZCell(0)); + Real XC[MaxMaxEdges] = {}; + Real YC[MaxMaxEdges] = {}; + Real ZC[MaxMaxEdges] = {}; + + Real Thetat_prev = 0; + Real Thetav_prev = 0; + + for (int I = 0; I <= NEdges; ++I) { + const auto J = CellList[I]; + XC[I] = XCell(J) / sphereRadius; + YC[I] = YCell(J) / sphereRadius; + ZC[I] = ZCell(J) / sphereRadius; + } + if (ZC[0] == 1.0_Real) + ThetaAbs = Pii / 2._Real; + else if (1 - ZC[0] < 1.0e-6) + ThetaAbs = Pii / 2._Real - (1 - ZC[0]) * std::atan2(YC[0], XC[0]); + else + ThetaAbs = + Pii / 2._Real - sphere_angle(XC[0], YC[0], ZC[0], XC[1], YC[1], + ZC[1], 0._Real, 0._Real, 1._Real); + + for (int I = 0; I < NEdges; ++I) { + int Ip1 = I + 1, Ip2 = I + 2; + if (NEdges < Ip2) + Ip2 = 1; + // angles from cell center to neighbor centers (thetav) + const Real Thetav = sphere_angle(XC[0], YC[0], ZC[0], XC[Ip1], YC[Ip1], + ZC[Ip1], XC[Ip2], YC[Ip2], ZC[Ip2]); + Real Dl_sphere = sphereRadius * arc_length(XC[0], YC[0], ZC[0], + XC[Ip1], YC[Ip1], ZC[Ip1]); + + Dl_sphere /= length_scale; + // Thetat = 0. this defines the x direction, + // cell center 0 -> this defines the x direction, longitude line + const Real Thetat = I ? Thetat_prev + Thetav_prev : ThetaAbs; + Thetat_prev = Thetat; + Thetav_prev = Thetav; + + XP[I] = Kokkos::cos(Thetat) * Dl_sphere; + YP[I] = Kokkos::sin(Thetat) * Dl_sphere; + + const I4 IEdge = EdgesOnCell[I]; + Real XV[2] = {}, YV[2] = {}, ZV[2] = {}, EC[3] = {}; + for (int J = 0; J < 2; ++J) { + const I4 vert = VerticesOnEdge(IEdge, J); + XV[J] = XVertex(vert) / sphereRadius; + YV[J] = YVertex(vert) / sphereRadius; + ZV[J] = ZVertex(vert) / sphereRadius; + } + arc_bisect(XV[0], YV[0], ZV[0], XV[1], YV[1], ZV[1], EC[0], EC[1], + EC[2]); + Real ThTmp = sphere_angle(XC[0], YC[0], ZC[0], XC[Ip1], YC[Ip1], + ZC[Ip1], EC[0], EC[1], EC[2]); + ThTmp += Thetat; + Angle2D[I] = ThTmp; + } + } +}; + +class MasksAndCoefficients { + + KOKKOS_INLINE_FUNCTION static void swap(Array2DI4 &vec, const int m, + const int n) { + for (int i : {0, 1}) { + const I4 j = vec(i, m); + vec(i, m) = vec(i, n); + vec(i, n) = j; + } + } + // Sort the second dimension (values) of vec based on the first (keys): + // Array1DI4 keys = Kokkos::subview(vec, 0, Kokkos::ALL); + // Array1DI4 values = Kokkos::subview(vec, 1, Kokkos::ALL); + KOKKOS_INLINE_FUNCTION static int partition(Array2DI4 &vec, const int low, + const int high) { + // Selecting last element as the pivot + const I4 pivot = vec(0, high); + int i = (low - 1); + for (int j = low; j < high; ++j) { + if (vec(0, j) <= pivot) { + i++; + swap(vec, i, j); + } + } + // Put pivot to its position + swap(vec, i + 1, high); + // Return the point of partition + return (i + 1); + } + + KOKKOS_INLINE_FUNCTION static void sort_by_key(Array2DI4 &vec, const I4 low, + const I4 high) { + // Base case: This part will be executed till the starting + // index low is lesser than the ending index high + if (low < high) { + // pi is Partitioning Index, vec[pi] is now at + // right place + const I4 pi = partition(vec, low, high); + // Separately sort elements before and after the + // Partition Index pi + sort_by_key(vec, low, pi - 1); + sort_by_key(vec, pi + 1, high); + } + } + KOKKOS_INLINE_FUNCTION static bool is_sorted(Array2DI4 &vec) { + bool sorted = true; + const I4 N = vec.extent(1) - 1; + for (I4 I = 0; I < N && sorted; ++I) + if (vec(0, I + 1) < vec(0, I)) + sorted = false; + return sorted; + } + KOKKOS_INLINE_FUNCTION static I4 search(Array1DI4 &vec, const I4 x) { + I4 low = 0, high = vec.extent(0) - 1; + while (low <= high) { + const I4 mid = low + (high - low) / 2; + if (vec(mid) == x) + return mid; + else if (vec(mid) < x) + low = mid + 1; + else + high = mid - 1; + } + // If we reach here, then element was not present + return -1; + } + KOKKOS_INLINE_FUNCTION static bool found_in_list(const Array1DI4 &List, + const I4 N, const I4 X) { + bool found = false; + for (I4 I = 0; I < N && !found; ++I) + if (X == List(I)) + found = true; + return found; + } + + public: + MasksAndCoefficients(HorzMesh const *Mesh, const Array3DReal DerivTwo, + Array1DI4 NAdvCellsForEdge, Array2DI4 AdvCellsForEdge, + Array1DI4 AdvMaskHighOrder, Array2DReal AdvCoefs, + Array2DReal AdvCoefs3rd); + + KOKKOS_FUNCTION void operator()(const int IEdge) const { + + // Array1DI4 PatchCellList = Kokkos::subview(PatchCellLists, IEdge, + // Kokkos::ALL); for (I4 I=0; I< PatchCellList.extent(0)) + // PatchCellList[I] = -1; + NAdvCellsForEdge(IEdge) = 0; + Array1DI4 CellIndex = Kokkos::subview(CellIndx, IEdge, Kokkos::ALL); + Array2DI4 CellIndexSorted = + Kokkos::subview(CellIndxSorted, IEdge, Kokkos::ALL, Kokkos::ALL); + for (unsigned I = 0; I < CellIndexSorted.extent(0); ++I) + for (unsigned K = 0; K < CellIndexSorted.extent(1); ++K) + CellIndexSorted(I, K) = MaxI4; + const int Cell1 = CellsOnEdge(IEdge, 0); + const int Cell2 = CellsOnEdge(IEdge, 1); + // at boundaries, must stay at low order + AdvMaskHighOrder(IEdge) = 1; + for (int K = 0; K < NEdgesOnCell(Cell1); ++K) + if (CellsOnCell(Cell1, K) == NCellsGlobal) + AdvMaskHighOrder(IEdge) = 0; + + for (int K = 0; K < NEdgesOnCell(Cell2); ++K) + if (CellsOnCell(Cell2, K) == NCellsGlobal) + AdvMaskHighOrder(IEdge) = 0; + // do only if this edge flux is needed to update owned cells + if (Cell1 < NCellsAll && Cell2 < NCellsAll) { + // Insert cellsOnEdge to list of advection cells + // insert_into_list(PatchCellList,Cell1); + // insert_into_list(PatchCellList,Cell2); + CellIndex(0) = Cell1; + CellIndex(1) = Cell2; + CellIndexSorted(0, 0) = CellID(Cell1); + CellIndexSorted(1, 0) = Cell1; + CellIndexSorted(0, 1) = CellID(Cell2); + CellIndexSorted(1, 1) = Cell2; + int N = 2; + // Build unique list of cells used for advection on edge + // by expanding to the extended neighbor cells + for (int I = 0; I < NEdgesOnCell(Cell1); ++I) { + const I4 CellOnCell = CellsOnCell(Cell1, I); + if (!found_in_list(CellIndex, N, CellOnCell)) { + CellIndex(N) = CellOnCell; + CellIndexSorted(0, N) = CellID(CellOnCell); + CellIndexSorted(1, N) = CellOnCell; + ++N; + } + } + for (int I = 0; I < NEdgesOnCell(Cell2); ++I) { + const I4 CellOnCell = CellsOnCell(Cell2, I); + if (!found_in_list(CellIndex, N, CellOnCell)) { + CellIndex(N) = CellOnCell; + CellIndexSorted(0, N) = CellID(CellOnCell); + CellIndexSorted(1, N) = CellOnCell; + ++N; + } + } + // sort the cell indices by cellID + sort_by_key(CellIndexSorted, 0, CellIndexSorted.extent(1) - 1); + + Array1DI4 keys = Kokkos::subview(CellIndexSorted, 0, Kokkos::ALL); + Array1DI4 values = Kokkos::subview(CellIndexSorted, 1, Kokkos::ALL); + // store local cell indices for high-order calculations + NAdvCellsForEdge(IEdge) = N; + for (int ICell = 0; ICell < N; ++ICell) + AdvCellsForEdge(IEdge, ICell) = values(ICell); + // equation 7 in Skamarock, W. C., & Gassmann, A. (2011): + // F(u,psi)_{i+1/2} = u_{i+1/2} * + // [ 1/2 (psi_{i+1} + psi_i) term 1 + // - 1/12(dx^2psi_{i+1} + dx^2psi_i) term 2 + // + sign(u) beta/12 (dx^2psi_{i+1} - dx^2psi_i)] term 3 + // (note minus sign) + // + // advCoefs accounts for terms 1 and 2 in SG11 equation 7. + // Term 1 is the 2nd-order flux-function term. advCoefs + // accounts for this with the "+ 0.5" lines below. In the + // advection routines that use these coefficients, the + // 2nd-order flux loop is then skipped. Term 2 is the + // 4th-order flux-function term. advCoefs accounts for + // term 3, the beta term. beta > 0 corresponds to the + // third-order flux function. The - sign in the derivTwo + // accumulation is for the i+1 part of term 3, while + // the + sign is for the i part. + + for (int I = 0; I < NAdvCellsMax; ++I) { + AdvCoefs(I, IEdge) = 0._Real; + AdvCoefs3rd(I, IEdge) = 0._Real; + } + // pull together third and fourth order contributions to the flux first + // from cell1 + Array1DI4 keys_edge = Kokkos::subview( + keys, Kokkos::make_pair(0, NAdvCellsForEdge(IEdge))); + if (const I4 I = search(keys_edge, CellID(Cell1)); -1 != I) { + AdvCoefs(I, IEdge) += DerivTwo(0, 0, IEdge); + AdvCoefs3rd(I, IEdge) += DerivTwo(0, 0, IEdge); + } + for (int ICell = 0; ICell < NEdgesOnCell(Cell1); ++ICell) { + if (const I4 I = + search(keys_edge, CellID(CellsOnCell(Cell1, ICell))); + -1 != I) { + AdvCoefs(I, IEdge) += DerivTwo(ICell + 1, 0, IEdge); + AdvCoefs3rd(I, IEdge) += DerivTwo(ICell + 1, 0, IEdge); + } + } + // pull together third and fourth order contributions to the flux first + // from cell2 + if (const I4 I = search(keys_edge, CellID(Cell2)); -1 != I) { + AdvCoefs(I, IEdge) += DerivTwo(0, 1, IEdge); + AdvCoefs3rd(I, IEdge) -= DerivTwo(0, 1, IEdge); + } + for (int ICell = 0; ICell < NEdgesOnCell(Cell2); ++ICell) { + if (const I4 I = + search(keys_edge, CellID(CellsOnCell(Cell2, ICell))); + -1 != I) { + AdvCoefs(I, IEdge) += DerivTwo(ICell + 1, 1, IEdge); + AdvCoefs3rd(I, IEdge) -= DerivTwo(ICell + 1, 1, IEdge); + } + } + for (int ICell = 0; ICell < NAdvCellsForEdge(IEdge); ++ICell) { + AdvCoefs(ICell, IEdge) *= -DcEdge(IEdge) * DcEdge(IEdge) / 12._Real; + AdvCoefs3rd(ICell, IEdge) *= + -DcEdge(IEdge) * DcEdge(IEdge) / 12._Real; + } + // 2nd order centered contribution place this in the main flux weights + if (const I4 I = search(keys_edge, CellID(Cell1)); -1 != I) { + AdvCoefs(I, IEdge) += 0.5_Real; + } + if (const I4 I = search(keys_edge, CellID(Cell2)); -1 != I) { + AdvCoefs(I, IEdge) += 0.5_Real; + } + // multiply by edge length - thus the flux is just dt*ru times the + // results of the vector-vector multiply + for (int ICell = 0; ICell < NAdvCellsForEdge(IEdge); ++ICell) { + AdvCoefs(ICell, IEdge) *= DvEdge(IEdge); + AdvCoefs3rd(ICell, IEdge) *= DvEdge(IEdge); + } + } + } + + private: + const I4 MaxI4 = std::numeric_limits::max(); + const I4 NCellsGlobal; + const I4 NCellsAll; + const I4 NAdvCellsMax; + Array2DI4 PatchCellLists; + Array1DI4 NAdvCellsForEdge; + Array2DI4 AdvCellsForEdge; + Array1DI4 NEdgesOnEdge; + Array1DI4 NEdgesOnCell; + Array2DI4 CellIndx; + Array3DI4 CellIndxSorted; + Array1DI4 CellID; + Array1DI4 AdvMaskHighOrder; + Array2DI4 EdgesOnEdge; + Array2DI4 CellsOnCell; + Array2DI4 CellsOnEdge; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array2DReal AdvCoefs; + Array2DReal AdvCoefs3rd; + Array3DReal DerivTwo; +}; } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/HorzUtil.h b/components/omega/src/ocn/HorzUtil.h new file mode 100644 index 000000000000..cc27927f8834 --- /dev/null +++ b/components/omega/src/ocn/HorzUtil.h @@ -0,0 +1,197 @@ +#ifndef OMEGA_HORZUTIL_H +#define OMEGA_HORZUTIL_H + +#include "DataTypes.h" +#include "Logging.h" +namespace OMEGA { + +KOKKOS_INLINE_FUNCTION Real distance(const Real x, const Real y, const Real z) { + const Real dist = Kokkos::sqrt(x * x + y * y + z * z); + return dist; +} + +KOKKOS_INLINE_FUNCTION Real sphere_angle(const Real ax, const Real ay, + const Real az, const Real bx, + const Real by, const Real bz, + const Real cx, const Real cy, + const Real cz) { + // Computes the angle between arcs AB and AC, given points A, B, and C + // Equation numbers w.r.t. + // http://mathworld.wolfram.com/SphericalTrigonometry.html + const Real one = 1._Real; + const Real neg_one = -1._Real; + const auto a = + Kokkos::acos(Kokkos::max(Kokkos::min(bx * cx + by * cy + bz * cz, one), + neg_one)); // Eqn. (3) + const auto b = + Kokkos::acos(Kokkos::max(Kokkos::min(ax * cx + ay * cy + az * cz, one), + neg_one)); // Eqn. (2) + const auto c = + Kokkos::acos(Kokkos::max(Kokkos::min(ax * bx + ay * by + az * bz, one), + neg_one)); // Eqn. (1) + const auto ABx = bx - ax; + const auto ABy = by - ay; + const auto ABz = bz - az; + + const auto ACx = cx - ax; + const auto ACy = cy - ay; + const auto ACz = cz - az; + + const auto Dx = (ABy * ACz) - (ABz * ACy); + const auto Dy = -((ABx * ACz) - (ABz * ACx)); + const auto Dz = (ABx * ACy) - (ABy * ACx); + + const auto s = 0.5_Real * (a + b + c); + const auto sin_angle = Kokkos::sqrt(Kokkos::min( + one, Kokkos::max(0._Real, + (Kokkos::sin(s - b) * Kokkos::sin(s - c)) / + (Kokkos::sin(b) * Kokkos::sin(c))))); // Eqn. (28) + Real sa = 2._Real * + Kokkos::asin(Kokkos::max(Kokkos::min(sin_angle, one), neg_one)); + if ((Dx * ax + Dy * ay + Dz * az) < 0.0) + sa *= neg_one; + return sa; +} +KOKKOS_INLINE_FUNCTION Real arc_length(const Real ax, const Real ay, + const Real az, const Real bx, + const Real by, const Real bz) { + // Returns the length of the great circle arc from A=(ax, ay, az) to + // B=(bx, by, bz). It is assumed that both A and B lie on the surface of the + // same sphere centered at the origin. + const auto cx = bx - ax; + const auto cy = by - ay; + const auto cz = bz - az; + const auto r = Kokkos::sqrt(ax * ax + ay * ay + az * az); + const auto c = Kokkos::sqrt(cx * cx + cy * cy + cz * cz); + const auto al = r * 2.0 * Kokkos::asin(c / (2.0 * r)); + return al; +} +template +static KOKKOS_FUNCTION void mat_t_mul(Real AtA[NA][NA], const Real A[][NA], + const I4 Rows) { + // Compute AtA = A^t * A for the first N rows of A + const I4 Cols = NA; + for (int I = 0; I < Cols; ++I) + for (int J = 0; J < Cols; ++J) + for (int K = 0; K < Rows; ++K) + AtA[I][J] += A[K][I] * A[K][J]; // = A^t(I,K) * A(K,j) +} +template +KOKKOS_INLINE_FUNCTION void elgs(Real A[NA][NA], I4 Indx[NA]) { + // subroutine to perform the partial-pivoting gaussian elimination. + // a(n,n) is the original matrix in the input and transformed matrix + // plus the pivoting element ratios below the diagonal in the output. + // indx(n) records the pivoting order. tao pang 2001. + // + // A is a square matrix NA by NA but only the first N rows,columns are + // pivoted. + const I4 N = NA; + for (int I = 0; I < N; ++I) + Indx[I] = I; + // find the rescaling factors, one from each row + Real C[NA] = {}; + for (int I = 0; I < N; ++I) + for (int J = 0; J < N; ++J) + C[I] = Kokkos::max(C[I], Kokkos::abs(A[I][J])); + + // search the pivoting (largest) element from each column + for (int J = 0; J < N - 1; ++J) { + Real pi1 = 0._Real; + I4 K = 0; + for (int I = J; I < N; ++I) { + const Real pi = Kokkos::abs(A[Indx[I]][J]) / C[Indx[I]]; + if (pi1 < pi) { + pi1 = pi; + K = I; + } + } + // interchange the rows via indx(n) to record pivoting order + const I4 ITmp = Indx[J]; + Indx[J] = Indx[K]; + Indx[K] = ITmp; + for (int I = J + 1; I < N; ++I) { + const Real pj = A[Indx[I]][J] / A[Indx[J]][J]; + // record pivoting ratios below the diagonal + A[Indx[I]][J] = pj; + // modify other elements accordingly + for (int K = J + 1; K < N; ++K) + A[Indx[I]][K] -= pj * A[Indx[J]][K]; + } + } +} +template +KOKKOS_INLINE_FUNCTION void migs(Real A[NA][NA], Real X[NA][NA]) { + // subroutine to invert matrix a(n,n) with the inverse stored + // in x(n,n) in the output. tao pang 2001. + // + // Matrix sizes are NA by NA but only the first N rowx & columns are + // inverted. + I4 Indx[NA] = {}; + Real B[NA][NA] = {}; + for (int I = 0; I < NA; ++I) + B[I][I] = 1._Real; + elgs(A, Indx); + for (int I = 0; I < NA - 1; ++I) + for (int J = I + 1; J < NA; ++J) + for (int K = 0; K < NA; ++K) + B[Indx[J]][K] -= A[Indx[J]][I] * B[Indx[I]][K]; + for (int I = 0; I < NA; ++I) { + X[NA - 1][I] = B[Indx[NA - 1]][I] / A[Indx[NA - 1]][NA - 1]; + for (int J = NA - 2; 0 <= J; --J) { + X[J][I] = B[Indx[J]][I]; + for (int K = J + 1; K < NA; ++K) + X[J][I] -= A[Indx[J]][K] * X[K][I]; + X[J][I] /= A[Indx[J]][J]; + } + } +} +template +KOKKOS_INLINE_FUNCTION void mat_mul_t(Array2DReal C, const Real B[NA][NA], + const Real A[][NA], const I4 Rows) { + // Compute C = B * A^t for first N columns of C + // C assumed to be initialized to 0. + // The number of rows in C is inferred from the number of column in A and B. + const I4 Cols = NA; + for (int I = 0; I < Cols; ++I) + for (int J = 0; J < Rows; ++J) + C(I, J) = 0; + for (int I = 0; I < Cols; ++I) + for (int J = 0; J < Rows; ++J) + for (int K = 0; K < Cols; ++K) + C(I, J) += B[I][K] * A[J][K]; // Transpose A +} +template +KOKKOS_INLINE_FUNCTION void poly_fit_2(const Real A[][NA], Array2DReal B, + const I4 Rows) { + Real AtA[NA][NA] = {}; + Real AtAInv[NA][NA] = {}; + mat_t_mul(AtA, A, Rows); + migs(AtA, AtAInv); + mat_mul_t(B, AtAInv, A, Rows); +} + +KOKKOS_INLINE_FUNCTION void arc_bisect(const Real ax, const Real ay, + const Real az, const Real bx, + const Real by, const Real bz, Real &cx, + Real &cy, Real &cz) { + + // Returns the point C=(cx, cy, cz) that bisects the great circle arc from + // A=(ax, ay, az) to B=(bx, by, bz). It is assumed that A and B lie on the + // surface of a sphere centered at the origin. + + cx = 0.5_Real * (ax + bx); + cy = 0.5_Real * (ay + by); + cz = 0.5_Real * (az + bz); + if (cx == 0. && cy == 0. && cz == 0.) { + printf("arc_bisect: A and B are diametrically opposite"); + } else { + const Real d = Kokkos::sqrt(cx * cx + cy * cy + cz * cz); + const Real r = Kokkos::sqrt(ax * ax + ay * ay + az * az); + cx = r * cx / d; + cy = r * cy / d; + cz = r * cz / d; + } +} + +} // namespace OMEGA +#endif diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 7cf8a7a401a8..8fdd8ad9a8a8 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -13,6 +13,7 @@ #include "Error.h" #include "Pacer.h" #include "Tracers.h" +#include namespace OMEGA { @@ -69,7 +70,8 @@ void Tendencies::init() { CustomThickTend, CustomVelTend); DefaultTendencies->readTendConfig(&TendConfig); - + if (DefaultTendencies->TracerHighOrderHorzAdv.Enabled) + DefaultTendencies->TracerHighOrderHorzAdv.init(); } // end init //------------------------------------------------------------------------------ @@ -164,11 +166,21 @@ void Tendencies::readTendConfig( CHECK_ERROR_ABORT(Err, "Tendencies: DivFactor not found in TendConfig"); } - Err += TendConfig->get("TracerHorzAdvTendencyEnable", - this->TracerHorzAdv.Enabled); - CHECK_ERROR_ABORT( - Err, "Tendencies: TracerHorzAdvTendencyEnable not found in TendConfig"); - + if (TendConfig->existsVar("TracerHorzAdvTendencyEnable")) { + I4 Order = 1; + if (TendConfig->existsVar("TracerHorzAdvTendencyOrder")) + TendConfig->get("TracerHorzAdvTendencyOrder", Order); + if (Order == 1) + this->TracerHorzAdv.Enabled = true; + else if (Order == 2) + this->TracerHighOrderHorzAdv.Enabled = true; + else { + const std::string msg = + "TracerHorzAdvTendencyOrder: Only values are 1 and 2, found " + + std::to_string(Order); + ABORT_ERROR(msg); + } + } Err += TendConfig->get("TracerDiffTendencyEnable", this->TracerDiffusion.Enabled); CHECK_ERROR_ABORT( @@ -224,6 +236,7 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), + TracerHighOrderHorzAdv(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend) { @@ -500,10 +513,13 @@ void Tendencies::computeTracerTendenciesOnly( ) { OMEGA_SCOPE(LocTracerTend, TracerTend); OMEGA_SCOPE(LocTracerHorzAdv, TracerHorzAdv); + OMEGA_SCOPE(LocTracerHighOrderHorzAdv, TracerHighOrderHorzAdv); OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); Pacer::start("Tend:computeTracerTendenciesOnly", 1); @@ -523,6 +539,8 @@ void Tendencies::computeTracerTendenciesOnly( // compute tracer horizotal advection const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; const Array3DReal &HTracersEdge = AuxState->TracerAux.HTracersEdge; + const Array2DReal &FluxLayerThickEdge = + AuxState->LayerThicknessAux.FluxLayerThickEdge; if (LocTracerHorzAdv.Enabled) { Pacer::start("Tend:tracerHorzAdv", 2); parallelForOuter( @@ -540,6 +558,34 @@ void Tendencies::computeTracerTendenciesOnly( }); Pacer::stop("Tend:tracerHorzAdv", 2); } + if (LocTracerHighOrderHorzAdv.Enabled) { + Pacer::start("Tend:TracerHighOrderHorzAdv", 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) { + LocTracerHighOrderHorzAdv(L, IEdge, KChunk, TracerArray, + FluxLayerThickEdge, + 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) { + LocTracerHighOrderHorzAdv(LocTracerTend, L, ICell, KChunk); + }); + }); + Pacer::stop("Tend:TracerHighOrderHorzAdv", 2); + } // compute tracer diffusion const Array2DReal &MeanLayerThickEdge = @@ -711,9 +757,6 @@ void Tendencies::computeAllTendencies( int VelTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { - - Pacer::start("Tend:computeAllTendencies", 1); - AuxState->computeAll(State, TracerArray, ThickTimeLevel, VelTimeLevel); computeThicknessTendenciesOnly(State, AuxState, ThickTimeLevel, VelTimeLevel, Time); @@ -721,9 +764,6 @@ void Tendencies::computeAllTendencies( Time); computeTracerTendenciesOnly(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, Time); - - Pacer::stop("Tend:computeAllTendencies", 1); - } // end all tendency compute } // end namespace OMEGA diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 851ce1fc646c..9804483b9b54 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -65,9 +65,10 @@ class Tendencies { VelocityHyperDiffOnEdge VelocityHyperDiff; WindForcingOnEdge WindForcing; BottomDragOnEdge BottomDrag; - TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; + TracerHorzAdvOnCell TracerHorzAdv; + TracerHighOrderHorzAdvOnCell TracerHighOrderHorzAdv; // Methods to compute tendency groups void computeThicknessTendencies(const OceanState *State, diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 57631a856d08..cbd9f2a0fd8e 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -12,6 +12,7 @@ #include "AuxiliaryState.h" #include "DataTypes.h" #include "HorzMesh.h" +#include "HorzOperators.h" #include "OceanState.h" #include "Tracers.h" @@ -80,6 +81,26 @@ TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +TracerHighOrderHorzAdvOnCell::TracerHighOrderHorzAdvOnCell( + const HorzMesh *Mesh, const VertCoord *VCoord) + : HorzontalMesh(Mesh), + NAdvCellsForEdge("NumberOfCellsContribToAdvectionAtEdge", + Mesh->NEdgesOwned), + AdvCellsForEdge("IndexOfCellsContributingToAdvection", Mesh->NEdgesOwned, + Mesh->MaxEdges2 + 2), + AdvMaskHighOrder("MaskForHighOrderAdvectionTerms", Mesh->NEdgesAll), + AdvCoefs("CommonAdvectionCoefficients", Mesh->MaxEdges2 + 2, + Mesh->NEdgesAll), + AdvCoefs3rd("CommonAdvectionCoeffsForHighOrder", Mesh->MaxEdges2 + 2, + Mesh->NEdgesAll), + HighOrderFlxHorz("HigherOrderHorizontalFlux", Tracers::getNumTracers(), + Mesh->NEdgesAll, VCoord->NVertLayers / VecLength), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), + DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) { + deepCopy(HighOrderFlxHorz, 0); +} + TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), @@ -100,6 +121,30 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +void TracerHighOrderHorzAdvOnCell::init() { + const HorzMesh *Mesh = this->HorzontalMesh; + const auto MaxEdges2 = Mesh->MaxEdges2; + const auto NEdgesAll = Mesh->NEdgesAll; + const auto NCellsAll = Mesh->NCellsAll; + const auto NEdgesOwned = Mesh->NEdgesOwned; + // Allocate Kokkos arrays in member data + + SecondDerivativeOnCell secondDerivativeOnCell(Mesh); + Array3DReal DerivTwo("DerivTwo", MaxEdges2 + 2, 2, NEdgesAll); + parallelFor( + {NCellsAll}, + KOKKOS_LAMBDA(int ICell) { secondDerivativeOnCell(DerivTwo, ICell); }); + // Compute masks and coefficients + Kokkos::fence(); + MasksAndCoefficients masksAndCoefficients(Mesh, DerivTwo, NAdvCellsForEdge, + AdvCellsForEdge, AdvMaskHighOrder, + AdvCoefs, AdvCoefs3rd); + Kokkos::fence(); + parallelFor( + {NEdgesOwned}, + KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); }); + Kokkos::fence(); +} } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index cc9aa841b60c..c800442131e0 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -17,8 +17,7 @@ #include "OceanState.h" #include "VertCoord.h" -#include -#include +#include // for std::copysign namespace OMEGA { @@ -26,7 +25,7 @@ namespace OMEGA { /// arrays class ThicknessFluxDivOnCell { public: - bool Enabled; + bool Enabled = false; /// constructor declaration ThicknessFluxDivOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); @@ -80,7 +79,7 @@ class ThicknessFluxDivOnCell { /// momentum equation class PotentialVortHAdvOnEdge { public: - bool Enabled; + bool Enabled = false; /// constructor declaration PotentialVortHAdvOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); @@ -131,7 +130,7 @@ class PotentialVortHAdvOnEdge { /// Gradient of kinetic energy defined on edges, for momentum equation class KEGradOnEdge { public: - bool Enabled; + bool Enabled = false; /// constructor declaration KEGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); @@ -166,7 +165,7 @@ class KEGradOnEdge { /// acceleration, for momentum equation class SSHGradOnEdge { public: - bool Enabled; + bool Enabled = false; /// constructor declaration SSHGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); @@ -201,7 +200,7 @@ class SSHGradOnEdge { /// Laplacian horizontal mixing, for momentum equation class VelocityDiffusionOnEdge { public: - bool Enabled; + bool Enabled = false; Real ViscDel2; @@ -252,7 +251,7 @@ class VelocityDiffusionOnEdge { /// Biharmonic horizontal mixing, for momentum equation class VelocityHyperDiffOnEdge { public: - bool Enabled; + bool Enabled = false; Real ViscDel4; Real DivFactor; @@ -305,7 +304,7 @@ class VelocityHyperDiffOnEdge { /// Wind forcing class WindForcingOnEdge { public: - bool Enabled; + bool Enabled = false; Real LocRhoSw; /// constructor declaration @@ -333,7 +332,7 @@ class WindForcingOnEdge { /// Bottom drag class BottomDragOnEdge { public: - bool Enabled; + bool Enabled = false; Real Coeff; /// constructor declaration @@ -369,21 +368,20 @@ class BottomDragOnEdge { // Tracer horizontal advection term class TracerHorzAdvOnCell { public: - bool Enabled; + bool Enabled = false; TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); - KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, - I4 KChunk, const Array2DReal &NormVelEdge, + KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, const I4 L, + const I4 ICell, const I4 KChunk, + const Array2DReal &NormVelEdge, const Array3DReal &HTracersOnEdge) const { const I4 KStartCell = chunkStart(KChunk, MinLayerCell(ICell)); const I4 KLenCell = chunkLength(KChunk, KStartCell, MaxLayerCell(ICell)); const I4 KEndCell = KStartCell + KLenCell - 1; - const Real InvAreaCell = 1._Real / AreaCell(ICell); - + const Real InvAreaCell = 1._Real / AreaCell(ICell); Real HAdvTmp[VecLength] = {0}; - for (int J = 0; J < NEdgesOnCell(ICell); ++J) { const I4 JEdge = EdgesOnCell(ICell, J); const I4 KStartEdge = Kokkos::max(KStartCell, MinLayerEdgeBot(JEdge)); @@ -417,10 +415,92 @@ class TracerHorzAdvOnCell { Array1DI4 MaxLayerEdgeTop; }; +// Tracer high order horizontal advection term +class TracerHighOrderHorzAdvOnCell { + public: + bool Enabled = false; + // coefficient for blending high-order terms + Real coef3rdOrder = 0.25; + + TracerHighOrderHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + void init(); + + KOKKOS_FUNCTION void operator()(const I4 L, const I4 IEdge, const I4 KChunk, + const Array3DReal &TracerCell, + const Array2DReal &FluxLayerThickEdge, + 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; + if (AdvMaskHighOrder(IEdge)) { + for (int I = 0; I < NAdvCellsForEdge(IEdge); ++I) { + const I4 ICell = AdvCellsForEdge(IEdge, I); + for (int K = KStart; K < KEnd; ++K) { + const Real normalThicknessFlux = + FluxLayerThickEdge(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); + } + } + } else { + for (int K = KStart; K < KEnd; ++K) { + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + const Real normalThicknessFlux = + FluxLayerThickEdge(IEdge, K) * NormVelEdge(IEdge, K); + const Real tracerWgt = + DvEdge(IEdge) * 0.5_Real * normalThicknessFlux; + HighOrderFlxHorz(L, IEdge, K) += + tracerWgt * + (TracerCell(L, JCell1, K) + TracerCell(L, JCell0, 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 Real InvAreaCell = 1._Real / AreaCell(ICell); + for (int K = KStart; K < KEnd; ++K) + 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; + } + } + } + + private: + const HorzMesh *HorzontalMesh; + Array1DI4 NAdvCellsForEdge; + Array2DI4 AdvCellsForEdge; + Array1DI4 AdvMaskHighOrder; + Array2DReal AdvCoefs; + Array2DReal AdvCoefs3rd; + Array3DReal HighOrderFlxHorz; + + Array1DI4 NEdgesOnCell; + Array2DI4 EdgesOnCell; + Array2DI4 CellsOnEdge; + Array2DReal EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal AreaCell; +}; + // Tracer horizontal diffusion term class TracerDiffOnCell { public: - bool Enabled; + bool Enabled = false; Real EddyDiff2; @@ -483,7 +563,7 @@ class TracerDiffOnCell { // Tracer biharmonic horizontal mixing term class TracerHyperDiffOnCell { public: - bool Enabled; + bool Enabled = false; Real EddyDiff4; diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 0c218325b100..386459b13706 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -180,12 +180,12 @@ int main(int argc, char *argv[]) { LocCells = Mesh->NCellsOwned; Err = MPI_Allreduce(&LocCells, &SumCells, 1, MPI_INT32_T, MPI_SUM, Comm); - if (SumCells == DefDecomp->NCellsGlobal) { + if (SumCells == Mesh->NCellsGlobal) { LOG_INFO("HorzMeshTest: Sum cell ID test PASS"); } else { RetVal += 1; LOG_INFO("HorzMeshTest: Sum cell ID test FAIL {} {}", SumCells, - DefDecomp->NCellsGlobal); + Mesh->NCellsGlobal); } // Test that cell centers are on sphere @@ -247,12 +247,12 @@ int main(int argc, char *argv[]) { LocEdges = Mesh->NEdgesOwned; Err = MPI_Allreduce(&LocEdges, &SumEdges, 1, MPI_INT32_T, MPI_SUM, Comm); - if (SumEdges == DefDecomp->NEdgesGlobal) { + if (SumEdges == Mesh->NEdgesGlobal) { LOG_INFO("HorzMeshTest: Sum edge ID test PASS"); } else { RetVal += 1; LOG_INFO("HorzMeshTest: Sum edge ID test FAIL {} {}", SumEdges, - DefDecomp->NEdgesGlobal); + Mesh->NEdgesGlobal); } // Test that edge coordinates are on sphere @@ -767,6 +767,25 @@ int main(int argc, char *argv[]) { RetVal += 1; LOG_INFO("HorzMeshTest: vertex halo exhange FAIL"); } + + // Test that all Cell IDs are accounted for by + // summing the list of owned values by all tasks. The result should + // be the sum of the integers from 1 to NCellsGlobal + OMEGA::I4 RefSumIDs = 0, LocSumIDs = 0, SumIDs = 0; + OMEGA::HostArray1DI4 CellIDH = OMEGA::createHostMirrorCopy(Mesh->CellID); + for (int n = 0; n < Mesh->NCellsGlobal; ++n) + RefSumIDs += n + 1; + for (int n = 0; n < Mesh->NCellsOwned; ++n) + LocSumIDs += CellIDH(n); + Err = MPI_Allreduce(&LocSumIDs, &SumIDs, 1, MPI_INT32_T, MPI_SUM, Comm); + + if (SumIDs == RefSumIDs) { + LOG_INFO("DecompTest: Sum cell ID test PASS"); + } else { + RetVal += 1; + LOG_INFO("DecompTest: Sum cell ID test FAIL {} {}", SumIDs, RefSumIDs); + } + // Finalize Omega objects HorzMesh::clear(); Dimension::clear(); diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index a9ec84331acf..eb8b4f44b624 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -422,6 +422,499 @@ int testInterpCellToEdge(Real RTol) { return Err; } +int testsecondderivativeoncellDeterminePlanerPatchGeometry(Real RTol) { + int Err = 0; + const int C[38][2] = { + {2200000, 2000000}, {2800000, 2000000}, {1300000, 2519615}, + {1900000, 2519615}, {3100000, 2519615}, {3700000, 2519615}, + {1000000, 3039230}, {2200000, 3039230}, {2800000, 3039230}, + {4000000, 3039230}, {1300000, 3558846}, {1900000, 3558846}, + {3100000, 3558846}, {3700000, 3558846}, {1000000, 4078461}, + {2200000, 4078461}, {2800000, 4078461}, {4000000, 4078461}, + {1300000, 4598076}, {1900000, 4598076}, {3100000, 4598076}, + {3700000, 4598076}, {1000000, 5117691}, {2200000, 5117691}, + {2800000, 5117691}, {4000000, 5117691}, {1300000, 5637307}, + {1900000, 5637307}, {3100000, 5637307}, {3700000, 5637307}, + {1000000, 6156922}, {2200000, 6156922}, {2800000, 6156922}, + {4000000, 6156922}, {1300000, 6676537}, {1900000, 6676537}, + {3100000, 6676537}, {3700000, 6676537}}; + const int E[12][6] = {{0, 1, 4, 8, 7, 3}, {2, 3, 7, 11, 10, 6}, + {4, 5, 9, 13, 12, 8}, {7, 8, 12, 16, 15, 11}, + {10, 11, 15, 19, 18, 14}, {12, 13, 17, 21, 20, 16}, + {15, 16, 20, 24, 23, 19}, {18, 19, 23, 27, 26, 22}, + {20, 21, 25, 29, 28, 24}, {23, 24, 28, 32, 31, 27}, + {26, 27, 31, 35, 34, 30}, {28, 29, 33, 37, 36, 32}}; + const int CtoC = 1039230; + const R8 Pii = 3.141592653589793_Real; + + std::map, int> Edges; + for (int c = 0; c < 12; ++c) { + for (int i = 0; i < 6; ++i) { + const int j = (i + 1) % 6; + const std::pair edge(std::min(E[c][i], E[c][j]), + std::max(E[c][i], E[c][j])); + if (!Edges.count(edge)) + Edges[edge] = Edges.size(); + } + } + const int NEdges = Edges.size(); + std::vector> edgesOnCell(12); + std::vector> cellsOnEdge(NEdges); + std::vector angleEdge(NEdges); + std::vector dcEdge(NEdges, -1); + for (int c = 0; c < 12; ++c) { + for (int i = 0; i < 6; ++i) { + const int j = (i + 1) % 6; + const std::pair edge(std::min(E[c][i], E[c][j]), + std::max(E[c][i], E[c][j])); + const int e = Edges.at(edge); + edgesOnCell[c][i] = e; + if (cellsOnEdge[e].empty()) { + const int v0[2] = {C[E[c][i]][0], C[E[c][i]][1]}; + const int v1[2] = {C[E[c][j]][0], C[E[c][j]][1]}; + const int v[2] = {v1[0] - v0[0], v1[1] - v0[1]}; + const double theta = std::atan2(-v[0], v[1]); + angleEdge[e] = theta; + } + cellsOnEdge[e].push_back(c); + } + } + for (unsigned e = 0; e < cellsOnEdge.size(); ++e) { + if (1 < cellsOnEdge[e].size()) + dcEdge[e] = CtoC; + } + Array2DI4 EdgesOnCell("EdgesOnCell", 12, 6); + auto EdgesOnCellH = createHostMirrorCopy(EdgesOnCell); + for (int c = 0; c < 12; ++c) + for (int i = 0; i < 6; ++i) + EdgesOnCellH(c, i) = edgesOnCell[c][i]; + OMEGA::deepCopy(EdgesOnCell, EdgesOnCellH); + + Array2DI4 CellsOnEdge("CellsOnEdge", cellsOnEdge.size(), 2); + auto CellsOnEdgeH = createHostMirrorCopy(CellsOnEdge); + for (unsigned e = 0; e < cellsOnEdge.size(); ++e) + for (unsigned i = 0; i < cellsOnEdge[e].size(); ++i) + CellsOnEdgeH(e, i) = cellsOnEdge[e][i]; + OMEGA::deepCopy(CellsOnEdge, CellsOnEdgeH); + + Array1DReal AngleEdge("AndleEdge", angleEdge.size()); + auto AngleEdgeH = createHostMirrorCopy(AngleEdge); + for (unsigned e = 0; e < angleEdge.size(); ++e) + AngleEdgeH(e) = angleEdge[e]; + OMEGA::deepCopy(AngleEdge, AngleEdgeH); + + Array1DReal DcEdge("DcEdge", dcEdge.size()); + auto DcEdgeH = createHostMirrorCopy(DcEdge); + for (unsigned e = 0; e < dcEdge.size(); ++e) + DcEdgeH(e) = dcEdge[e]; + OMEGA::deepCopy(DcEdge, DcEdgeH); + + class SecondDerivativeOnCellTest : public SecondDerivativeOnCell { + public: + virtual ~SecondDerivativeOnCellTest() {} + KOKKOS_INLINE_FUNCTION static void + Test(const int ICell, const int NEdges, const Array1DI4 EdgesOnCell, + const Array2DI4 CellsOnEdge, const Array1DReal AngleEdge, + const Array1DReal DcEdge, Array1DReal XP, Array1DReal YP, + Array1DReal Angle2D) { + SecondDerivativeOnCell::DeterminePlanerPatchGeometry( + ICell, NEdges, EdgesOnCell, CellsOnEdge, AngleEdge, DcEdge, XP, YP, + Angle2D); + } + }; + for (int ICell = 3; ICell < 7; ICell += 3) { + const I4 NEdges = 6; + const I4 MaxEdges = 7; + Array1DReal XP("XP", MaxEdges); + Array1DReal YP("YP", MaxEdges); + Array1DReal Angle2D("Array2D", MaxEdges); + const Array1DI4 edgesOnCell = + Kokkos::subview(EdgesOnCell, ICell, Kokkos::ALL); + Kokkos::parallel_for( + 1, KOKKOS_LAMBDA(int /* dummy */) { + SecondDerivativeOnCellTest::Test(ICell, NEdges, edgesOnCell, + CellsOnEdge, AngleEdge, DcEdge, + XP, YP, Angle2D); + }); + auto XPH = createHostMirrorCopy(XP); + auto YPH = createHostMirrorCopy(YP); + auto Angle2DH = createHostMirrorCopy(Angle2D); + OMEGA::deepCopy(XPH, XP); + OMEGA::deepCopy(YPH, YP); + OMEGA::deepCopy(Angle2DH, Angle2D); + for (int i = 0; i < 6; ++i) { + const Real RTol = sizeof(Real) == 4 ? 1e-3 : 1e-5; + const double theta = -Pii / 2 + (5 == i ? -Pii / 3 : i * Pii / 3); + const double x = CtoC * std::cos(theta); + const double y = CtoC * std::sin(theta); + if (!isApprox(XPH[i], x, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDeterminePlanerPatchGeometry:x", x, + XP[i]); + } + if (!isApprox(YPH[i], y, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDeterminePlanerPatchGeometry:y", y, + YP[i]); + } + if (!isApprox(Angle2DH[i], theta, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDeterminePlanerPatchGeometry:theta", + theta, Angle2D[i]); + } + } + } + if (Err == 0) + LOG_INFO("testsecondderivativeoncellDeterminePlanerPatchGeometry: " + "Successful completion"); + return Err; +} + +int testsecondderivativeoncellLeastSquaresFit(Real RTol) { + int Err = 0; + + const int CtoC = 7; + const R8 Pii = 3.141592653589793_Real; + const I4 MaxMaxEdges = 10; + const int NEdges = 6; + class SecondDerivativeOnCellTest : public SecondDerivativeOnCell { + public: + virtual ~SecondDerivativeOnCellTest() {} + KOKKOS_INLINE_FUNCTION static void Test(const Array1DReal XP, + const Array1DReal YP, + const int NEdges, Array2DReal B) { + SecondDerivativeOnCell::LeastSquaresFit(XP, YP, NEdges, B); + } + }; + Array1DReal XP("XP", NEdges); + Array1DReal YP("YP", NEdges); + auto XPH = createHostMirrorCopy(XP); + auto YPH = createHostMirrorCopy(YP); + for (int i = 0; i < NEdges; ++i) { + const double theta = -Pii / 2 + i * Pii / 3; + XPH[i] = CtoC * std::cos(theta); + YPH[i] = CtoC * std::sin(theta); + } + OMEGA::deepCopy(XP, XPH); + OMEGA::deepCopy(YP, YPH); + Array2DReal B("B", MaxMaxEdges, MaxMaxEdges); + Kokkos::parallel_for( + 1, KOKKOS_LAMBDA(int /* dummy */) { + SecondDerivativeOnCellTest::Test(XP, YP, NEdges, B); + }); + + // From a Mathematica script. + const Real c0 = 0.0412393049421; + const Real c1 = 0.0238095238095; + const Real c2 = 0.00680272108844; + const Real c3 = 0.0117826585549; + const Real M[6][7] = { + {1., 0, 0, 0, 0, 0, 0}, + {0, 0, c0, c0, 0, -c0, -c0}, + {0, -0.047619047619, -c1, c1, 0.047619047619, c1, -c1}, + {-0.0204081632653, -0.00340136054422, c2, c2, -0.00340136054422, c2, c2}, + {0, 0, -c3, c3, 0, -c3, c3}, + {-0.0204081632653, 0.0102040816327, 0, 0, 0.0102040816327, 0, 0}}; + + auto BH = createHostMirrorCopy(B); + OMEGA::deepCopy(BH, B); + for (int i = 0; i < MaxMaxEdges; ++i) { + for (int j = 0; j < MaxMaxEdges; ++j) { + const Real RTol = sizeof(Real) == 4 ? 1e-3 : 1e-5; + const Real m = (i < 6 && j < 7) ? M[i][j] : 0; + if (i < 6 && j < 7 && !m) { + if (1e-12 < std::abs(BH(i, j))) { + Err++; + LOG_ERROR("{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellLeastSquaresFit", m, + BH(i, j)); + } + } else { + if (!isApprox(BH(i, j), m, RTol)) { + Err++; + LOG_ERROR("{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellLeastSquaresFit", m, + BH(i, j)); + } + } + } + } + if (Err == 0) + LOG_INFO( + "testsecondderivativeoncellLeastSquaresFit: Successful completion"); + return Err; +} + +int testsecondderivativeoncellDetermineSphericalPatchGeometry(Real RTol) { + int Err = 0; + const int C[38][2] = { + {2200000, 2000000}, {2800000, 2000000}, {1300000, 2519615}, + {1900000, 2519615}, {3100000, 2519615}, {3700000, 2519615}, + {1000000, 3039230}, {2200000, 3039230}, {2800000, 3039230}, + {4000000, 3039230}, {1300000, 3558846}, {1900000, 3558846}, + {3100000, 3558846}, {3700000, 3558846}, {1000000, 4078461}, + {2200000, 4078461}, {2800000, 4078461}, {4000000, 4078461}, + {1300000, 4598076}, {1900000, 4598076}, {3100000, 4598076}, + {3700000, 4598076}, {1000000, 5117691}, {2200000, 5117691}, + {2800000, 5117691}, {4000000, 5117691}, {1300000, 5637307}, + {1900000, 5637307}, {3100000, 5637307}, {3700000, 5637307}, + {1000000, 6156922}, {2200000, 6156922}, {2800000, 6156922}, + {4000000, 6156922}, {1300000, 6676537}, {1900000, 6676537}, + {3100000, 6676537}, {3700000, 6676537}}; + const int E[12][6] = {{0, 1, 4, 8, 7, 3}, {2, 3, 7, 11, 10, 6}, + {4, 5, 9, 13, 12, 8}, {7, 8, 12, 16, 15, 11}, + {10, 11, 15, 19, 18, 14}, {12, 13, 17, 21, 20, 16}, + {15, 16, 20, 24, 23, 19}, {18, 19, 23, 27, 26, 22}, + {20, 21, 25, 29, 28, 24}, {23, 24, 28, 32, 31, 27}, + {26, 27, 31, 35, 34, 30}, {28, 29, 33, 37, 36, 32}}; + const int CtoC = 1039230; + const int R = 100 * CtoC; + const R8 Pii = 3.141592653589793_Real; + + // Project coordinates to sphere. + + std::map, int> Edges; + for (int c = 0; c < 12; ++c) { + for (int i = 0; i < 6; ++i) { + const int j = (i + 1) % 6; + const std::pair edge(std::min(E[c][i], E[c][j]), + std::max(E[c][i], E[c][j])); + if (!Edges.count(edge)) + Edges[edge] = Edges.size(); + } + } + const int NEdges = Edges.size(); + std::vector> edgesOnCell(12); + std::vector> cellsOnEdge(NEdges); + std::vector angleEdge(NEdges); + std::vector> verticesOnEdge(NEdges); + std::vector dcEdge(NEdges, -1); + for (int c = 0; c < 12; ++c) { + for (int i = 0; i < 6; ++i) { + const int j = (i + 1) % 6; + const std::pair edge(std::min(E[c][i], E[c][j]), + std::max(E[c][i], E[c][j])); + const int e = Edges.at(edge); + edgesOnCell[c][i] = e; + if (cellsOnEdge[e].empty()) { + const int v0[2] = {C[E[c][i]][0], C[E[c][i]][1]}; + const int v1[2] = {C[E[c][j]][0], C[E[c][j]][1]}; + const int v[2] = {v1[0] - v0[0], v1[1] - v0[1]}; + const double theta = std::atan2(-v[0], v[1]); + angleEdge[e] = theta; + verticesOnEdge[e][0] = E[c][i]; + verticesOnEdge[e][1] = E[c][j]; + } + cellsOnEdge[e].push_back(c); + } + } + for (unsigned e = 0; e < cellsOnEdge.size(); ++e) { + if (1 < cellsOnEdge[e].size()) + dcEdge[e] = CtoC; + } + + double T[2] = {}; + for (int I = 0; I < 6; ++I) + for (int J = 0; J < 2; ++J) + T[J] += C[E[3][I]][J]; + + T[0] /= 6; + T[1] /= 6; + double X[38][2] = {}; + for (int I = 0; I < 38; ++I) + for (int J = 0; J < 2; ++J) + X[I][J] = C[I][J] - T[J]; + + Array1DReal XCell("XCell", 12); + Array1DReal YCell("YCell", 12); + Array1DReal ZCell("ZCell", 12); + auto XCellH = createHostMirrorCopy(XCell); + auto YCellH = createHostMirrorCopy(YCell); + auto ZCellH = createHostMirrorCopy(ZCell); + std::vector xcellh(12, 0); + std::vector ycellh(12, 0); + std::vector zcellh(12, 0); + + for (int I = 0; I < 12; ++I) { + T[0] = T[1] = 0; + for (int J = 0; J < 6; ++J) { + for (int K = 0; K < 2; ++K) + T[K] += X[E[I][J]][K]; + } + T[0] /= 6; + T[1] /= 6; + const double d = std::sqrt(T[0] * T[0] + T[1] * T[1]); + xcellh[I] = R * std::sin(d / R) * T[0] / d; + ycellh[I] = R * std::sin(d / R) * T[1] / d; + zcellh[I] = R * std::cos(d / R); + } + for (int I = 0; I < 12; ++I) { + XCellH(I) = xcellh[I]; + YCellH(I) = ycellh[I]; + ZCellH(I) = zcellh[I]; + } + OMEGA::deepCopy(XCell, XCellH); + OMEGA::deepCopy(YCell, YCellH); + OMEGA::deepCopy(ZCell, ZCellH); + + Array1DReal XVertex("XVertex", 38); + Array1DReal YVertex("YVertex", 38); + Array1DReal ZVertex("ZVertex", 38); + auto XVertexH = createHostMirrorCopy(XVertex); + auto YVertexH = createHostMirrorCopy(YVertex); + auto ZVertexH = createHostMirrorCopy(ZVertex); + + for (int I = 0; I < 38; ++I) { + const double x = X[I][0], y = X[I][1]; + const double d = std::sqrt(x * x + y * y); + XVertexH(I) = R * std::sin(d / R) * X[I][0] / d; + YVertexH(I) = R * std::sin(d / R) * X[I][1] / d; + ZVertexH(I) = R * std::cos(d / R); + } + OMEGA::deepCopy(XVertex, XVertexH); + OMEGA::deepCopy(YVertex, YVertexH); + OMEGA::deepCopy(ZVertex, ZVertexH); + + Array2DI4 EdgesOnCell("EdgesOnCell", 12, 6); + auto EdgesOnCellH = createHostMirrorCopy(EdgesOnCell); + for (int c = 0; c < 12; ++c) + for (int i = 0; i < 6; ++i) + EdgesOnCellH(c, i) = edgesOnCell[c][i]; + OMEGA::deepCopy(EdgesOnCell, EdgesOnCellH); + + Array2DI4 VerticesOnEdge("VerticesOnEdge", NEdges, 2); + auto VerticesOnEdgeH = createHostMirrorCopy(VerticesOnEdge); + for (int c = 0; c < NEdges; ++c) + for (int i = 0; i < 2; ++i) + VerticesOnEdgeH(c, i) = verticesOnEdge[c][i]; + OMEGA::deepCopy(VerticesOnEdge, VerticesOnEdgeH); + + Array2DI4 CellsOnEdge("CellsOnEdge", cellsOnEdge.size(), 2); + auto CellsOnEdgeH = createHostMirrorCopy(CellsOnEdge); + for (unsigned e = 0; e < cellsOnEdge.size(); ++e) + for (unsigned i = 0; i < cellsOnEdge[e].size(); ++i) + CellsOnEdgeH(e, i) = cellsOnEdge[e][i]; + OMEGA::deepCopy(CellsOnEdge, CellsOnEdgeH); + + Array1DReal AngleEdge("AndleEdge", angleEdge.size()); + auto AngleEdgeH = createHostMirrorCopy(AngleEdge); + for (unsigned e = 0; e < angleEdge.size(); ++e) + AngleEdgeH(e) = angleEdge[e]; + OMEGA::deepCopy(AngleEdge, AngleEdgeH); + + Array1DReal DcEdge("DcEdge", dcEdge.size()); + auto DcEdgeH = createHostMirrorCopy(DcEdge); + for (unsigned e = 0; e < dcEdge.size(); ++e) + DcEdgeH(e) = dcEdge[e]; + OMEGA::deepCopy(DcEdge, DcEdgeH); + + class SecondDerivativeOnCellTest : public SecondDerivativeOnCell { + public: + virtual ~SecondDerivativeOnCellTest() {} + KOKKOS_INLINE_FUNCTION static void + Test(const int NEdges, const Array1DI4 EdgesOnCell, + const Array2DI4 VerticiesOnEdge, const Array1DReal XCell, + const Array1DReal YCell, const Array1DReal ZCell, + const Array1DReal XVertex, const Array1DReal YVertex, + const Array1DReal ZVertex, const Array1DI4 CellList, Array1DReal XP, + Array1DReal YP, Array1DReal Angle2D, Real &ThetaAbs) { + SecondDerivativeOnCell::DetermineSphericalPatchGeometry( + NEdges, EdgesOnCell, VerticiesOnEdge, XCell, YCell, ZCell, XVertex, + YVertex, ZVertex, CellList, XP, YP, Angle2D, ThetaAbs); + } + }; + { + const I4 NEdges = 6; + Array1DReal XP("XP", NEdges); + Array1DReal YP("YP", NEdges); + Array1DReal Angle2D("Angle2D", NEdges); + Array1DI4 CellList("CellList", NEdges + 1); + auto CellListH = createHostMirrorCopy(CellList); + { + const I4 List[NEdges + 1] = {3, 0, 2, 5, 6, 4, 1}; + for (int I = 0; I < NEdges + 1; ++I) + CellListH[I] = List[I]; + OMEGA::deepCopy(CellList, CellListH); + } + Array1DI4 edgesOnCell = Kokkos::subview(EdgesOnCell, 3, Kokkos::ALL); + Kokkos::parallel_for( + 1, KOKKOS_LAMBDA(int /* dummy */) { + Real ThetaAbs = {}; + SecondDerivativeOnCellTest::Test( + NEdges, edgesOnCell, VerticesOnEdge, XCell, YCell, ZCell, + XVertex, YVertex, ZVertex, CellList, XP, YP, Angle2D, + ThetaAbs); + }); + auto XPH = createHostMirrorCopy(XP); + auto YPH = createHostMirrorCopy(YP); + auto Angle2DH = createHostMirrorCopy(Angle2D); + for (int i = 0; i < 6; ++i) { + const Real RTol = sizeof(Real) == 4 ? 1e-3 : 1e-5; + T[0] = T[1] = 0; + for (int J = 0; J < 6; ++J) + for (int K = 0; K < 2; ++K) + T[K] += X[E[CellListH[i + 1]][J]][K]; + T[0] /= 6; + T[1] /= 6; + + const double phi = Pii / 2 + i * 2 * Pii / 6; + // Rotate the mesh by Pi. This is what is done in + // DetermineSphericalPatchGeometry + const double x = -T[0]; // same as CtoC*std::cos(phi); + const double y = -T[1]; // same as CtoC*std::sin(phi); + + if (!isApprox(1 + XPH[i], 1 + x, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDetermineSphericalPatchGeometry:x", + x, XPH[i]); + } + if (!isApprox(1 + YPH[i], 1 + y, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDetermineSphericalPatchGeometry:y", + y, YPH[i]); + } + if (!isApprox(Angle2DH[i], phi, RTol)) { + Err++; + LOG_ERROR( + "{}: FAIL, expected {}, got {}", + "testsecondderivativeoncellDetermineSphericalPatchGeometry:phi", + phi, Angle2DH[i]); + } + } + } + if (Err == 0) + LOG_INFO("testsecondderivativeoncellDetermineSphericalPatchGeometry: " + "Successful completion"); + return Err; +} + +int testsecondderivativeoncellconstructor(Real RTol) { + int Err = 0; + const auto &Mesh = HorzMesh::getDefault(); + // Compute numerical result + const auto MaxEdges2 = Mesh->MaxEdges2; + const auto NEdgesAll = Mesh->NEdgesAll; + const auto NCellsOwned = Mesh->NCellsOwned; + + Array3DReal DerivTwo("DerivTwo", MaxEdges2, 2, NEdgesAll); + SecondDerivativeOnCell DerivativeOnCell(Mesh); + parallelFor( + {NCellsOwned}, + KOKKOS_LAMBDA(int ICell) { DerivativeOnCell(DerivTwo, ICell); }); + if (Err == 0) + LOG_INFO("OperatorsTest: TestSecondDerivativeOnCell PASS"); + return Err; +} //------------------------------------------------------------------------------ // The initialization routine for Operators testing int initOperatorsTest(const std::string &MeshFile) { @@ -473,6 +966,10 @@ int operatorsTest(const std::string &MeshFile = DefaultMeshFile) { Err += testCurl(RTol); Err += testRecon(RTol); Err += testInterpCellToEdge(RTol); + Err += testsecondderivativeoncellconstructor(RTol); + Err += testsecondderivativeoncellDeterminePlanerPatchGeometry(RTol); + Err += testsecondderivativeoncellLeastSquaresFit(RTol); + Err += testsecondderivativeoncellDetermineSphericalPatchGeometry(RTol); if (Err == 0) { LOG_INFO("OperatorsTest: Successful completion"); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 3817ca584da9..4aac2a0bb1e2 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -15,6 +15,7 @@ // //===-----------------------------------------------------------------------===/ #include "TendencyTerms.h" +#include "Config.h" #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" @@ -285,12 +286,12 @@ struct TestSetupSphere { } KOKKOS_FUNCTION Real scalarC(Real Lon, Real Lat) const { - return -(Radius / 2) * std::sqrt(3 / 2 / Pi) * std::cos(Lat) * + return -(Radius / 2) * std::sqrt(3. / 2. / Pi) * std::cos(Lat) * std::cos(Lon); } KOKKOS_FUNCTION Real tracerHyperDiff(Real Lon, Real Lat) const { - return std::sqrt(3 / 2 / Pi) * std::cos(Lat) * std::cos(Lon) / Radius; + return std::sqrt(3. / 2. / Pi) * std::cos(Lat) * std::cos(Lon) / Radius; } KOKKOS_FUNCTION Real windForcingX(Real Lon, Real Lat, @@ -1013,7 +1014,7 @@ void initTendTest(const std::string &MeshFile, int NVertLayers) { initLogging(DefEnv); // Open config file - Config("Omega"); + Config::Initialize(); Config::readAll("omega.yml"); IO::init(DefComm); diff --git a/configs/Default.yml b/configs/Default.yml index 4b9d791b9efd..9e6cd313a519 100644 --- a/configs/Default.yml +++ b/configs/Default.yml @@ -27,6 +27,7 @@ Omega: ViscDel4: 1.2e11 DivFactor: 1.0 TracerHorzAdvTendencyEnable: true + TracerHorzAdvTendencyOrder: 2 TracerDiffTendencyEnable: true EddyDiff2: 10.0 TracerHyperDiffTendencyEnable: true