diff --git a/components/omega/doc/devGuide/SfcCoupling.md b/components/omega/doc/devGuide/SfcCoupling.md new file mode 100644 index 000000000000..a3c5938d42a5 --- /dev/null +++ b/components/omega/doc/devGuide/SfcCoupling.md @@ -0,0 +1,83 @@ +(omega-dev-sfc-coupling)= + +# Surface Coupling + +The `SfcCoupling` class manages the variables exchanged to (`o2x`) and from +(`x2o`) the coupler. It handles import/export of the raw coupler data, +application of imported fields to the [`Forcing`](#omega-dev-forcing) object, +and accumulation of export fields over the coupling interval. It is possible +to have multiple `SfcCoupling` instances; every instance has a name and is +tracked in a static C++ map called `AllSfcCoupling`. + +## Initialization + +The static method: +```c++ +OMEGA::SfcCoupling::init(CouplingInitParams); +``` +initializes the default `SfcCoupling` given a `CouplingInitParams` struct +(number of import/export fields, name-to-index maps, coupling time step, and +coupling `Layout`, i.e. `MCT` or `MOAB`) supplied by the coupler. A pointer to +the default instance can be retrieved at any time using: +```c++ +OMEGA::SfcCoupling* DefSfcCoupling = OMEGA::SfcCoupling::getDefault(); +``` + +## Creation of non-default surface coupling objects + +A non-default instance can be created with the static `create` method, +which returns a pointer to the newly created object: +```c++ +OMEGA::SfcCoupling* NewSfcCoupling = OMEGA::SfcCoupling::create( + Name, Mesh, NImportFields, NExportFields, ImportIdxMap, ExportIdxMap, + Stepper, CouplingTimeStep, Layout); +``` +Given its name, a pointer to a named instance can be obtained at any time: +```c++ +OMEGA::SfcCoupling* NewSfcCoupling = OMEGA::SfcCoupling::get(Name); +``` + +## Data exchange + +Before import/export, unmanaged views must be attached to the raw coupler +data pointers: +```c++ +SfcCoupling.attachData(CplToOcnData, OcnToCplData); +``` +To copy data from the coupler into the `CplToOcn` fields, and from the +`OcnToCpl` fields back out to the coupler: +```c++ +SfcCoupling.importFromCoupler(); +SfcCoupling.exportToCoupler(); +``` +Imported fields are applied to a `Forcing` object with: +```c++ +SfcCoupling.applyImportFields(ForcingPtr); +``` +Export fields are accumulated (running average) each ocean time step with: +```c++ +SfcCoupling.updateExportFields(State, TracerArray); +``` + +## Units + +`OcnToCplFields`'s averaged temperature is Conservative Temperature (deg C) +on device (`AvgSfcTemperature`), invariant at all times. `copyToHost()` +converts it to in-situ temperature in Kelvin (via TEOS-10, when in use) +into the host mirror `AvgSfcTemperatureH`, once per coupling interval. +All other averaged fields keep the same units and value on host and device. + +To keep this invariant, `OcnToCplFields`'s averaged device arrays are +private: they can only be written via `updateAverages()` (called each +ocean time step from `updateExportFields()`) and read via `copyToHost()`. + +## Removal of surface coupling objects + +To erase a specific named instance use `erase`: +```c++ +OMEGA::SfcCoupling::erase(Name); +``` +To clear all instances do: +```c++ +OMEGA::SfcCoupling::clear(); +``` diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 8d1b38501df8..aee6135d35a2 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -54,6 +54,7 @@ userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv userGuide/Forcing +userGuide/SfcCoupling ``` ```{toctree} @@ -100,6 +101,7 @@ devGuide/Timing devGuide/VerticalMixingCoeff devGuide/VertAdv devGuide/Forcing +devGuide/SfcCoupling ``` ```{toctree} diff --git a/components/omega/doc/userGuide/SfcCoupling.md b/components/omega/doc/userGuide/SfcCoupling.md new file mode 100644 index 000000000000..acb1f3d7d624 --- /dev/null +++ b/components/omega/doc/userGuide/SfcCoupling.md @@ -0,0 +1,7 @@ +(omega-user-sfc-coupling)= + +## Surface Coupling + +The `SfcCoupling` class manages the exchange of surface variables with the +coupler (MCT or MOAB) between Omega and the other components of E3SM. There are +no user-configurable options. diff --git a/components/omega/src/ocn/SfcCoupling.cpp b/components/omega/src/ocn/SfcCoupling.cpp new file mode 100644 index 000000000000..42f362474c3d --- /dev/null +++ b/components/omega/src/ocn/SfcCoupling.cpp @@ -0,0 +1,385 @@ +#include "SfcCoupling.h" +#include "Eos.h" +#include "GlobalConstants.h" +#include "Logging.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "Tracers.h" +#include "VertCoord.h" + +namespace OMEGA { + +// create the static class member +SfcCoupling *SfcCoupling::DefaultSfcCoupling = nullptr; +std::map> SfcCoupling::AllSfcCoupling; + +// Initalize the surface coupling. Assumes the default HorzMesh and +// TimeStepper have been initialized +int SfcCoupling::init(const CouplingInitParams &CouplingInitParams) { + + int Err = 0; // default successful return code + + // Retrieve the default horizontal mesh and timestepper + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); + auto *DefTimeStepper = TimeStepper::getDefault(); + + TimeInterval OcnTimeStep = DefTimeStepper->getTimeStep(); + TimeInterval CplTimeStep = CouplingInitParams.CouplingTimeStep; + + R8 CplTimeStepSeconds, OcnTimeStepSeconds; + CplTimeStep.get(CplTimeStepSeconds, TimeUnits::Seconds); + OcnTimeStep.get(OcnTimeStepSeconds, TimeUnits::Seconds); + + if (CplTimeStepSeconds < OcnTimeStepSeconds) { + LOG_ERROR("Coupling interval is: {} (seconds)", CplTimeStepSeconds); + LOG_ERROR("Ocean timestep is: {} (seconds)", OcnTimeStepSeconds); + ABORT_ERROR( + "The ocean timestep cannot be longer than the coupling interval."); + } + + if (std::fmod(CplTimeStepSeconds, OcnTimeStepSeconds) > 1e-10) { + LOG_ERROR("Coupling interval is: {} (seconds)", CplTimeStepSeconds); + LOG_ERROR("Ocean timestep is: {} (seconds)", OcnTimeStepSeconds); + ABORT_ERROR( + "Coupling interval must be evenly divisible by the ocean timestep."); + } + + // Create the default surface coupling object and set pointer to it + SfcCoupling::DefaultSfcCoupling = SfcCoupling::create( + "Default", DefHorzMesh, CouplingInitParams.NImportFields, + CouplingInitParams.NExportFields, CouplingInitParams.ImportIdxMap, + CouplingInitParams.ExportIdxMap, DefTimeStepper, CplTimeStep, + CouplingInitParams.Layout); + + return Err; +} + +// Construct a new surface coupling object +SfcCoupling::SfcCoupling(const std::string &Name_, const HorzMesh *Mesh, + const int NImportFields_, const int NExportFields_, + const std::map &ImportIdxMap, + const std::map &ExportIdxMap, + TimeStepper *Stepper, + const TimeInterval &CouplingTimeStep, + const CouplingLayout &Layout) + : Name(Name_), NImportFields(NImportFields_), NExportFields(NExportFields_), + ImportIdxMap(ImportIdxMap), ExportIdxMap(ExportIdxMap), + CplToOcn(Name_, Mesh), OcnToCpl(Name_, Mesh), Layout(Layout) { + + // Retrieve mesh cell count + NCellsOwned = Mesh->NCellsOwned; + + NAccumSteps = 0; + + // Allocate variables on stack for creating the CouplingAlarm + std::string AlarmName = "CouplingAlarm"; + Clock *StepperClock = Stepper->getClock(); + TimeInstant StartTime = Stepper->getStartTime(); + + // Avoid alarm name collisions on the shared clock for non-default instances + if (Name_ != "Default") + AlarmName += Name_; + + CouplingAlarm = Alarm(AlarmName, CouplingTimeStep, StartTime); + StepperClock->attachAlarm(&CouplingAlarm); +} + +// Create a new surface coupling object by calling the constructor and storing +// it in the AllSfcCoupling map +SfcCoupling *SfcCoupling::create( + const std::string &Name, const HorzMesh *Mesh, const int NImportFields, + const int NExportFields, const std::map &ImportIdxMap, + const std::map &ExportIdxMap, TimeStepper *Stepper, + const TimeInterval &CouplingTimeStep, const CouplingLayout &Layout) { + + // Check to see if a surface coupling of the same name already exists + if (AllSfcCoupling.find(Name) != AllSfcCoupling.end()) { + LOG_ERROR("Attempted to create a SfcCoupling with name {}, but a " + "SfcCoupling with that name already exists", + Name); + return nullptr; + } + + // create a new surface coupling on the heap and store it in the map of + // unique_ptrs, which will manage its lifetime + auto *NewSfcCoupling = + new SfcCoupling(Name, Mesh, NImportFields, NExportFields, ImportIdxMap, + ExportIdxMap, Stepper, CouplingTimeStep, Layout); + AllSfcCoupling.emplace(Name, NewSfcCoupling); + + return NewSfcCoupling; +} // end SfcCoupling create + +// Get the default surface coupling object +SfcCoupling *SfcCoupling::getDefault() { + return SfcCoupling::DefaultSfcCoupling; +} + +// Get a surface coupling object by name +SfcCoupling *SfcCoupling::get(const std::string Name) { + // look for an instance of this name + auto it = AllSfcCoupling.find(Name); + + // if found, return the pointer + if (it != AllSfcCoupling.end()) { + return it->second.get(); + + // otherwise print error message and return nullptr + } else { + LOG_ERROR("SfcCoupling::get: Attempted to retrieve non-existent " + "surface coupling object:"); + LOG_ERROR("{} has not been defined or has been removed", Name); + return nullptr; + } +} + +// Destructor +SfcCoupling::~SfcCoupling() {} + +// Remove surface coupling object by name +void SfcCoupling::erase(const std::string Name) { AllSfcCoupling.erase(Name); } + +// Remove all surface coupling objects +void SfcCoupling::clear() { + AllSfcCoupling.clear(); + DefaultSfcCoupling = nullptr; // prevent dangling pointer +} + +// Getter for private member NAccumSteps +I4 SfcCoupling::getNAccumSteps() const { return NAccumSteps; } + +// Create views of the raw coupling data arrays +void SfcCoupling::attachData(const Real *CplToOcnData, Real *OcnToCplData) { + + // Kokkos::LayoutStride index math uses a runtime stride value, rather than + // a compile-time-optimized stride value. Can switch to ifdefs if this + // becomes a performance bottleneck + Kokkos::LayoutStride CplToOcnLayout, OcnToCplLayout; + + if (Layout == CouplingLayout::MCT) { + /// MCT layout: (NCellsOwned, NImportFields) - field idx strides faster + CplToOcnLayout = + Kokkos::LayoutStride(NImportFields, 1, NCellsOwned, NImportFields); + OcnToCplLayout = + Kokkos::LayoutStride(NExportFields, 1, NCellsOwned, NExportFields); + } else if (Layout == CouplingLayout::MOAB) { + /// MOAB layout: (NImportFields, NCellsOwned) - cell idx strides faster + CplToOcnLayout = + Kokkos::LayoutStride(NImportFields, NCellsOwned, NCellsOwned, 1); + OcnToCplLayout = + Kokkos::LayoutStride(NExportFields, NCellsOwned, NCellsOwned, 1); + } else { + ABORT_ERROR("SfcCoupling::attachData: Unknown coupling layout"); + } + + CplToOcnView = decltype(CplToOcnView)(CplToOcnData, CplToOcnLayout); + OcnToCplView = decltype(OcnToCplView)(OcnToCplData, OcnToCplLayout); +} + +void SfcCoupling::importFromCoupler() { + + if (CplToOcnView.data() == nullptr) { + ABORT_ERROR( + "CplToOcnView is not attached to data. The SfcCoupling::attachData " + "method must be called before importing data from the coupler."); + } + + // Get import field indices for surface stress components + int TauxIdx = ImportIdxMap.at("Foxx_taux"); + int TauyIdx = ImportIdxMap.at("Foxx_tauy"); + + // Copy Kokkos view handles + auto CplToOcnView_ = CplToOcnView; + auto SfcStressZonal_ = CplToOcn.SfcStressZonal; + auto SfcStressMerid_ = CplToOcn.SfcStressMerid; + + /// TODO: Shouldn't be making direct calls to Kokkos here. + /// How often is threading used? Becuase this will be a serial loop + /// unless threading is used. But this has to be run on the host. + auto Policy = Kokkos::RangePolicy>( + 0, NCellsOwned); + Kokkos::parallel_for("importFromCoupler", Policy, [=](int Idx) { + SfcStressZonal_(Idx) = CplToOcnView_(TauxIdx, Idx); + SfcStressMerid_(Idx) = CplToOcnView_(TauyIdx, Idx); + }); +} + +void SfcCoupling::exportToCoupler() { + + if (OcnToCplView.data() == nullptr) { + ABORT_ERROR( + "OcnToCplView is not attached to data. The SfcCoupling::attachData " + "method must be called before exporting data to the coupler."); + } + + // Copy the OcnToCpl fields to their host mirrors + OcnToCpl.copyToHost(); + + int TempIdx = ExportIdxMap.at("So_t"); + int SalinIdx = ExportIdxMap.at("So_s"); + int VelUIdx = ExportIdxMap.at("So_u"); + int VelVIdx = ExportIdxMap.at("So_v"); + int SshIdx = ExportIdxMap.at("So_ssh"); + + // Copy Kokkos view handles + auto OcnToCplView_ = OcnToCplView; + auto AvgSfcTemperature_ = OcnToCpl.AvgSfcTemperatureH; + auto AvgSfcSalinity_ = OcnToCpl.AvgSfcSalinityH; + auto AvgSfcVelocityZonal_ = OcnToCpl.AvgSfcVelocityZonalH; + auto AvgSfcVelocityMerid_ = OcnToCpl.AvgSfcVelocityMeridH; + auto InstSshCellH_ = OcnToCpl.InstSshCellH; + + /// TODO: Shouldn't be making direct calls to Kokkos here. + auto Policy = Kokkos::RangePolicy>( + 0, NCellsOwned); + Kokkos::parallel_for("exportToCoupler", Policy, [=](int Idx) { + OcnToCplView_(TempIdx, Idx) = AvgSfcTemperature_(Idx); + OcnToCplView_(SalinIdx, Idx) = AvgSfcSalinity_(Idx); + OcnToCplView_(VelUIdx, Idx) = AvgSfcVelocityZonal_(Idx); + OcnToCplView_(VelVIdx, Idx) = AvgSfcVelocityMerid_(Idx); + OcnToCplView_(SshIdx, Idx) = InstSshCellH_(Idx); + }); + + OcnToCpl.resetFields(); // Reset fields to 0 for the next coupling interval + NAccumSteps = 0; // Reset step counter for the next coupling interval +} +void SfcCoupling::applyImportFields(Forcing *Forcing) { + + // Copy the SfcCoupling host arrays into the Forcing device arrays. + // Copy is only done over the owned cells, since thats all the SfcCoupling + // data is defined over. Forcing will be responsible for halo exchanges. + deepCopy(ownedSubView(Forcing->SfcStressForcing.ZonalStressCell), + CplToOcn.SfcStressZonal); + deepCopy(ownedSubView(Forcing->SfcStressForcing.MeridStressCell), + CplToOcn.SfcStressMerid); +}; + +void SfcCoupling::updateExportFields(const OceanState *State, + const Array3DReal &TracerArray) { + + OcnToCpl.updateFields(State, TracerArray, NAccumSteps, NCellsOwned); + + NAccumSteps++; +} + +CplToOcnFields::CplToOcnFields(const std::string &Suffix, const HorzMesh *Mesh) + : SfcStressZonal("SfcStressZonal" + Suffix, Mesh->NCellsOwned), + SfcStressMerid("SfcStressMeridional" + Suffix, Mesh->NCellsOwned) {} + +OcnToCplFields::OcnToCplFields(const std::string &Suffix, const HorzMesh *Mesh) + : AvgSfcTemperature("AvgSfcTemperature" + Suffix, Mesh->NCellsOwned), + AvgSfcSalinity("AvgSfcSalinity" + Suffix, Mesh->NCellsOwned), + AvgSfcVelocityZonal("AvgSfcVelocityZonal" + Suffix, Mesh->NCellsOwned), + AvgSfcVelocityMerid("AvgSfcVelocityMeridional" + Suffix, + Mesh->NCellsOwned), + InstSshCellH("InstSshCellH" + Suffix, Mesh->NCellsOwned), + InSituTempScratch("InSituTempScratch" + Suffix, Mesh->NCellsOwned) { + + // Kokkok views created with a label are zero-initialized by default. + // We reset the fields here anyway to be explicit about the fact that the + // OcnToCpl fields need to begin a coupling interval with all zeros. + resetFields(); + + AvgSfcTemperatureH = createHostMirrorCopy(AvgSfcTemperature); + AvgSfcSalinityH = createHostMirrorCopy(AvgSfcSalinity); + AvgSfcVelocityZonalH = createHostMirrorCopy(AvgSfcVelocityZonal); + AvgSfcVelocityMeridH = createHostMirrorCopy(AvgSfcVelocityMerid); +} + +void OcnToCplFields::updateFields(const OceanState *State, + const Array3DReal &TracerArray, + const I4 NAccumSteps, const I4 NCellsOwned) { + + I4 TemperatureIdx, SalinityIdx; + Tracers::getIndex(TemperatureIdx, "Temperature"); + Tracers::getIndex(SalinityIdx, "Salinity"); + + auto Temperature = + Kokkos::subview(TracerArray, TemperatureIdx, Kokkos::ALL, Kokkos::ALL); + auto Salinity = + Kokkos::subview(TracerArray, SalinityIdx, Kokkos::ALL, Kokkos::ALL); + + VertCoord *DefVertCoord = VertCoord::getDefault(); + + OMEGA_SCOPE(LocMinLayerCell, DefVertCoord->MinLayerCell); + OMEGA_SCOPE(LocAvgSfcSalinity, AvgSfcSalinity); + OMEGA_SCOPE(LocAvgSfcTemp, AvgSfcTemperature); + OMEGA_SCOPE(LocAvgSfcVelZonal, AvgSfcVelocityZonal); + OMEGA_SCOPE(LocAvgSfcVelMerid, AvgSfcVelocityMerid); + + // TODO: Implement vector reconsturction for velocity field. + constexpr Real ConstSfcVelocity = 1e-4; + + parallelFor( + {NCellsOwned}, KOKKOS_LAMBDA(int ICell) { + const int KSfc = LocMinLayerCell(ICell); + + // Update the averaged fields using Welford's online algorithm + LocAvgSfcTemp(ICell) = updateAverage( + LocAvgSfcTemp(ICell), Temperature(ICell, KSfc), NAccumSteps); + + LocAvgSfcSalinity(ICell) = updateAverage( + LocAvgSfcSalinity(ICell), Salinity(ICell, KSfc), NAccumSteps); + + LocAvgSfcVelZonal(ICell) = updateAverage( + LocAvgSfcVelZonal(ICell), ConstSfcVelocity, NAccumSteps); + + LocAvgSfcVelMerid(ICell) = updateAverage( + LocAvgSfcVelMerid(ICell), ConstSfcVelocity, NAccumSteps); + }); +} + +// Conversion from conservative [C] to in-situ [K] temperature is done on +// the device, as well as conversion from practical to absolute salinity. +// This allows the expensive unit conversions to be called on device, once +// per coupling interval. Conversion is written to scratch buffers, to +// keep units consistent in time. Special care is paid to guard external +// access to the device arrays, so that inconsistent units between the device +// and host mirrors are not exposed to the rest of the code. +void OcnToCplFields::copyToHost() { + + // Convert averaged Conservative Temp to in-situ (approx by potential + // temp at P=0) Kelvin, once per coupling interval, on device. Written + // into a scratch buffer so AvgSfcTemperature always stays deg C. + // A local Teos10Eos is constructed here rather than reusing Eos's + // instance: calcPtFromCt() needs no config-derived state (fixed + // polynomial coefficients only), so this avoids exposing Eos internals. + EosType LocEosChoice = Eos::getInstance()->EosChoice; + Teos10Eos LocTeos10(VertCoord::getDefault()); + OMEGA_SCOPE(LocAvgSfcTemp, AvgSfcTemperature); + OMEGA_SCOPE(LocAvgSfcSalinity, AvgSfcSalinity); + OMEGA_SCOPE(LocInSituTemp, InSituTempScratch); + + parallelFor( + {(int)AvgSfcTemperature.extent(0)}, KOKKOS_LAMBDA(int Cell) { + const Real Ct = LocAvgSfcTemp(Cell); + const Real Sa = LocAvgSfcSalinity(Cell); + const Real Pt = LocEosChoice == EosType::Teos10Eos + ? LocTeos10.calcPtFromCt(Sa, Ct) + : Ct; + LocInSituTemp(Cell) = Pt + TkFrz; // C to K temperature conversion + }); + + deepCopy(AvgSfcTemperatureH, InSituTempScratch); + deepCopy(AvgSfcSalinityH, AvgSfcSalinity); + deepCopy(AvgSfcVelocityZonalH, AvgSfcVelocityZonal); + deepCopy(AvgSfcVelocityMeridH, AvgSfcVelocityMerid); + + // SSH is an instantaneous field, so we don't bother with a device mirror of + // our own. Instead, copy from the VertCoord, which owns SSH, host array. + VertCoord *DefVertCoord = VertCoord::getDefault(); + + auto SSHCellOwned = Kokkos::subview( + DefVertCoord->SshCell, std::make_pair(0, (int)InstSshCellH.extent(0))); + + deepCopy(InstSshCellH, SSHCellOwned); +} + +// OcnToCpl fields need to begin a coupling interval with all values set to 0. +void OcnToCplFields::resetFields() { + deepCopy(AvgSfcTemperature, 0.0_Real); + deepCopy(AvgSfcSalinity, 0.0_Real); + deepCopy(AvgSfcVelocityZonal, 0.0_Real); + deepCopy(AvgSfcVelocityMerid, 0.0_Real); +} +} // namespace OMEGA diff --git a/components/omega/src/ocn/SfcCoupling.h b/components/omega/src/ocn/SfcCoupling.h new file mode 100644 index 000000000000..85a63bd44860 --- /dev/null +++ b/components/omega/src/ocn/SfcCoupling.h @@ -0,0 +1,218 @@ +#ifndef OMEGA_SURFACECOUPLING_H +#define OMEGA_SURFACECOUPLING_H +//===-- ocn/SfcCouling.h - surface coupling ----------------*- C++ -*-===// +// +/// \file +/// \brief Contains the coupling variables exchanged with the coupler +/// +/// The SfcCouling class contains the variables exchanged with the coupler +/// for a sub-domain of the global horizontal mesh. +// +//===----------------------------------------------------------------------===// + +#include "DataTypes.h" +#include "Forcing.h" +#include "HorzMesh.h" +#include "OceanState.h" +#include "TimeMgr.h" +#include "TimeStepper.h" + +#include + +namespace OMEGA { + +// Welford's online algorithm for computing a running average +KOKKOS_INLINE_FUNCTION Real updateAverage(const Real OldAvg, + const Real NewValue, + const I4 NAccumSteps) { + return OldAvg + (NewValue - OldAvg) / (NAccumSteps + 1); +} + +enum class CouplingLayout { MCT, MOAB }; + +// Parameters needed to initialize a SfcCoupling object. The information +// needed to initialize these parameters is provided by the coupler. +struct CouplingInitParams { + int NImportFields; + int NExportFields; + std::map ImportIdxMap; + std::map ExportIdxMap; + TimeInterval CouplingTimeStep; + CouplingLayout Layout; +}; + +// x2o: Coupler to Ocean +class CplToOcnFields { + public: + // x2o fields only need to be stored on the host. + // The SfcCoupling::applyImportFields() method will handle copying the + // data to the device. + HostArray1DReal SfcStressZonal; ///< Foxx_taux [N m^-2] + HostArray1DReal SfcStressMerid; ///< Foxx_tauy [N m^-2] + + CplToOcnFields(const std::string &Suffix, const HorzMesh *Mesh); +}; + +// o2x: Ocean to Coupler +class OcnToCplFields { + public: + ///< So_t [K], in-situ approx (potential temp at P=0) + HostArray1DReal AvgSfcTemperatureH; + + /// TODO: Export practical salinity (unitless) to coupler + ///< So_s [g kg^-1], absolute salinity + HostArray1DReal AvgSfcSalinityH; + + ///< So_u [m s^-1] + HostArray1DReal AvgSfcVelocityZonalH; + + ///< So_v [m s^-1] + HostArray1DReal AvgSfcVelocityMeridH; + + ///< So_ssh [m] + /// instantaneous field, so no device mirror is needed + HostArray1DReal InstSshCellH; + + // Accumulate one ocean timestep's contribution to the running averages + void updateFields(const OceanState *State, const Array3DReal &TracerArray, + I4 NAccumSteps, I4 NCellsOwned); + + // Copy device arrays into their host mirrors and do unit conversion. + void copyToHost(); + + // Reset all fields to 0 + void resetFields(); + + OcnToCplFields(const std::string &Suffix, const HorzMesh *Mesh); + + private: + // Device arrays are private to prevent any code from mirroring to host + // except through copyToHost(), so that the different unit bewteen the + // device and host mirrors of temperature and salinity are not exposed to + // the rest of the code. + Array1DReal AvgSfcTemperature; // [C], conservative temperature + Array1DReal AvgSfcSalinity; // [g kg^-1], absolute salinity + Array1DReal AvgSfcVelocityZonal; + Array1DReal AvgSfcVelocityMerid; + + // Scratch buffer for the in-situ Kelvin conversion in copyToHost() + Array1DReal InSituTempScratch; // [K], in-situ approx (potential temp at P=0) +}; + +/// A class for interfacing with the coupler + +/// The SfcCoupling class provides a container for the variables exchanged +/// to (o2x) and from (x2o) the coupler. +class SfcCoupling { + + private: + static SfcCoupling *DefaultSfcCoupling; + + static std::map> AllSfcCoupling; + + // Construct a new local coupling object + SfcCoupling(const std::string &Name_, const HorzMesh *Mesh, + const int NImportFields_, const int NExportFields_, + const std::map &ImportIdxMap, + const std::map &ExportIdxMap, + TimeStepper *Stepper, const TimeInterval &CouplingTimeStep, + const CouplingLayout &Layout); + + // Forbid copy and move construction + SfcCoupling(const SfcCoupling &) = delete; + SfcCoupling(SfcCoupling &&) = delete; + + // Create subview that only include the owned cells + template auto ownedSubView(const View &V) const { + return Kokkos::subview(V, std::make_pair(0, NCellsOwned)); + } + + // Number of ocn timesteps acccumulated over the coupling interval + I4 NAccumSteps; + + CouplingLayout Layout; ///< Coupling layout (MCT or MOAB) + + // Map of import/export variable names to index in the raw data arrays + std::map ImportIdxMap; + std::map ExportIdxMap; + + public: + std::string Name; + + I4 NCellsOwned; ///< Number of cells owned by this task + + // The values below will be larger than InportIdx.size() and + // ExportIdxMap.size() because omega does not ingest all cpl fields (e.g. + // BGC, landice), yet... + I4 NImportFields; ///< Num of fields in the x2o pointer array + I4 NExportFields; ///< Num of fields in the o2x pointer array + + // Coupling Variable containers + CplToOcnFields CplToOcn; ///< Coupler to Ocean (x2o) + OcnToCplFields OcnToCpl; ///< Ocean to Coupler (o2x) + + Alarm CouplingAlarm; ///< Alarm for coupling interval + + /// View of Coupler to Ocean (x2o) raw data + Kokkos::View> + CplToOcnView; + + /// View of Ocean to Coupler (o2x) raw data + Kokkos::View> + OcnToCplView; + + // Methods + + /// Create a new surface coupling by calling the constructor and put it + /// in the AllSfcCoupling map + static SfcCoupling *create(const std::string &Name, const HorzMesh *Mesh, + const int NImportFields, const int NExportFields, + const std::map &ImportIdxMap, + const std::map &ExportIdxMap, + TimeStepper *Stepper, + const TimeInterval &CouplingTimeStep, + const CouplingLayout &CouplingLayout); + + /// Initialize SfcCoupling + static int init(const CouplingInitParams &CouplingInitParams); + + /// Destructor - deallocates all memory and deletes an SfcCoupling + ~SfcCoupling(); + + /// Dealllocates arrays + static void clear(); + + /// Remove surface coupling object by name + static void erase(const std::string InName); + + /// Get the default surface coupling object + static SfcCoupling *getDefault(); + + /// Get a surface coupling object by name + static SfcCoupling *get(const std::string name); + + /// Getter for the number of ocean timesteps accumulated over the coupling + /// interval + I4 getNAccumSteps() const; + + /// Create views of the coupling data arrays + void attachData(const Real *CplToOcnData, Real *OcnToCplData); + + /// Import data from the unmanaged view of x2o pointer into OcnToCpl object + void importFromCoupler(); + + /// Export data from OcnToCpl object into the unmanaged view of o2x pointer + void exportToCoupler(); + + /// Apply the imported data to the Forcing object + void applyImportFields(Forcing *Forcing); + + /// Update the export fields + void updateExportFields(const OceanState *State, + const Array3DReal &TracerArray); +}; + +} // end namespace OMEGA +#endif // defined OMEGA_SURFACECOUPLING_H diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index b30760fea72b..535e7c71aefd 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -555,3 +555,14 @@ add_omega_test( ocn/StateValidationTest.cpp "-n;8" ) + +##################### +# SurfaceCouling test +##################### + +add_omega_test( + SFCCOUPLING_TEST + testSfcCoupling.exe + ocn/SfcCouplingTest.cpp + "-n;8" +) diff --git a/components/omega/test/ocn/SfcCouplingTest.cpp b/components/omega/test/ocn/SfcCouplingTest.cpp new file mode 100644 index 000000000000..78964f622364 --- /dev/null +++ b/components/omega/test/ocn/SfcCouplingTest.cpp @@ -0,0 +1,591 @@ +#include "SfcCoupling.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Eos.h" +#include "Field.h" +#include "Forcing.h" +#include "GlobalConstants.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanTestCommon.h" +#include "OmegaKokkos.h" +#include "Pacer.h" +#include "TimeStepper.h" +#include "VertCoord.h" +#include "mpi.h" + +using namespace OMEGA; + +struct TestSetup { + + std::map ImportIdxMap = {{"Foxx_taux", 3}, + {"Foxx_tauy", 8}}; + std::map ExportIdxMap = { + {"So_t", 2}, {"So_s", 4}, {"So_u", 6}, {"So_v", 9}, {"So_ssh", 1}}; +}; + +CouplingInitParams mockCouplingInitParams( + CouplingLayout Layout = CouplingLayout::MCT, + std::optional CouplingTimeStep = std::nullopt) { + + TestSetup Setup; + + auto *DefTimeStepper = TimeStepper::getDefault(); + TimeInterval CouplingTimeStep_ = + CouplingTimeStep.value_or(DefTimeStepper->getTimeStep()); + + CouplingInitParams CouplingParams{.NImportFields = 10, + .NExportFields = 10, + .ImportIdxMap = Setup.ImportIdxMap, + .ExportIdxMap = Setup.ExportIdxMap, + .CouplingTimeStep = CouplingTimeStep_, + .Layout = Layout}; + + return CouplingParams; +} + +std::string toString(const CouplingLayout &Layout) { + switch (Layout) { + case CouplingLayout::MCT: + return "MCT"; + case CouplingLayout::MOAB: + return "MOAB"; + default: + return "Unknown"; + } +} + +HostArray1DReal makeCellVarryingArray(const std::string &Name, const int NCells, + const Real Base) { + HostArray1DReal Array(Name, NCells); + for (int Cell = 0; Cell < NCells; Cell++) { + Array(Cell) = Base + Cell; + } + return Array; +} + +// Shared index formula for packing/unpacking from raw coupler buffer +int flatIdx(const CouplingLayout Layout, const int Cell, const int Field, + const int NCells, const int NFields) { + if (Layout == CouplingLayout::MCT) + return Cell * NFields + Field; + else // MOAB + return Field * NCells + Cell; +} + +int initSfcCouplingTest(const std::string &MeshFile) { + + int Err = 0; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + LOG_INFO("------ Surface Coupling unit tests ------"); + + // Open config file + Config("Omega"); + Config::readAll("omega.yml"); + + // Force ConstantEos so calcPtFromCt() is a no-op, letting tests below use + // simple mocked Temperature/Salinity without the TEOS-10 polynomial + Config *OmegaConfig = Config::getOmegaConfig(); + Config EosConfig("Eos"); + OmegaConfig->get(EosConfig); + EosConfig.set("EosType", std::string("constant")); + + // Initialize time stepper and get model clock + TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); + + IO::init(DefComm); + Decomp::init(MeshFile); + + // Initialize streams + Field::init(ModelClock); + IOStream::init(ModelClock); + + // TODO: Need to initialize halo? SfcCoupling has no lateral exchanges + int HaloErr = Halo::init(); + if (HaloErr != 0) { + Err++; + LOG_ERROR("SfcCouplingTest: Error initializing default halo"); + } + + HorzMesh::init(ModelClock); + + VertCoord::init(); + + int StateErr = OceanState::init(); + if (StateErr != 0) { + Err++; + LOG_ERROR("SfcCouplingTest: Error initializing defualt ocean state"); + } + + Eos::init(); + + Forcing::init(); + Tracers::init(); + + return Err; +} + +int testImportFromCoupler(const CouplingLayout Layout) { + + int Err = 0; + + auto CouplingParams = mockCouplingInitParams(Layout); + Err += SfcCoupling::init(CouplingParams); + + SfcCoupling *DefCoupling = SfcCoupling::getDefault(); + + int NCells = DefCoupling->NCellsOwned; + int NImports = DefCoupling->NImportFields; + int NExports = DefCoupling->NExportFields; + + int TauxIdx = CouplingParams.ImportIdxMap.at("Foxx_taux"); + int TauyIdx = CouplingParams.ImportIdxMap.at("Foxx_tauy"); + + std::vector CplToOcnData(NCells * NImports, 0.0); + std::vector OcnToCplData(NCells * NExports, 0.0); + + HostArray1DReal ExpectedSfcStressZonal = + makeCellVarryingArray("ExpectedSfcStressZonal", NCells, Real(TauxIdx)); + HostArray1DReal ExpectedSfcStressMerid = + makeCellVarryingArray("ExpectedSfcStressMerid", NCells, Real(TauyIdx)); + + for (int Cell = 0; Cell < NCells; Cell++) { + CplToOcnData[flatIdx(Layout, Cell, TauxIdx, NCells, NImports)] = + ExpectedSfcStressZonal(Cell); + CplToOcnData[flatIdx(Layout, Cell, TauyIdx, NCells, NImports)] = + ExpectedSfcStressMerid(Cell); + } + + DefCoupling->attachData(CplToOcnData.data(), OcnToCplData.data()); + DefCoupling->importFromCoupler(); + + auto ImportPass = arraysEqual(DefCoupling->CplToOcn.SfcStressZonal, + ExpectedSfcStressZonal) && + arraysEqual(DefCoupling->CplToOcn.SfcStressMerid, + ExpectedSfcStressMerid); + + if (ImportPass) { + LOG_INFO("SfcCouplingTest: importFromCoupler with {} layout PASS", + toString(Layout)); + } else { + Err++; + LOG_ERROR("SfcCouplingTest: importFromCoupler with {} layout FAIL", + toString(Layout)); + } + + SfcCoupling::clear(); + + return Err; +} + +int testApplyImportFields() { + + int Err = 0; + + auto CouplingParams = mockCouplingInitParams(); + Err += SfcCoupling::init(CouplingParams); + + Forcing *DefForcing = Forcing::getDefault(); + SfcCoupling *DefCoupling = SfcCoupling::getDefault(); + + int NCells = DefCoupling->NCellsOwned; + int NImports = DefCoupling->NImportFields; + int NExports = DefCoupling->NExportFields; + + int TauxIdx = CouplingParams.ImportIdxMap.at("Foxx_taux"); + int TauyIdx = CouplingParams.ImportIdxMap.at("Foxx_tauy"); + + std::vector CplToOcnData(NCells * NImports, 0.0); + std::vector OcnToCplData(NCells * NExports, 0.0); + + int Offset = 27; + HostArray1DReal ExpectedSfcStressZonal = makeCellVarryingArray( + "ExpectedSfcStressZonal", NCells, Real(TauxIdx + Offset)); + HostArray1DReal ExpectedSfcStressMerid = makeCellVarryingArray( + "ExpectedSfcStressMerid", NCells, Real(TauyIdx + Offset)); + + // Copy the expected values into the CplToOcn fields directly + deepCopy(DefCoupling->CplToOcn.SfcStressZonal, ExpectedSfcStressZonal); + deepCopy(DefCoupling->CplToOcn.SfcStressMerid, ExpectedSfcStressMerid); + + DefCoupling->applyImportFields(DefForcing); + + auto SfcStressZonalOwned = Kokkos::subview( + DefForcing->SfcStressForcing.ZonalStressCell, std::pair(0, NCells)); + auto SfcStressMeridOwned = Kokkos::subview( + DefForcing->SfcStressForcing.MeridStressCell, std::pair(0, NCells)); + + auto ApplyPass = arraysEqual(SfcStressZonalOwned, ExpectedSfcStressZonal) && + arraysEqual(SfcStressMeridOwned, ExpectedSfcStressMerid); + + if (ApplyPass) { + LOG_INFO("SfcCouplingTest: applyImportFields PASS"); + } else { + Err++; + LOG_ERROR("SfcCouplingTest: applyImportFields FAIL"); + } + + SfcCoupling::clear(); + + return Err; +} + +int testUpdateExportFields(const I4 NSteps) { + + int Err = 0; + + auto *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); + + // Reset the shared clock + ModelClock->setCurrentTime(DefStepper->getStartTime()); + + // Coupling interval spans NSteps ocean timesteps + auto CouplingParams = mockCouplingInitParams( + CouplingLayout::MCT, DefStepper->getTimeStep() * NSteps); + + Err += SfcCoupling::init(CouplingParams); + + SfcCoupling *DefCoupling = SfcCoupling::getDefault(); + OceanState *DefState = OceanState::getDefault(); + VertCoord *DefVertCoord = VertCoord::getDefault(); + + int NCells = DefCoupling->NCellsOwned; + + Real TempBase = static_cast(CouplingParams.ExportIdxMap.at("So_t")); + Real SalinBase = static_cast(CouplingParams.ExportIdxMap.at("So_s")); + + I4 TempIdx, SalinIdx; + Tracers::getIndex(TempIdx, "Temperature"); + Tracers::getIndex(SalinIdx, "Salinity"); + + while (!DefCoupling->CouplingAlarm.isRinging()) { + Real CurrStep = static_cast(DefCoupling->getNAccumSteps()); + + HostArray2DReal TempH = Tracers::getHostByIndex(0, TempIdx); + HostArray2DReal SalinH = Tracers::getHostByIndex(0, SalinIdx); + + for (int Cell = 0; Cell < NCells; Cell++) { + int KSfc = DefVertCoord->MinLayerCellH(Cell); + + TempH(Cell, KSfc) = TempBase + Cell + CurrStep; + SalinH(Cell, KSfc) = SalinBase + Cell + CurrStep; + } + Tracers::copyToDevice(0); + + DefCoupling->updateExportFields(DefState, Tracers::getAll(0)); + + ModelClock->advance(); + } + + // Sanity check: alarm should ring after NSteps + if (DefCoupling->getNAccumSteps() != NSteps) { + Err++; + LOG_ERROR("SfcCouplingTest: updateExportFields FAIL - " + "NAccumSteps = {}, expected {}", + DefCoupling->getNAccumSteps(), NSteps); + } + + // copyToHost() is the only sanctioned way to read the averages; it also + // converts temp to in-situ Kelvin (identity CT->PT w/ ConstantEos) + DefCoupling->OcnToCpl.copyToHost(); + auto &AvgTempH = DefCoupling->OcnToCpl.AvgSfcTemperatureH; + auto &AvgSalinH = DefCoupling->OcnToCpl.AvgSfcSalinityH; + + // Only accumulated round-off error (no discretization/truncation error), + // so use a tighter tolerance than tests comparing against analytic + // solutions of numerical operators. + Real RTol = sizeof(Real) == 4 ? 1e-5 : 1e-10; + Real StepOffset = static_cast(NSteps - 1) / 2.0; + HostArray1DReal ExpectedTemp = makeCellVarryingArray( + "ExpectedTemp", NCells, TempBase + StepOffset + TkFrz); + HostArray1DReal ExpectedSalin = + makeCellVarryingArray("ExpectedSalin", NCells, SalinBase + StepOffset); + + I4 TempErr = 0; + I4 SalinErr = 0; + + for (int Cell = 0; Cell < NCells; Cell++) { + if (!isApprox(AvgTempH(Cell), ExpectedTemp(Cell), RTol)) { + TempErr++; + } + + if (!isApprox(AvgSalinH(Cell), ExpectedSalin(Cell), RTol)) { + SalinErr++; + } + } + + if (TempErr == 0) { + LOG_INFO("SfcCouplingTest: updateExportFields PASS - " + "AvgSfcTemperature within relative tolerance of {}", + RTol); + } else { + Err += TempErr; + LOG_ERROR("SfcCouplingTest: updateExportFields FAIL - " + "AvgSfcTemperature outside relative tolerance of {} for {} " + "cells", + RTol, TempErr); + } + + if (SalinErr == 0) { + LOG_INFO("SfcCouplingTest: updateExportFields PASS - " + "AvgSfcSalinity within relative tolerance of {}", + RTol); + } else { + Err += SalinErr; + LOG_ERROR("SfcCouplingTest: updateExportFields FAIL - " + "AvgSfcSalinity outside relative tolerance of {} for {} " + "cells", + RTol, SalinErr); + } + + // reset model clock to the start time for any subsequent tests + ModelClock->setCurrentTime(DefStepper->getStartTime()); + + SfcCoupling::clear(); + return Err; +} + +int testExportToCoupler(const CouplingLayout Layout) { + + int Err = 0; + + auto CouplingParams = mockCouplingInitParams(Layout); + Err += SfcCoupling::init(CouplingParams); + + SfcCoupling *DefCoupling = SfcCoupling::getDefault(); + OceanState *DefState = OceanState::getDefault(); + VertCoord *DefVertCoord = VertCoord::getDefault(); + + int NCells = DefCoupling->NCellsOwned; + int NImports = DefCoupling->NImportFields; + int NExports = DefCoupling->NExportFields; + + std::vector CplToOcnData(NCells * NImports, 0.0); + std::vector OcnToCplData(NCells * NExports, 0.0); + + int TempIdx = CouplingParams.ExportIdxMap.at("So_t"); + int SalinIdx = CouplingParams.ExportIdxMap.at("So_s"); + int SshIdx = CouplingParams.ExportIdxMap.at("So_ssh"); + + DefCoupling->attachData(CplToOcnData.data(), OcnToCplData.data()); + + HostArray1DReal ExpectedTemp = + makeCellVarryingArray("ExpectedTemp", NCells, Real(TempIdx)); + HostArray1DReal ExpectedSalin = + makeCellVarryingArray("ExpectedSalin", NCells, Real(SalinIdx)); + HostArray1DReal ExpectedSsh = + makeCellVarryingArray("ExpectedSsh", NCells, Real(SshIdx)); + + // Seed the averages through the public API: at NAccumSteps == 0, + // Welford's update reduces to AvgSfc* := mocked value, so one call sets + // the averages to exactly ExpectedTemp/ExpectedSalin. + I4 TempTracerIdx, SalinTracerIdx; + Tracers::getIndex(TempTracerIdx, "Temperature"); + Tracers::getIndex(SalinTracerIdx, "Salinity"); + + HostArray2DReal TempH = Tracers::getHostByIndex(0, TempTracerIdx); + HostArray2DReal SalinH = Tracers::getHostByIndex(0, SalinTracerIdx); + for (int Cell = 0; Cell < NCells; Cell++) { + int KSfc = DefVertCoord->MinLayerCellH(Cell); + TempH(Cell, KSfc) = ExpectedTemp(Cell); + SalinH(Cell, KSfc) = ExpectedSalin(Cell); + } + Tracers::copyToDevice(0); + + DefCoupling->updateExportFields(DefState, Tracers::getAll(0)); + + auto SshCellOwned = + Kokkos::subview(DefVertCoord->SshCell, std::pair(0, NCells)); + deepCopy(SshCellOwned, ExpectedSsh); + + DefCoupling->exportToCoupler(); + + // Check 1: exportToCoupler properly packs into OcnToCplView. Velocity + // is skipped here: its averaging is a hardcoded stub pending real vector + // reconstruction (see OcnToCplFields::updateAverages), not yet + // meaningful to check. + // copyToHost() converts temp to Kelvin (identity CT->PT w/ ConstantEos) + int PackErr = 0; + for (int Cell = 0; Cell < NCells; Cell++) { + if (OcnToCplData[flatIdx(Layout, Cell, TempIdx, NCells, NExports)] != + ExpectedTemp(Cell) + TkFrz) { + PackErr++; + } + if (OcnToCplData[flatIdx(Layout, Cell, SalinIdx, NCells, NExports)] != + ExpectedSalin(Cell)) { + PackErr++; + } + if (OcnToCplData[flatIdx(Layout, Cell, SshIdx, NCells, NExports)] != + ExpectedSsh(Cell)) { + PackErr++; + } + } + + if (PackErr == 0) { + LOG_INFO("SfcCouplingTest: exportToCoupler with {} layout PASS", + toString(Layout)); + } else { + Err += PackErr; + LOG_ERROR("SfcCouplingTest: exportToCoupler with {} layout FAIL - " + "{} packing errors", + toString(Layout), PackErr); + } + + // Check 2: resetFields() zeroed the running-average accumulators + HostArray1DReal ZeroArray("ZeroArray", NCells); + deepCopy(ZeroArray, 0.0); + + // AvgSfcTemperatureH holds Kelvin, so a reset (0 degC) mirrors to TkFrz + HostArray1DReal ZeroTempArray("ZeroTempArray", NCells); + deepCopy(ZeroTempArray, TkFrz); + + // Refresh OcnToCpl's own host mirrors post-reset, rather than creating + // separate test-only mirrors + DefCoupling->OcnToCpl.copyToHost(); + + int ResetErr = 0; + if (!arraysEqual(DefCoupling->OcnToCpl.AvgSfcTemperatureH, ZeroTempArray)) + ResetErr++; + if (!arraysEqual(DefCoupling->OcnToCpl.AvgSfcSalinityH, ZeroArray)) + ResetErr++; + if (!arraysEqual(DefCoupling->OcnToCpl.AvgSfcVelocityZonalH, ZeroArray)) + ResetErr++; + if (!arraysEqual(DefCoupling->OcnToCpl.AvgSfcVelocityMeridH, ZeroArray)) + ResetErr++; + + if (ResetErr == 0) { + LOG_INFO("SfcCouplingTest: exportToCoupler resetFields PASS"); + } else { + Err += ResetErr; + LOG_ERROR("SfcCouplingTest: exportToCoupler resetFields FAIL"); + } + + // Check 3: NAccumSteps was reset to 0 + if (DefCoupling->getNAccumSteps() != 0) { + Err++; + LOG_ERROR("SfcCouplingTest: exportToCoupler reset counter FAIL"); + } else { + LOG_INFO("SfcCouplingTest: exportToCoupler reset counter PASS"); + } + + SfcCoupling::clear(); + + return Err; +} + +int testEraseAndGet() { + + int Err = 0; + + TestSetup Setup; + auto *DefMesh = HorzMesh::getDefault(); + auto *DefStepper = TimeStepper::getDefault(); + TimeInterval TimeStep = DefStepper->getTimeStep(); + + // test creation of a non-default, named surface coupling object + SfcCoupling::create("AnotherSfcCoupling", DefMesh, 10, 12, + Setup.ImportIdxMap, Setup.ExportIdxMap, DefStepper, + TimeStep, CouplingLayout::MCT); + + if (SfcCoupling::get("AnotherSfcCoupling")) { + LOG_INFO("SfcCouplingTest: Non-default SfcCoupling retrieval PASS"); + } else { + Err++; + LOG_ERROR("SfcCouplingTest: Non-default SfcCoupling retrieval FAIL"); + } + + // test erase + SfcCoupling::erase("AnotherSfcCoupling"); + + // get() logs an expected error here; wrap with spdlog::set_level to + // silence it if the message becomes confusing in test logs + if (SfcCoupling::get("AnotherSfcCoupling")) { + Err++; + LOG_ERROR("SfcCouplingTest: erased SfcCoupling retrieval FAIL"); + } else { + LOG_INFO("SfcCouplingTest: erased SfcCoupling retrieval PASS"); + } + + return Err; +} + +void finalizeSfcCouplingTest() { + + Tracers::clear(); + Forcing::clear(); + OceanState::clear(); + VertCoord::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + TimeStepper::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); +} + +int sfcCouplingTest(const std::string &MeshFile = "OmegaMesh.nc") { + + int Err = initSfcCouplingTest(MeshFile); + + if (Err != 0) { + LOG_CRITICAL("SfcCouplingTest: Error initializing"); + } + + Err += testImportFromCoupler(CouplingLayout::MCT); + Err += testImportFromCoupler(CouplingLayout::MOAB); + + Err += testApplyImportFields(); + + Err += testUpdateExportFields(1); + Err += testUpdateExportFields(5); + + Err += testExportToCoupler(CouplingLayout::MCT); + Err += testExportToCoupler(CouplingLayout::MOAB); + + Err += testEraseAndGet(); + + if (Err == 0) { + LOG_INFO("SfcCouplingTest: Successful completion"); + } + + finalizeSfcCouplingTest(); + + return Err; +} + +int main(int argc, char *argv[]) { + + int RetVal = 0; + + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + RetVal += sfcCouplingTest(); + + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + if (RetVal >= 256) + RetVal = 255; + if (RetVal == 0) + LOG_INFO("------ SufaceCoupling unit tests successful ------"); + + return RetVal; +} // end of main