Rebase against develop.

Author: overfelt

components/omega/configs/Default.yml

@@ -32,7 +32,8 @@ Omega:
     Density0: 1026.0
     BottomDragTendencyEnable: false
     BottomDragCoeff: 0.0
-    TracerHorzAdvTendencyEnable: true
+    TracerHorzAdvTendencyEnable: false
+    TracerHighOrderHorzAdvTendencyEnable: true
     TracerDiffTendencyEnable: true
     EddyDiff2: 10.0
     TracerHyperDiffTendencyEnable: true

components/omega/doc/userGuide/TendencyTerms.md

@@ -42,7 +42,8 @@ the currently available tendency terms:
 | VelocityHyperDiffOnEdge | VelHyperDiffTendencyEnable | enable/disable term
 | | ViscDel4 | horizontal biharmonic mixing coefficient for normal velocity
 | | DivFactor | scale factor for the divergence term
-| TracerHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term
+| TracerHorzAdvOnCell | TracerHorzAdvTendencyEnable
+| | TracerHighOrderHorzAdvTendencyEnable | enable/disable term
 | TracerDiffOnCell | TracerDiffTendencyEnable | enable/disable term
 | | EddyDiff2 | horizontal diffusion coefficient
 | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term

components/omega/src/infra/OmegaKokkos.h

@@ -69,8 +69,8 @@ template <class T> struct ArrayRank {
    static constexpr bool Is5D = T::rank == 5;
 };
 
-using ExecSpace       = MemSpace::execution_space;
-using HostExecSpace   = HostMemSpace::execution_space;
+using ExecSpace     = MemSpace::execution_space;
+using HostExecSpace = HostMemSpace::execution_space;
 using TeamPolicy      = Kokkos::TeamPolicy<ExecSpace>;
 using TeamMember      = TeamPolicy::member_type;
 using ScratchMemSpace = ExecSpace::scratch_memory_space;

components/omega/src/infra/TimeMgr.cpp

