diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index eb393b4d0982..549a6fd3bec3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -31,6 +31,9 @@ Omega: VerticalTracerFluxOrder: 3 WindStress: InterpType: Isotropic + SurfaceRestoring: + TracersToRestore: [Temperature, Salinity] + PistonVelocity: 1.585e-5 VertCoord: MovementWeightType: Uniform PressureGrad: @@ -46,6 +49,7 @@ Omega: ViscDel4: 1.2e11 DivFactor: 1.0 WindForcingTendencyEnable: false + SurfaceTracerRestoringEnable: false BottomDragTendencyEnable: false BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true diff --git a/components/omega/doc/devGuide/AuxiliaryVariables.md b/components/omega/doc/devGuide/AuxiliaryVariables.md index d75a8f2b3f84..bd43757a42e3 100644 --- a/components/omega/doc/devGuide/AuxiliaryVariables.md +++ b/components/omega/doc/devGuide/AuxiliaryVariables.md @@ -104,3 +104,9 @@ The following auxiliary variable groups are currently implemented: | WindForcingAuxVars | ZonalStressCell || || MeridStressCell || || NormalStressEdge || +| SurfTracerRestAuxVars | TracersMonthlySurfClimoCell || + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-dev-forcing). diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md new file mode 100644 index 000000000000..5d2d52f20d6b --- /dev/null +++ b/components/omega/doc/devGuide/Forcing.md @@ -0,0 +1,79 @@ +(omega-dev-forcing)= + +# Forcing + +This page describes design and implementation details for forcing-related +pathways in Omega, currently this includes: + +- Wind forcing +- Surface tracer restoring + +## Wind forcing design + +### Wind forcing data flow + +1. External fields provide: + - `WindStressZonal` + - `WindStressMeridional` +2. Auxiliary-state compute builds `NormalStressEdge` from those fields. +3. Tendency term applies wind-stress forcing to edge-normal velocity tendency. + +### Wind forcing key classes/components + +- `WindForcingAuxVars` + - Stores wind-stress cell fields and computed `NormalStressEdge` + - Applies configured interpolation choice (`InterpType`) +- `AuxiliaryState::computeMomAux` + - Calls `WindForcingAuxVars::computeVarsOnEdge` +- `WindForcingOnEdge` tendency term + - Adds contribution proportional to normal stress and inverse layer + thickness in the surface layer + +### Wind forcing config coupling + +- `Omega.WindStress.InterpType` + - mapped to `InterpCellToEdgeOption` +- `Omega.Tendencies.WindForcingTendencyEnable` + - gates execution of wind forcing tendency kernel + +## Surface tracer restoring design + +### Surface tracer restoring data flow + +1. External fields provide target values: `TracersMonthlySurfClimoCell` (values and units should match the state variables) +2. Auxiliary-state compute forms restoring differences: `SurfTracerRestoringDiffsCell = target - tracer_surface` +3. Tendency term applies restoring only at surface layer and only for tracers selected from `SurfaceRestoring.TracersToRestore`. + +### Surface tracer restoring key classes/components + +- `SurfTracerRestAuxVars` + - Inputs: `TracersMonthlySurfClimoCell`, tracer state array + - Output: `SurfTracerRestoringDiffsCell` + - Uses `MinLayerCell` to select surface layer index +- `SurfaceTracerRestoringOnCell` tendency term + - Applies `PistonVelocity * SurfTracerRestoringDiffsCell` at surface +- `Tendencies` + - Parses `SurfaceRestoring.TracersToRestore` and resolves tracer indices + - Builds `TracerIdsToRestore` and `NTracersToRestore` + - Applies tracer-selection logic at call site in + `computeTracerTendenciesOnly` + - Aborts if restoring is enabled but no tracer IDs are available + +### Surface tracer restoring config coupling + +- `Omega.SurfaceRestoring.PistonVelocity` + - tendency scaling +- `Omega.SurfaceRestoring.TracersToRestore` + - tracer-level enable list used to build `TracerIdsToRestore` +- `Omega.Tendencies.SurfaceTracerRestoringEnable` + - gates restoring tendency execution + +## Notes + +- If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. +- If restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error. +- It is assumed that the incoming `TracersMonthlySurfClimoCell` fields (values and units) match the Omega state variables (i.e. conservative temperature and absolute salinity for Teos-10). If not, a pre-processing conversion should be implemented. +- Surface tracer restoring is active everywhere if enabled. A flag to turn it off under sea ice will need to be added in later development if this feature is desired. +- Unlike MPAS-Ocean, a `MaxDiff` clamping is not applied here. This check should instead be implemented in Ocean Validate when that is available. +- A global scaling to ensure zero-sum has not been implemented for the surface tracer restoring, but should be added in later development. +- At this stage, there is no temporal interpolation applied to the restoring targets, the raw `TracersMonthlySurfClimoCell` snapshot is used. diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index e665fb03ed76..08005b1c00be 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -41,3 +41,9 @@ implemented: - `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` +- `SurfaceTracerRestoringOnCell` + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-dev-forcing). diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 220965ff1fdb..6712bb3e7f8d 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -53,6 +53,7 @@ userGuide/PGrad userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv +userGuide/Forcing ``` ```{toctree} @@ -98,6 +99,7 @@ devGuide/PGrad devGuide/Timing devGuide/VerticalMixingCoeff devGuide/VertAdv +devGuide/Forcing ``` ```{toctree} diff --git a/components/omega/doc/userGuide/AuxiliaryVariables.md b/components/omega/doc/userGuide/AuxiliaryVariables.md index 4f6fbc18c616..f6e377c388f9 100644 --- a/components/omega/doc/userGuide/AuxiliaryVariables.md +++ b/components/omega/doc/userGuide/AuxiliaryVariables.md @@ -32,3 +32,10 @@ The following auxiliary variables are currently available: | ZonalStressCell | zonal component of wind stress on cells | MeridStressCell | meridional component of wind stress on cells | NormalStressEdge | normal component of wind stress on edge +| SurfTracerRestoringDiffsCell | surface tracer restoring differences on cells +| TracersMonthlySurfClimoCell | monthly climatology values to restore to for surface tracer on cells + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-user-forcing). diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md new file mode 100644 index 000000000000..469693958b41 --- /dev/null +++ b/components/omega/doc/userGuide/Forcing.md @@ -0,0 +1,86 @@ +(omega-user-forcing)= + +# Forcing + +This page documents the user-facing configuration and behavior for current forcing in Omega: + +- Wind forcing +- Surface tracer restoring + +## Wind forcing + +Wind forcing adds momentum tendency from surface wind stress. + +### Wind forcing configuration + +Wind forcing behavior is controlled by two configuration blocks: + +```yaml +Omega: + WindStress: + InterpType: Isotropic + + Tendencies: + WindForcingTendencyEnable: true +``` + +- `WindStress.InterpType` + - `Isotropic`: isotropic cell-to-edge interpolation for wind stress + - `Anisotropic`: anisotropic interpolation option +- `Tendencies.WindForcingTendencyEnable`: switch to enable wind forcing tendency + +### Required input fields + +Wind forcing uses auxiliary wind-stress fields: + +- `WindStressZonal` +- `WindStressMeridional` + +These are used to form edge-normal stress (`NormalStressEdge`) that enters +momentum tendencies. + +## Surface tracer restoring + +Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean +surface for selected tracers. This is implemented to mitigate drifts in chosen tracers +(most often salinity) by nudging the model's simulated tracer values towards observed climatological values. +This process prevents oceanic regimes from shifting away from reality due to errors in surface freshwater +forcing (in the case of salinity restoring). Currently, it is applied everywhere when enabled. + +### Surface tracer restoring configuration + +Surface tracer restoring is controlled by two configuration blocks: + +```yaml +Omega: + SurfaceRestoring: + TracersToRestore: [Temperature, Salinity] + PistonVelocity: 1.585e-5 + + Tendencies: + SurfaceTracerRestoringEnable: true +``` + +- `TracersToRestore`: list of tracer names that restoring is applied to +- `PistonVelocity`: restoring rate coefficient +- `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring + +When restoring is enabled, Omega resolves `TracersToRestore` into an internal +list of tracer IDs and applies restoring only to tracers in that list. + +### Restoring target fields + +Surface restoring uses auxiliary fields: + +- `TracersMonthlySurfClimoCell`: restoring target climatological values +- `SurfTracerRestoringDiffsCell`: computed target-minus-state differences + +The restoring tendency is computed at the surface layer only and is limited by +the configured `PistonVelocity` and target-minus-state difference. + +## Notes + +- If a tracer is not listed in `TracersToRestore`, no restoring tendency is + applied to that tracer. +- If surface restoring is enabled but no tracer IDs are available at tendency + compute-time, Omega aborts with an error. diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 5757e98ebe0f..4dbcfcb13519 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -20,6 +20,7 @@ tendency terms are currently implemented: | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers | WindForcingOnEdge | forcing by wind stress, defined on edges | BottomDragOnEdge | bottom drag, defined on edges +| SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells Among the internal data stored by each functor is a `bool` which can enable or disable the contribution of that particular term to the tendency. These flags @@ -54,7 +55,7 @@ the currently available tendency terms: | WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term | | BottomDragCoeff | bottom drag coefficient - +| SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term ## Second Order Horizontal Advection Algorithm @@ -136,3 +137,8 @@ is about order 1.7 as shown in {numref}`tracer-higher-order-convergence`: :width: 600 px Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate ``` + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-user-forcing). diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 4de2a47884f0..30efe9b5ffc5 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -3,6 +3,7 @@ #include "Field.h" #include "Logging.h" #include "Pacer.h" +#include "Tendencies.h" #include "TimeStepper.h" namespace OMEGA { @@ -28,6 +29,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), WindForcingAux(stripDefault(Name), Mesh), + SurfTracerRestAux(stripDefault(Name), Mesh, NTracers), TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { GroupName = "AuxiliaryState"; @@ -43,6 +45,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux.registerFields(GroupName, AuxMeshName); VelocityDel2Aux.registerFields(GroupName, AuxMeshName); WindForcingAux.registerFields(GroupName, AuxMeshName); + SurfTracerRestAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); } @@ -54,6 +57,7 @@ AuxiliaryState::~AuxiliaryState() { VorticityAux.unregisterFields(); VelocityDel2Aux.unregisterFields(); WindForcingAux.unregisterFields(); + SurfTracerRestAux.unregisterFields(); TracerAux.unregisterFields(); FieldGroup::destroy(GroupName); @@ -406,6 +410,16 @@ I4 AuxiliaryState::exchangeHalo() { Err += MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell); + // Performing halo exchange on individual tracers because full halo exchange + // on a 2D array assumes the first dimension is the vertical + const I4 NTracers = + SurfTracerRestAux.TracersMonthlySurfClimoCell.extent_int(0); + for (I4 LTracer = 0; LTracer < NTracers; ++LTracer) { + auto TracerSurfClimoCell = Kokkos::subview( + SurfTracerRestAux.TracersMonthlySurfClimoCell, LTracer, Kokkos::ALL); + Err += MeshHalo->exchangeFullArrayHalo(TracerSurfClimoCell, OnCell); + } + return Err; } // end exchangeHalo diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index f3c14be20ac3..1228f06962d4 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -12,6 +12,7 @@ #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" #include "auxiliaryVars/LayerThicknessAuxVars.h" +#include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" @@ -42,6 +43,7 @@ class AuxiliaryState { VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; WindForcingAuxVars WindForcingAux; + SurfTracerRestAuxVars SurfTracerRestAux; ~AuxiliaryState(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index a22e3b3b0667..02b9f9010d78 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -139,6 +139,15 @@ int initOmegaModules(MPI_Comm Comm) { PressureGrad::init(); Eos::init(); Tendencies::init(); + + // Validate SurfaceTracerRestoring configuration + Tendencies *DefTend = Tendencies::getDefault(); + if (DefTend->SurfaceTracerRestoring.Enabled && + DefTend->SurfaceTracerRestoring.NTracersToRestore == 0) { + ABORT_ERROR("OceanInit: SurfaceTracerRestoring is enabled but " + "TracersToRestore is empty"); + } + TimeStepper::init2(); Err = OceanState::init(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 3c3aacb8dde5..7788bda4363d 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -251,6 +251,61 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("SurfaceTracerRestoringEnable", + this->SurfaceTracerRestoring.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: SurfaceTracerRestoringEnable not found in TendConfig"); + if (this->SurfaceTracerRestoring.Enabled) { + Config SurfRestConfig("SurfaceRestoring"); + Err += OmegaConfig->get(SurfRestConfig); + Err += SurfRestConfig.get("PistonVelocity", + this->SurfaceTracerRestoring.PistonVelocity); + CHECK_ERROR_ABORT( + Err, + "Tendencies: PistonVelocity not found in SurfaceRestoringConfig"); + + std::vector TracersToRestore; + SurfRestConfig.get("TracersToRestore", TracersToRestore); + + // Enable restoring for specified individual tracers + I4 NumInvalidTracers = 0; + std::vector TracerIdsToRestoreVec; + for (const auto &TracerName : TracersToRestore) { + I4 TracerIndex = -1; + Tracers::getIndex(TracerIndex, TracerName); + if (TracerIndex == -1) { + LOG_ERROR("Tendencies: Tracer {} in TracersToRestore is not " + "defined", + TracerName); + NumInvalidTracers++; + } else { + bool IsDuplicate = false; + for (const auto ExistingTracerIndex : TracerIdsToRestoreVec) { + if (ExistingTracerIndex == TracerIndex) { + IsDuplicate = true; + break; + } + } + if (!IsDuplicate) { + TracerIdsToRestoreVec.push_back(TracerIndex); + } + } + } + if (NumInvalidTracers > 0) { + ABORT_ERROR("Tendencies: {} invalid tracer(s) in TracersToRestore " + "configuration", + NumInvalidTracers); + } + + this->SurfaceTracerRestoring.NTracersToRestore = + static_cast(TracerIdsToRestoreVec.size()); + this->SurfaceTracerRestoring.TracerIdsToRestore = Array1DI4( + "TracerIdsToRestore", this->SurfaceTracerRestoring.NTracersToRestore); + deepCopy(this->SurfaceTracerRestoring.TracerIdsToRestore, + HostArray1DI4(TracerIdsToRestoreVec.data(), + TracerIdsToRestoreVec.size())); + } } //------------------------------------------------------------------------------ @@ -325,7 +380,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), - CustomThicknessTend(InCustomThicknessTend), + SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { // Tendency arrays @@ -653,6 +708,7 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerHorzAdv, TracerHorzAdv); OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); + OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -746,8 +802,8 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:tracerHyperDiff", 2); } + // compute tracer tendency from vertical advection Pacer::start("Tend:computeTracerVAdvTend", 2); - // compute tracer tendencies from vertical advection Array2DReal ThicknessForVAdv; if (VAdv->VertAdvChoice == VertAdvOption::Standard) { ThicknessForVAdv = State->getLayerThickness(ThickTimeLevel); @@ -758,6 +814,25 @@ void Tendencies::computeTracerTendenciesOnly( TimeStep); Pacer::stop("Tend:computeTracerVAdvTend", 2); + // compute tracer surface restoring + const Array2DReal &TracersMonthlySurfClimo = + AuxState->SurfTracerRestAux.TracersMonthlySurfClimoCell; + const I4 NTracersToRestore = LocSurfaceTracerRestoring.NTracersToRestore; + const auto &TracerIdsToRestore = + LocSurfaceTracerRestoring.TracerIdsToRestore; + if (LocSurfaceTracerRestoring.Enabled && NTracersToRestore > 0) { + Pacer::start("Tend:surfaceTracerRestoring", 2); + parallelFor( + {NTracersToRestore, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int R, int ICell) { + const int KMin = MinLayerCell(ICell); + const int L = TracerIdsToRestore(R); + LocSurfaceTracerRestoring(LocTracerTend, L, ICell, KMin, + TracersMonthlySurfClimo, TracerArray); + }); + Pacer::stop("Tend:surfaceTracerRestoring", 2); + } + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index d5687f9274a6..96651eeacf7a 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -71,6 +71,7 @@ class Tendencies { TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; + SurfaceTracerRestoringOnCell SurfaceTracerRestoring; std::string Name; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 74dbcd033649..27b19e9f7a68 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -111,6 +111,9 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( + const HorzMesh *Mesh) {} + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const auto MaxEdges2 = Mesh->MaxEdges2; diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 55a108e8d3f5..4a1d7cbe3968 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -570,5 +570,30 @@ class TracerHyperDiffOnCell { Array1DI4 MaxLayerEdgeTop; }; +/// Surface tracer restoring term +class SurfaceTracerRestoringOnCell { + public: + bool Enabled; + Real PistonVelocity = 1.585e-5; ///< piston velocity + I4 NTracersToRestore = 0; ///< number of tracers to restore + Array1DI4 TracerIdsToRestore; ///< tracer IDs to restore + /// Need to add under sea ice restoring option when that is available + + /// constructor declaration + SurfaceTracerRestoringOnCell(const HorzMesh *Mesh); + + /// The functor takes the cell index and the array for the tracer surface + /// restoring values, outputs tendency array + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KMin, + const Array2DReal &TracersMonthlySurfClimoCell, + const Array3DReal &TracerCell) const { + + Tend(L, ICell, KMin) += + PistonVelocity * + (TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, KMin)); + } +}; + } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp new file mode 100644 index 000000000000..9a0219ff8204 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -0,0 +1,56 @@ +#include "SurfTracerRestAuxVars.h" +#include "DataTypes.h" +#include "Field.h" + +#include + +namespace OMEGA { + +SurfTracerRestAuxVars::SurfTracerRestAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, + const I4 NTracers) + : TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell" + + AuxStateSuffix, + NTracers, Mesh->NCellsSize) {} + +void SurfTracerRestAuxVars::registerFields( + const std::string &AuxGroupName, // name of Auxiliary field group + const std::string &MeshName // name of horizontal mesh +) const { + + // Create fields + const Real FillValue = -9.99e30; + int NDims = 2; + std::vector DimNames(NDims); + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + auto TracersMonthlySurfClimoCellField = + Field::create(TracersMonthlySurfClimoCell.label(), + "monthly surface tracer climatology", // long name/describe + "tracer units", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar for undefined entries + NDims, // number of dimensions + DimNames); // dim names + + // Add fields to FieldGroup + FieldGroup::addFieldToGroup(TracersMonthlySurfClimoCell.label(), + AuxGroupName); + + // Attach data + TracersMonthlySurfClimoCellField->attachData( + TracersMonthlySurfClimoCell); +} + +void SurfTracerRestAuxVars::unregisterFields() const { + Field::destroy(TracersMonthlySurfClimoCell.label()); +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h new file mode 100644 index 000000000000..c5bfc32bb43d --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h @@ -0,0 +1,27 @@ +#ifndef OMEGA_AUX_SURF_REST_H +#define OMEGA_AUX_SURF_REST_H + +#include "DataTypes.h" +#include "HorzMesh.h" +#include "HorzOperators.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +class SurfTracerRestAuxVars { + public: + Array2DReal TracersMonthlySurfClimoCell; + + SurfTracerRestAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, const I4 NTracers); + + void registerFields(const std::string &AuxGroupName, + const std::string &MeshName) const; + void unregisterFields() const; +}; + +} // namespace OMEGA +#endif diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 980a7d759dcb..a63154561c8b 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -36,6 +36,7 @@ #include #include +#include using namespace OMEGA; @@ -58,6 +59,7 @@ struct TestSetupPlane { 0.00290978146207349032}; ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, 0.00523080740758275625}; + ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedWindForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.033848740052302935, 0.01000133508329411}; @@ -197,6 +199,7 @@ struct TestSetupSphere { 0.005105510870642706}; ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, 0.0007118574326665881}; + ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedWindForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, 0.0014897009917655022}; @@ -1011,6 +1014,121 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { return Err; } // end testTracerHyperDiffOnCell +int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { + + I4 Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + + if (NTracers < 3) { + LOG_ERROR("TendencyTermsTest: SurfaceTracerRestoring requires at least 3 " + "tracers, found {}", + NTracers); + return 1; + } + + // Test multiple cases with different combinations of tracers being restored + const char *CaseLabels[3] = {"SurfaceTracerRestoringSalinityOnly", + "SurfaceTracerRestoringTemperatureOnly", + "SurfaceTracerRestoringTempSaltDebug"}; + + const std::vector> CaseTracerIds = { + {1}, // Salinity Only + {0}, // Temperature Only + {0, 1, 2}, // Temperature, Salinity, and a Debug Tracer + }; + + // Loop over cases, computing exact result, numerical result, and + // errors for each. + for (int Case = 0; Case < 3; ++Case) { + + Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); + Array2DReal InputField("InputField", Mesh->NCellsSize, 1); + Array2DReal TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell", + NTracers, Mesh->NCellsSize); + Array3DReal TracersOnCell("TracersOnCell", NTracers, Mesh->NCellsSize, + NVertLayers); + Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); + + deepCopy(ExactSurfRest, 0); + deepCopy(InputField, 0); + deepCopy(TracersMonthlySurfClimoCell, 0); + deepCopy(TracersOnCell, 0); + deepCopy(NumSurfRest, 0); + + // Set Input Field values. Use a combination of scalarB and vectorX to + // ensure the full surface tracer restoring logic is exercised. + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { + return Setup.scalarB(X, Y) + Setup.vectorX(X, Y); + }, + InputField, Geom, Mesh, OnCell); + + // Set TracersOnCell values (use scalarB for simplicity, but could be + // any field). + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + TracersOnCell, Geom, Mesh, OnCell); + parallelFor( + {NTracers, Mesh->NCellsSize, NVertLayers}, + KOKKOS_LAMBDA(int L, int ICell, int K) { + TracersOnCell(L, ICell, K) += 0.04_Real * K; + }); + + SurfaceTracerRestoringOnCell SurfRestOnC(Mesh); + SurfRestOnC.PistonVelocity = 1.585e-5; + + // Build host-selected tracer IDs for restoring and copy to device. + const I4 NTracersToRestore = static_cast(CaseTracerIds[Case].size()); + Array1DI4 TracerIdsToRestore("TracerIdsToRestore", NTracersToRestore); + HostArray1DI4 TracerIdsToRestoreH("TracerIdsToRestoreH", + NTracersToRestore); + for (I4 R = 0; R < NTracersToRestore; ++R) { + TracerIdsToRestoreH(R) = CaseTracerIds[Case][R]; + } + deepCopy(TracerIdsToRestore, TracerIdsToRestoreH); + + // Compute exact result using the same logic as + // SurfaceTracerRestoringOnCell, but iterating + // selected tracer IDs to match restoring implementation. + parallelFor( + {NTracersToRestore, Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int R, int ICell) { + const int L = TracerIdsToRestore(R); + TracersMonthlySurfClimoCell(L, ICell) = InputField(ICell, 0); + const Real Diff = TracersMonthlySurfClimoCell(L, ICell) - + TracersOnCell(L, ICell, 0); + + ExactSurfRest(L, ICell, 0) = SurfRestOnC.PistonVelocity * Diff; + }); + + // Compute numerical result + parallelFor( + {NTracersToRestore, Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int R, int ICell) { + const int L = TracerIdsToRestore(R); + SurfRestOnC(NumSurfRest, L, ICell, 0, TracersMonthlySurfClimoCell, + TracersOnCell); + }); + + ErrorMeasures SurfRestErrors; + Err += computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, + OnCell); + + Err += checkErrors("TendencyTermsTest", CaseLabels[Case], SurfRestErrors, + Setup.ExpectedSurfTrRestErrors, RTol); + } + + if (Err == 0) { + LOG_INFO("TendencyTermsTest: SurfaceTracerRestoring PASS"); + } + + return Err; +} // end testSurfaceTracerRestoringOnCell + void initTendTest(const std::string &MeshFile, int NVertLayers) { Error Err; @@ -1091,6 +1209,8 @@ int tendencyTermsTest(const std::string &MeshFile = DefaultMeshFile) { Err += testTracerHyperDiffOnCell(NVertLayers, NTracers, RTol); + Err += testSurfaceTracerRestoringOnCell(NVertLayers, NTracers, RTol); + if (Err == 0) { LOG_INFO("TendencyTermsTest: Successful completion"); } diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index aa2e2e89bbe2..27d4cdb75c67 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -231,20 +231,21 @@ int initTimeStepperTest(const std::string &mesh) { } // Disable all other tendencies - TestTendencies->ThicknessFluxDiv.Enabled = false; - TestTendencies->PotentialVortHAdv.Enabled = false; - TestTendencies->KEGrad.Enabled = false; - TestTendencies->SSHGrad.Enabled = false; - TestTendencies->VelocityDiffusion.Enabled = false; - TestTendencies->VelocityHyperDiff.Enabled = false; - TestTendencies->TracerHorzAdv.Enabled = false; - TestTendencies->TracerDiffusion.Enabled = false; - TestTendencies->TracerHyperDiff.Enabled = false; - TestTendencies->WindForcing.Enabled = false; - TestTendencies->BottomDrag.Enabled = false; - DefVAdv->ThickVertAdvEnabled = false; - DefVAdv->VelVertAdvEnabled = false; - DefVAdv->TracerVertAdvEnabled = false; + TestTendencies->ThicknessFluxDiv.Enabled = false; + TestTendencies->PotentialVortHAdv.Enabled = false; + TestTendencies->KEGrad.Enabled = false; + TestTendencies->SSHGrad.Enabled = false; + TestTendencies->VelocityDiffusion.Enabled = false; + TestTendencies->VelocityHyperDiff.Enabled = false; + TestTendencies->TracerHorzAdv.Enabled = false; + TestTendencies->TracerDiffusion.Enabled = false; + TestTendencies->TracerHyperDiff.Enabled = false; + TestTendencies->WindForcing.Enabled = false; + TestTendencies->SurfaceTracerRestoring.Enabled = false; + TestTendencies->BottomDrag.Enabled = false; + DefVAdv->ThickVertAdvEnabled = false; + DefVAdv->VelVertAdvEnabled = false; + DefVAdv->TracerVertAdvEnabled = false; return Err; }