Merge pull request #170 from hyungyukang/omega/add-manufactured-tende…

…ncies

Add manufactured solution tendency terms to Omega

components/omega/configs/Default.yml

@@ -30,6 +30,8 @@ Omega:
     EddyDiff2: 10.0
     TracerHyperDiffTendencyEnable: true
     EddyDiff4: 0.0
+    UseCustomTendency: false
+    ManufacturedSolutionTendency: false
   Tracers:
     Base: [Temperature, Salinity]
     Debug: [Debug1, Debug2, Debug3]
@@ -102,3 +104,7 @@ Omega:
       EndTime: 0001-06-30_00:00:00
       Contents:
         - Tracers
+  ManufacturedSolution:
+    WavelengthX: 5.0e6
+    WavelengthY: 4.33013e6
+    Amplitude: 1.0

components/omega/src/ocn/CustomTendencyTerms.cpp

@@ -0,0 +1,228 @@
+//===-- ocn/CustomTendencyTerms.cpp - Custom tendency terms -----*- C++ -*-===//
+//
+// The customized tendency terms can be added to the tendency terms based
+// based on an option 'UseCustomTendency' in Tendencies Config group.
+// This file contains functions for initializing customized tendency terms.
+//
+//===----------------------------------------------------------------------===//
+
+#include "CustomTendencyTerms.h"
+#include "Config.h"
+#include "TimeStepper.h"
+
+namespace OMEGA {
+
+//===-----------------------------------------------------------------------===/
+// Initialize the manufactured solution tendency terms.
+//===-----------------------------------------------------------------------===/
+int ManufacturedSolution::init() {
+   int Err;
+
+   // Get ManufacturedSolConfig group
+   Config *OmegaConfig = Config::getOmegaConfig();
+   Config ManufacturedSolConfig("ManufacturedSolution");
+   Err = OmegaConfig->get(ManufacturedSolConfig);
+   if (Err != 0) {
+      LOG_CRITICAL("ManufacturedSolution:: ManufacturedSolution group "
+                   "not found in Config");
+      return Err;
+   }
+
+   // Get TendConfig group
+   Config TendConfig("Tendencies");
+   Err = OmegaConfig->get(TendConfig);
+   if (Err != 0) {
+      LOG_CRITICAL("ManufacturedSolution:: Tendencies group "
+                   "not found in Config");
+      return Err;
+   }
+
+   // Get manufactured solution parameters from Config
+   R8 WavelengthX;
+   R8 WavelengthY;
+   R8 Amplitude;
+
+   Err = ManufacturedSolConfig.get("WavelengthX", WavelengthX);
+   if (Err != 0) {
+      LOG_ERROR("ManufacturedSolution:: WavelengthX not found in "
+                "ManufacturedSolConfig");
+      return Err;
+   }
+
+   Err = ManufacturedSolConfig.get("WavelengthY", WavelengthY);
+   if (Err != 0) {
+      LOG_ERROR("ManufacturedSolution:: WavelengthY not found in "
+                "ManufacturedSolConfig");
+      return Err;
+   }
+
+   Err = ManufacturedSolConfig.get("Amplitude", Amplitude);
+   if (Err != 0) {
+      LOG_ERROR("ManufacturedSolution:: Amplitude not found in "
+                "ManufacturedSolConfig");
+      return Err;
+   }
+
+   // Get Tendendices parameters for del2 and del4 source terms
+   Err = TendConfig.get("VelDiffTendencyEnable",
+                        ManufacturedVelTend.VelDiffTendencyEnable);
+   Err += TendConfig.get("VelHyperDiffTendencyEnable",
+                         ManufacturedVelTend.VelHyperDiffTendencyEnable);
+   Err += TendConfig.get("ViscDel2", ManufacturedVelTend.ViscDel2);
+   Err += TendConfig.get("ViscDel4", ManufacturedVelTend.ViscDel4);
+
+   if (Err != 0) {
+      LOG_ERROR("ManufacturedSolution::Error reading Tendencies config");
+      return Err;
+   }
+
+   // Get the reference time to compute the model elapsed time
+   /// Get model clock from time stepper
+   TimeStepper *DefStepper             = TimeStepper::getDefault();
+   Clock *ModelClock                   = DefStepper->getClock();
+   ManufacturedThickTend.ReferenceTime = ModelClock->getCurrentTime();
+   ManufacturedVelTend.ReferenceTime   = ManufacturedThickTend.ReferenceTime;
+
+   // Get BottomDepth for the resting thickness
+   /// This test case assumes that the restingThickness is horizontally uniform
+   /// and that only one vertical level is used so only one set of indices is
+   /// used here.
+   HorzMesh *DefHorzMesh = HorzMesh::getDefault();
+   R8 H0                 = DefHorzMesh->BottomDepthH(0);
+
+   // Define and compute common constants
+   R8 Grav    = 9.80665_Real;                          // Gravity acceleration
+   R8 Pii     = 3.141592653589793_Real;                // Pi
+   R8 Kx      = 2.0_Real * Pii / WavelengthX;          // Wave in X-dir
+   R8 Ky      = 2.0_Real * Pii / WavelengthY;          // Wave in Y-dir
+   R8 AngFreq = sqrt(H0 * Grav * (Kx * Kx + Ky * Ky)); // Angular frequency
+
+   // Assign constants for thickness tendency function
+   ManufacturedThickTend.H0      = H0;
+   ManufacturedThickTend.Eta0    = Amplitude;
+   ManufacturedThickTend.Kx      = Kx;
+   ManufacturedThickTend.Ky      = Ky;
+   ManufacturedThickTend.AngFreq = AngFreq;
+
+   // Assign constants for velocity tendency function
+   ManufacturedVelTend.Grav    = Grav;
+   ManufacturedVelTend.Eta0    = Amplitude;
+   ManufacturedVelTend.Kx      = Kx;
+   ManufacturedVelTend.Ky      = Ky;
+   ManufacturedVelTend.AngFreq = AngFreq;
+
+   return Err;
+
+} // end ManufacturedSolution init
+
+//===--------------------------------------------------------------------===/
+// Manufactured tendency term for the thickness equation
+//===--------------------------------------------------------------------===/
+void ManufacturedSolution::ManufacturedThicknessTendency::operator()(
+    Array2DReal ThicknessTend, const OceanState *State,
+    const AuxiliaryState *AuxState, int ThickTimeLevel, int VelTimeLevel,
+    TimeInstant Time) const {
+
+   // Get elapsed time since reference time
+   R8 ElapsedTimeSec;
+   TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
+   ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);
+
+   auto *Mesh       = HorzMesh::getDefault();
+   auto NVertLevels = ThicknessTend.extent_int(1);
+
+   Array1DReal XCell = Mesh->XCell;
+   Array1DReal YCell = Mesh->YCell;
+
+   OMEGA_SCOPE(LocH0, H0);
+   OMEGA_SCOPE(LocEta0, Eta0);
+   OMEGA_SCOPE(LocKx, Kx);
+   OMEGA_SCOPE(LocKy, Ky);
+   OMEGA_SCOPE(LocAngFreq, AngFreq);
+
+   parallelFor(
+       {Mesh->NCellsAll, NVertLevels}, KOKKOS_LAMBDA(int ICell, int KLevel) {
+          R8 X     = XCell(ICell);
+          R8 Y     = YCell(ICell);
+          R8 Phase = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
+          ThicknessTend(ICell, KLevel) +=
+              LocEta0 *
+              (-LocH0 * (LocKx + LocKy) * sin(Phase) - LocAngFreq * cos(Phase) +
+               LocEta0 * (LocKx + LocKy) * cos(2.0_Real * Phase));
+       });
+
+} // end void ManufacturedThicknessTendency
+
+//===--------------------------------------------------------------------===/
+// Manufactured tendency term for the momentum equation
+//===--------------------------------------------------------------------===/
+void ManufacturedSolution::ManufacturedVelocityTendency::operator()(
+    Array2DReal NormalVelTend, const OceanState *State,
+    const AuxiliaryState *AuxState, int ThickTimeLevel, int VelTimeLevel,
+    TimeInstant Time) const {
+
+   // Get elapsed time since reference time
+   R8 ElapsedTimeSec;
+   TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
+   ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);
+
+   auto *Mesh       = HorzMesh::getDefault();
+   auto NVertLevels = NormalVelTend.extent_int(1);
+
+   Array1DReal FEdge     = Mesh->FEdge;
+   Array1DReal XEdge     = Mesh->XEdge;
+   Array1DReal YEdge     = Mesh->YEdge;
+   Array1DReal AngleEdge = Mesh->AngleEdge;
+
+   OMEGA_SCOPE(LocGrav, Grav);
+   OMEGA_SCOPE(LocEta0, Eta0);
+   OMEGA_SCOPE(LocKx, Kx);
+   OMEGA_SCOPE(LocKy, Ky);
+   OMEGA_SCOPE(LocAngFreq, AngFreq);
+   OMEGA_SCOPE(LocViscDel2, ViscDel2);
+   OMEGA_SCOPE(LocViscDel4, ViscDel4);
+   OMEGA_SCOPE(LocVelDiffTendencyEnable, VelDiffTendencyEnable);
+   OMEGA_SCOPE(LocVelHyperDiffTendencyEnable, VelHyperDiffTendencyEnable);
+
+   R8 LocKx2 = LocKx * LocKx;
+   R8 LocKy2 = LocKy * LocKy;
+   R8 LocKx4 = LocKx2 * LocKx2;
+   R8 LocKy4 = LocKy2 * LocKy2;
+
+   parallelFor(
+       {Mesh->NEdgesAll, NVertLevels}, KOKKOS_LAMBDA(int IEdge, int KLevel) {
+          R8 X = XEdge(IEdge);
+          R8 Y = YEdge(IEdge);
+
+          R8 Phase       = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
+          R8 SourceTerm0 = LocAngFreq * sin(Phase) - 0.5_Real * LocEta0 *
+                                                         (LocKx + LocKy) *
+                                                         sin(2.0_Real * Phase);
+
+          R8 U = LocEta0 *
+                 ((-FEdge(IEdge) + LocGrav * LocKx) * cos(Phase) + SourceTerm0);
+          R8 V = LocEta0 *
+                 ((FEdge(IEdge) + LocGrav * LocKy) * cos(Phase) + SourceTerm0);
+
+          // Del2 and del4 source terms
+          if (LocVelDiffTendencyEnable) {
+             U += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
+             V += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
+          }
+          if (LocVelHyperDiffTendencyEnable) {
+             U -= LocViscDel4 * LocEta0 *
+                  ((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
+             V -= LocViscDel4 * LocEta0 *
+                  ((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
+          }
+
+          R8 NormalCompSourceTerm =
+              cos(AngleEdge(IEdge)) * U + sin(AngleEdge(IEdge)) * V;
+          NormalVelTend(IEdge, KLevel) += NormalCompSourceTerm;
+       });
+
+} // end void ManufacturedVelocityTendency
+
+} // end namespace OMEGA
+
+//=-------------------------------------------------------------------------===/

components/omega/src/ocn/CustomTendencyTerms.h

@@ -0,0 +1,81 @@
+#ifndef OMEGA_CUSTOMTENDENCYTERMS_H
+#define OMEGA_CUSTOMTENDENCYTERMS_H
+//===-- ocn/CustomTendencyTerms.h - Custom tendency terms -------*- C++ -*-===//
+//
+/// \file
+/// \brief Contains customized tendency terms for the thickness and momentum
+///        equations
+///
+/// For details on the manufactured solution class, see Bishnu et al. (2024)
+/// (https://doi.org/10.1029/2022MS003545) and the manufactured solution test
+/// case in Polaris. The Polaris package leverages this feature to validate
+/// an expected order of convergence of Omega.
+//
+//===----------------------------------------------------------------------===//
+
+#include "HorzMesh.h"
+#include "TendencyTerms.h"
+#include "TimeMgr.h"
+
+namespace OMEGA {
+
+//===-----------------------------------------------------------------------===/
+// A class for the manufactured solution tendency terms
+//===-----------------------------------------------------------------------===/
+class ManufacturedSolution {
+
+ public:
+   //===--------------------------------------------------------------------===/
+   // Manufactured tendency term for the thickness equation
+   //===--------------------------------------------------------------------===/
+   struct ManufacturedThicknessTendency {
+
+      // Constants defined in 'init'
+      TimeInstant ReferenceTime;
+      R8 H0;
+      R8 Eta0;
+      R8 Kx;
+      R8 Ky;
+      R8 AngFreq;
+
+      void operator()(Array2DReal ThicknessTend, const OceanState *State,
+                      const AuxiliaryState *AuxState, int ThickTimeLevel,
+                      int VelTimeLevel, TimeInstant Time) const;
+   }; // end struct ManufacturedThicknessTendency
+
+   //===--------------------------------------------------------------------===/
+   // Manufactured tendency term for the momentum equation
+   //===--------------------------------------------------------------------===/
+   struct ManufacturedVelocityTendency {
+
+      // Constants defined in 'init'
+      TimeInstant ReferenceTime;
+      R8 Grav;
+      R8 Eta0;
+      R8 Kx;
+      R8 Ky;
+      R8 AngFreq;
+      R8 ViscDel2;
+      R8 ViscDel4;
+      bool VelDiffTendencyEnable;
+      bool VelHyperDiffTendencyEnable;
+
+      void operator()(Array2DReal NormalVelTend, const OceanState *State,
+                      const AuxiliaryState *AuxState, int ThickTimeLevel,
+                      int VelTimeLevel, TimeInstant Time) const;
+
+   }; // end struct ManufacturedVelocityTendency
+
+   // Instances of manufactured tendencies
+   ManufacturedThicknessTendency ManufacturedThickTend;
+   ManufacturedVelocityTendency ManufacturedVelTend;
+
+   int init();
+
+}; // end class ManufacturedSolution
+
+} // end namespace OMEGA
+
+//===-----------------------------------------------------------------------===/
+
+#endif

components/omega/src/ocn/HorzMesh.cpp

@@ -690,6 +690,7 @@ void HorzMesh::copyToDevice() {
    WeightsOnEdge     = createDeviceMirrorCopy(WeightsOnEdgeH);
    FVertex           = createDeviceMirrorCopy(FVertexH);
    BottomDepth       = createDeviceMirrorCopy(BottomDepthH);
+   FEdge             = createDeviceMirrorCopy(FEdgeH);
    XCell             = createDeviceMirrorCopy(XCellH);
    YCell             = createDeviceMirrorCopy(YCellH);
    XEdge             = createDeviceMirrorCopy(XEdgeH);

components/omega/src/ocn/Tendencies.cpp

@@ -9,6 +9,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "Tendencies.h"
+#include "CustomTendencyTerms.h"
 #include "Tracers.h"
 
 namespace OMEGA {
@@ -35,8 +36,51 @@ int Tendencies::init() {
       return Err;
    }
 
+   // Check if use the customized tendencies
+   bool UseCustomTendency = false;
+   Err = TendConfig.get("UseCustomTendency", UseCustomTendency);
+   if (Err != 0) {
+      LOG_ERROR("Tendencies:: UseCustomTendency not found in Config");
+      return Err;
+   }
+
+   /// Instances of custom tendencies - empty by default
+   CustomTendencyType CustomThickTend;
+   CustomTendencyType CustomVelTend;
+
+   if (UseCustomTendency) {
+      // Check if use manufactured tendency terms
+      bool ManufacturedTend = false;
+      I4 ManufacturedTendErr =
+          TendConfig.get("ManufacturedSolutionTendency", ManufacturedTend);
+
+      if (ManufacturedTendErr != 0 && ManufacturedTend) {
+         LOG_CRITICAL("Tendencies: ManufacturedSolutionTendency "
+                      "not found in TendConfig");
+         return ManufacturedTendErr;
+      }
+
+      if (ManufacturedTend) {
+         ManufacturedSolution ManufacturedSol;
+         I4 ManufacturedInitErr = ManufacturedSol.init();
+
+         if (ManufacturedInitErr != 0) {
+            LOG_CRITICAL("Error in initializing the manufactured solution "
+                         "tendency terms");
+            return ManufacturedInitErr;
+         }
+
+         CustomThickTend = ManufacturedSol.ManufacturedThickTend;
+         CustomVelTend   = ManufacturedSol.ManufacturedVelTend;
+
+      } // if ManufacturedTend
+
+   } // end if UseCustomTendency
+
+   // Ceate default tendencies
    Tendencies::DefaultTendencies =
-       create("Default", DefHorzMesh, NVertLevels, NTracers, &TendConfig);
+       create("Default", DefHorzMesh, NVertLevels, NTracers, &TendConfig,
+              CustomThickTend, CustomVelTend);
 
    Err = DefaultTendencies->readTendConfig(&TendConfig);