diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index ccdf768b1949..0a6d964f81f0 100644 --- a/components/omega/configs/Default.yml +++ b/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 diff --git a/components/omega/src/ocn/CustomTendencyTerms.cpp b/components/omega/src/ocn/CustomTendencyTerms.cpp new file mode 100644 index 000000000000..02ac9d1a16cb --- /dev/null +++ b/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 + +//=-------------------------------------------------------------------------===/ diff --git a/components/omega/src/ocn/CustomTendencyTerms.h b/components/omega/src/ocn/CustomTendencyTerms.h new file mode 100644 index 000000000000..d2be1a64c6e9 --- /dev/null +++ b/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 diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index d161702c2179..dc7c12f74ec3 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/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); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index a60894f7875b..df3f9e8f5e34 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/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);