@@ -250,6 +250,7 @@ void TimeFrac::setSeconds(R8 Seconds) { // [in] floating point seconds
    I8 N; // numerator
    I8 D; // denominator
    R8 F = 0.0;
+   I4 I = 0;
    do {
       // Compute next convergent N/D
       A = (I8)R;
@@ -382,8 +383,12 @@ void TimeFrac::getHMS(I4 &Hours,           // [out] integer hours
 R8 TimeFrac::getSeconds(void) const { // \result Time in real seconds
 
    // check for divide-by-zero
-   if (Denom == 0)
-      ABORT_ERROR("TimeMgr: TimeFrac::getSeconds encountered 0 denominator.");
+   I4 Err{0};
+   if (Denom == 0) {
+      Err = 1;
+      LOG_ERROR("TimeMgr: encountered 0 denominator.");
+      return 0.0;
+   }
 
    // convert integer fractional seconds to real result
    return (R8)Whole + (R8)Numer / (R8)Denom;
@@ -397,8 +402,12 @@ R8 TimeFrac::getSeconds(void) const { // \result Time in real seconds
 R8 TimeFrac::getHours(void) const { // \result Time in real hours
 
    // check for divide-by-zero
-   if (Denom == 0)
-      ABORT_ERROR("TimeMgr: TimeFrac::getHours encountered 0 denominator.");
+   I4 Err{0};
+   if (Denom == 0) {
+      Err = 1;
+      LOG_ERROR("TimeMgr: encountered 0 denominator.");
+      return 0.0;
+   }
 
    // make local copy for manipulation
    TimeFrac TmpTime = *this;
@@ -420,8 +429,12 @@ R8 TimeFrac::getHours(void) const { // \result Time in real hours
 R8 TimeFrac::getMinutes(void) const { // \result Time in real minutes
 
    // check for divide-by-zero
-   if (Denom == 0)
-      ABORT_ERROR("TimeMgr: TimeFrac::getMinutes encountered 0 denominator.");
+   int Err{0};
+   if (Denom == 0) {
+      Err = 1;
+      LOG_ERROR("TimeMgr: encountered 0 denominator.");
+      return 0.0;
+   }
 
    // make local copy for manipulation
    TimeFrac TmpTime = *this;
@@ -1012,13 +1025,11 @@ std::unique_ptr<Calendar> Calendar::OmegaCal = nullptr;
 //-------------------------------------------------------------------------
 // Calendar::get - retrieves pointer to model calendar
 Calendar *Calendar::get() {
-
-   // Check for valid calendar
-   if (!isDefined())
-      ABORT_ERROR("TimeMgr: Attempt to retrieve undefined calendar");
-
-   return Calendar::OmegaCal.get();
-
+   if (isDefined()) {
+      return Calendar::OmegaCal.get();
+   } else {
+      LOG_CRITICAL("Attempt to retrieve undefined calendar");
+   }
 } // end retrieve calendar pointer
 
 //------------------------------------------------------------------------------
@@ -1044,15 +1055,16 @@ CalendarKind Calendar::getKind() {
 std::vector<I4> Calendar::getDaysPerMonth() {
 
    // check if calendar defined
-   if (!isDefined())
-      ABORT_ERROR("TimeMgr: Cannot retrieve DaysPerMonth"
-                  " - Calendar not initialized");
-
-   // also check for valid vector
-   if (Calendar::OmegaCal->DaysPerMonth.size() < 1)
-      ABORT_ERROR("TimeMgr: Invalid DaysPerMonth vector in Calendar");
-
-   return Calendar::OmegaCal->DaysPerMonth;
+   if (isDefined) {
+      // also check for valid vector
+      if (Calendar::OmegaCal->DaysPerMonth.size() > 0) {
+         return Calendar::OmegaCal->DaysPerMonth;
+      } else {
+         LOG_CRITICAL("Invalid DaysPerMonth vector");
+      }
+   } else {
+      LOG_CRITICAL("Cannot retrieve calendar props - Calendar not initialized");
+   }
 }
 
 I4 Calendar::getMonthsPerYear() {
@@ -1381,6 +1393,7 @@ TimeFrac Calendar::getElapsedTime(
    // initialize the basetime result, and common temps
    TimeFrac Result(0, 0, 1);
    I8 ResultWhole{0}; // whole seconds for result
+   I8 DayOfYear{0};   // day since beginning of year
    I8 JD{0};          // Julian Day used for many conversions
    I8 HourTmp{0};     // For half-day JD conversions
    I4 FebDays{0};     // For tracking leap days
@@ -1979,18 +1992,19 @@ TimeInterval::TimeInterval(
    Units       = InUnits;
    IsCalendar  = false;
    CalInterval = 0;
+   I4 Err{0};
 
    // Now set values based on input units
    I8 Length = InLength;
    switch (Units) {
    case TimeUnits::Seconds:
-      Interval.set(Length, 0, 1);
+      Err = Interval.set(Length, 0, 1);
       break;
    case TimeUnits::Minutes:
-      Interval.setHMS(0, Length, 0);
+      Err = Interval.setHMS(0, Length, 0);
       break;
    case TimeUnits::Hours:
-      Interval.setHMS(Length, 0, 0);
+      Err = Interval.setHMS(Length, 0, 0);
       break;
    case TimeUnits::Years:  // these three are all calendar
    case TimeUnits::Months: // intervals with the input length
@@ -2019,17 +2033,18 @@ TimeInterval::TimeInterval(
    Units       = InUnits;
    IsCalendar  = false;
    CalInterval = 0;
+   I4 Err{0};
 
    // Now set values based on input units
    switch (InUnits) {
    case TimeUnits::Seconds:
-      Interval.set(Length, 0, 1);
+      Err = Interval.set(Length, 0, 1);
       break;
    case TimeUnits::Minutes:
-      Interval.setHMS(0, Length, 0);
+      Err = Interval.setHMS(0, Length, 0);
       break;
    case TimeUnits::Hours:
-      Interval.setHMS(Length, 0, 0);
+      Err = Interval.setHMS(Length, 0, 0);
       break;
    case TimeUnits::Years:  // these three are all calendar
    case TimeUnits::Months: // intervals with the input length
@@ -2060,17 +2075,18 @@ TimeInterval::TimeInterval(
    Units       = InUnits;
    IsCalendar  = false;
    CalInterval = 0;
+   I4 Err{0};
 
    // Now set values based on input units
    switch (InUnits) {
    case TimeUnits::Seconds:
-      Interval.setSeconds(Length); // TimeFrac set
+      Err = Interval.setSeconds(Length); // TimeFrac set
       break;
    case TimeUnits::Minutes:
-      Interval.setMinutes(Length); // TimeFrac set
+      Err = Interval.setMinutes(Length); // TimeFrac set
       break;
    case TimeUnits::Hours:
-      Interval.setHours(Length); // TimeFrac set
+      Err = Interval.setHours(Length); // TimeFrac set
       break;
    case TimeUnits::Years:  // these three are all calendar
    case TimeUnits::Months: // intervals with the input length
@@ -3030,10 +3046,10 @@ bool TimeInterval::isPositive(void) {
    // make sure fraction is in proper form, particularly that the
    // whole and fraction have same sign and numerator and denomintator
    // have the same sign.
-   Interval.simplify();
+   I4 Err = Interval.simplify();
 
    // now retrieve fractional seconds and check components
-   Interval.get(W, N, D);
+   Err = Interval.get(W, N, D);
    if (W >= 0) { // whole part of seconds non-negative
       // only need check the numerator with simplified form
       if (N > 0)
@@ -3513,10 +3529,9 @@ std::string TimeInstant::getString(
    const char *FormatStr = Tmp.c_str();
 
    // now use sprintf to print to C string
-   I4 NChars =
-       sprintf(TimeStr, FormatStr, Year, Month, Day, Hour, Minute, Second);
-   if (NChars < 0)
-      ABORT_ERROR("TimeMgr: TimeInstant::getString error in sprintf");
+   Err = sprintf(TimeStr, FormatStr, Year, Month, Day, Hour, Minute, Second);
+   if (Err < 0)
+      LOG_ERROR("TimeMgr: TimeInstant::getString error in sprintf");
 
    // now convert C string to std::string for result
    std::string Result(TimeStr);

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 {

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 {
@@ -57,8 +56,9 @@ void HorzMesh::init() {
 HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
                    Decomp *MeshDecomp,      //< [in] Decomp for the new mesh
                    I4 InNVertLevels         //< [in} num vertical levels
-) {
-
+                   )
+    : CellID(MeshDecomp->CellID) ///< global cell ID for each local cell
+{
    MeshName = Name;
 
    // Retrieve mesh files name from Decomp
@@ -68,11 +68,12 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
    NVertLevels = InNVertLevels;
 
    // 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;
@@ -81,6 +82,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;

components/omega/src/ocn/HorzMesh.h

@@ -103,6 +103,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
@@ -111,7 +112,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
@@ -152,6 +154,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)

components/omega/src/ocn/HorzOperators.cpp

@@ -27,4 +27,45 @@ 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,
+    Array2DI4 AdvMaskHighOrder, Array2DReal AdvCoefs, Array2DReal AdvCoefs3rd)
+    : NCellsGlobal(Mesh->NCellsGlobal), NCellsAll(Mesh->NCellsAll),
+      NAdvCellsMax(Mesh->MaxEdges2), OmegaHash(Mesh->NEdgesAll),
+      NAdvCellsForEdge(NAdvCellsForEdge), AdvCellsForEdge(AdvCellsForEdge),
+      NEdgesOnEdge(Mesh->NEdgesOnEdge), NEdgesOnCell(Mesh->NEdgesOnCell),
+      CellIndx("CellIndx", Mesh->MaxEdges2 + 2),
+      CellIndxSorted("CellIndxSorted", 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