From 6b9c8cd3f5c05429cef540f3ebbfd8b56f322d41 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Wed, 27 May 2026 14:48:01 -0700 Subject: [PATCH 01/11] added Frazil class infrastructure --- components/omega/configs/Default.yml | 4 + components/omega/src/ocn/Frazil.cpp | 144 +++++++++++++++++++++++++ components/omega/src/ocn/Frazil.h | 84 +++++++++++++++ components/omega/src/ocn/OceanInit.cpp | 2 + 4 files changed, 234 insertions(+) create mode 100644 components/omega/src/ocn/Frazil.cpp create mode 100644 components/omega/src/ocn/Frazil.h diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 9db42e178466..0ae0ecbce9cc 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -76,6 +76,10 @@ Omega: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + Frazil: + FrazilType: teos + massLimit: 0.1 + phi: 0.75 VertMix: Background: Diffusivity: 1.0e-5 diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp new file mode 100644 index 000000000000..070eea23bc20 --- /dev/null +++ b/components/omega/src/ocn/Frazil.cpp @@ -0,0 +1,144 @@ +//===-- ocn/Frazil.cpp - Frazil Ice Formation -----------------*- C++ -*-===// +// +// The Frazil class manages frazil tendencies and accumulators. +// This initial scaffold wires allocation and configuration only. +// +//===----------------------------------------------------------------------===// + +#include "Frazil.h" +#include "Error.h" +#include "Logging.h" + +namespace OMEGA { + +Frazil *Frazil::DefaultFrazil = nullptr; +std::map> Frazil::AllFrazil; + +void Frazil::init() { + + if (!HorzMesh::getDefault() or !VertCoord::getDefault()) { + ABORT_ERROR("Frazil::init: HorzMesh and VertCoord must be initialized"); + } + + if (!DefaultFrazil) { + DefaultFrazil = create("Default"); + } +} + +Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) + : frazilChoice(FrazilType::TeosFrazil), massLimit(0.1_Real), phi(0.75_Real), + NCellsAll(Mesh->NCellsAll), + NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength) { + + FrazilTTend = + Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilSTend = + Array2DReal("FrazilSTend", Mesh->NCellsSize, VCoord->NVertLayers); + FrazilHTend = + Array2DReal("FrazilHTend", Mesh->NCellsSize, VCoord->NVertLayers); + + AccMIce = Array1DReal("AccMIce", Mesh->NCellsSize); + AccEIce = Array1DReal("AccEIce", Mesh->NCellsSize); + AccMLiq = Array1DReal("AccMLiq", Mesh->NCellsSize); + AccELiq = Array1DReal("AccELiq", Mesh->NCellsSize); + AccMSalt = Array1DReal("AccMSalt", Mesh->NCellsSize); + + deepCopy(FrazilTTend, 0.0_Real); + deepCopy(FrazilSTend, 0.0_Real); + deepCopy(FrazilHTend, 0.0_Real); + deepCopy(AccMIce, 0.0_Real); + deepCopy(AccEIce, 0.0_Real); + deepCopy(AccMLiq, 0.0_Real); + deepCopy(AccELiq, 0.0_Real); + deepCopy(AccMSalt, 0.0_Real); +} + +Frazil::~Frazil() {} + +Frazil *Frazil::create(const std::string &Name) { + if (AllFrazil.find(Name) != AllFrazil.end()) { + LOG_ERROR("Attempted to create Frazil {} but it already exists", Name); + return nullptr; + } + + auto *NewFrazil = + new Frazil(HorzMesh::getDefault(), VertCoord::getDefault()); + AllFrazil.emplace(Name, NewFrazil); + + Error Err; + Config *OmegaConfig = Config::getOmegaConfig(); + Config FrazilConfig("Frazil"); + Err += OmegaConfig->get(FrazilConfig); + CHECK_ERROR_ABORT(Err, "Frazil::create: Frazil group not found in Config"); + + std::string FrazilTypeStr; + Err += FrazilConfig.get("FrazilType", FrazilTypeStr); + CHECK_ERROR_ABORT(Err, + "Frazil::create: FrazilType not found in Frazil config"); + + if ((FrazilTypeStr == "Basic") or (FrazilTypeStr == "basic") or + (FrazilTypeStr == "BasicFrazil")) { + NewFrazil->frazilChoice = FrazilType::BasicFrazil; + } else if ((FrazilTypeStr == "Simple") or (FrazilTypeStr == "simple") or + (FrazilTypeStr == "SimpleFrazil")) { + NewFrazil->frazilChoice = FrazilType::SimpleFrazil; + } else if ((FrazilTypeStr == "teos") or (FrazilTypeStr == "Teos") or + (FrazilTypeStr == "TEOS") or (FrazilTypeStr == "Teos10") or + (FrazilTypeStr == "teos10") or (FrazilTypeStr == "TEOS10")) { + NewFrazil->frazilChoice = FrazilType::TeosFrazil; + } else { + ABORT_ERROR("Frazil::create: Unknown FrazilType requested"); + } + + Err += FrazilConfig.get("massLimit", NewFrazil->massLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: massLimit not found in Frazil config"); + + Err += FrazilConfig.get("phi", NewFrazil->phi); + CHECK_ERROR_ABORT(Err, "Frazil::create: phi not found in Frazil config"); + + if (Name == "Default") { + DefaultFrazil = NewFrazil; + } + + return NewFrazil; +} + +Frazil *Frazil::getDefault() { return DefaultFrazil; } + +Frazil *Frazil::get(const std::string &Name) { + auto it = AllFrazil.find(Name); + if (it != AllFrazil.end()) { + return it->second.get(); + } + + LOG_ERROR("Frazil::get: Attempted to retrieve non-existent Frazil {}", Name); + return nullptr; +} + +void Frazil::erase(std::string InName) { + auto *ToErase = get(InName); + AllFrazil.erase(InName); + if (ToErase == DefaultFrazil) { + DefaultFrazil = nullptr; + } +} + +void Frazil::clear() { + AllFrazil.clear(); + DefaultFrazil = nullptr; +} + +void Frazil::checkFrazil() { + // Placeholder for frazil consistency checks. +} + +void Frazil::computeFrazilFormation() { + // Placeholder for frazil formation implementation. +} + +void Frazil::computeFrazilMelt() { + // Placeholder for frazil melt implementation. +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h new file mode 100644 index 000000000000..d739fc9e1546 --- /dev/null +++ b/components/omega/src/ocn/Frazil.h @@ -0,0 +1,84 @@ +#ifndef OMEGA_FRAZIL_H +#define OMEGA_FRAZIL_H +//===-- ocn/Frazil.h - Frazil Ice Formation -------------------*- C++ -*-===// +// +// This header defines a scaffold for frazil-related tendencies and +// accumulators. Physics implementations are intentionally left empty. +// +//===----------------------------------------------------------------------===// + +#include "Config.h" +#include "DataTypes.h" +#include "HorzMesh.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include +#include +#include + +namespace OMEGA { + +enum class FrazilType { + BasicFrazil, ///< Placeholder basic frazil option + SimpleFrazil, ///< Placeholder simple frazil option + TeosFrazil ///< Placeholder TEOS frazil option +}; + +class Frazil { + public: + static void init(); + /// Creates a new frazil object and stores it in the AllFrazil map. + static Frazil *create(const std::string &Name); + + /// Retrieve frazil object by name. + static Frazil *get(const std::string &Name); + + /// Retrieve default frazil object. + static Frazil *getDefault(); + + /// Destructor + ~Frazil(); + + /// Deallocates arrays + static void clear(); + + /// Remove frazil object by name. + static void erase(std::string InName); ///< [in] name to remove + + Array2DReal FrazilTTend; + Array2DReal FrazilSTend; + Array2DReal FrazilHTend; + Array1DReal AccMIce; + Array1DReal AccEIce; + Array1DReal AccMLiq; + Array1DReal AccELiq; + Array1DReal AccMSalt; + + void checkFrazil(); + + private: + static Frazil *DefaultFrazil; + static std::map> AllFrazil; + + Frazil(const HorzMesh *Mesh, const VertCoord *VCoord); + + // Forbid copy and move construction/assignment. + Frazil(const Frazil &) = delete; + Frazil &operator=(const Frazil &) = delete; + Frazil(Frazil &&) = delete; + Frazil &operator=(Frazil &&) = delete; + + FrazilType frazilChoice; + Real massLimit; + Real phi; + I4 NCellsAll; + I4 NChunks; + + void computeFrazilFormation(); + void computeFrazilMelt(); +}; + +} // namespace OMEGA + +#endif diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index ceabc0ede033..64f975570b10 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -15,6 +15,7 @@ #include "Error.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -199,6 +200,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { Forcing::init(); AuxiliaryState::init(); Eos::init(); + Frazil::init(); PressureGrad::init(); Tendencies::init(); From eb972f6a39f0a0a8efd70d5f3c9187ab15aa1bde Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 28 May 2026 14:35:06 -0700 Subject: [PATCH 02/11] first implementation of FrazilFormation --- components/omega/src/ocn/Frazil.cpp | 82 +++++++++++++++++++++---- components/omega/src/ocn/Frazil.h | 94 +++++++++++++++++++++++++++-- 2 files changed, 159 insertions(+), 17 deletions(-) diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 070eea23bc20..5fde6eca5669 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -14,6 +14,12 @@ namespace OMEGA { Frazil *Frazil::DefaultFrazil = nullptr; std::map> Frazil::AllFrazil; +/// Constructor for FrazilFormation +FrazilFormation::FrazilFormation() {} + +/// Constructor for FrazilMelt +FrazilMelt::FrazilMelt() {} + void Frazil::init() { if (!HorzMesh::getDefault() or !VertCoord::getDefault()) { @@ -28,7 +34,8 @@ void Frazil::init() { Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) : frazilChoice(FrazilType::TeosFrazil), massLimit(0.1_Real), phi(0.75_Real), NCellsAll(Mesh->NCellsAll), - NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength) { + NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh), + VCoordPtr(VCoord), computeFrazilFormation(), computeFrazilMelt() { FrazilTTend = Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers); @@ -90,11 +97,12 @@ Frazil *Frazil::create(const std::string &Name) { ABORT_ERROR("Frazil::create: Unknown FrazilType requested"); } - Err += FrazilConfig.get("massLimit", NewFrazil->massLimit); + Err += FrazilConfig.get("massLimit", + NewFrazil->computeFrazilFormation.MassLimit); CHECK_ERROR_ABORT(Err, "Frazil::create: massLimit not found in Frazil config"); - Err += FrazilConfig.get("phi", NewFrazil->phi); + Err += FrazilConfig.get("phi", NewFrazil->computeFrazilFormation.Phi); CHECK_ERROR_ABORT(Err, "Frazil::create: phi not found in Frazil config"); if (Name == "Default") { @@ -129,16 +137,64 @@ void Frazil::clear() { DefaultFrazil = nullptr; } -void Frazil::checkFrazil() { - // Placeholder for frazil consistency checks. -} - -void Frazil::computeFrazilFormation() { - // Placeholder for frazil formation implementation. -} - -void Frazil::computeFrazilMelt() { - // Placeholder for frazil melt implementation. +void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &LayerH) { + OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell); + + OMEGA_SCOPE(LocComputeFrazilFormation, computeFrazilFormation); + OMEGA_SCOPE(LocComputeFrazilMelt, computeFrazilMelt); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocAccMIce, AccMIce); + OMEGA_SCOPE(LocAccEIce, AccEIce); + OMEGA_SCOPE(LocAccMLiq, AccMLiq); + OMEGA_SCOPE(LocAccELiq, AccELiq); + OMEGA_SCOPE(LocAccMSalt, AccMSalt); + + parallelFor( + {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + // Explicit accumulation order: bottom layer to top layer. + for (I4 K = KMax; K >= KMin; --K) { + const Real SAIn = SA(ICell, K); + const Real CTIn = CT(ICell, K); + const Real PIn = P(ICell, K); + const Real H = LayerH(ICell, K); + + const Real Tfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + Real Dt = 1800.0_Real; // hard-coded for now (30min in s) + + if (CTIn < Tfrz) { + LocComputeFrazilFormation( + SAIn, CTIn, PIn, H, Dt, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), + LocAccEIce(ICell), HTend, TTend, STend); + } else { + LocComputeFrazilMelt(LocAccMIce(ICell), LocAccMLiq(ICell), + LocAccMSalt(ICell), LocAccELiq(ICell), + LocAccEIce(ICell), HTend, TTend, STend); + } + + // LocAccMIce(ICell) += solidMass; + // LocAccMLiq(ICell) += liquidMass; + // LocAccMSalt(ICell) += liquidMass * SAnew; + // LocAccELiq(ICell) += liquidMass * Cp0Sw * CTnew; + // LocAccEIce(ICell) += + // solidMass * gsw_pot_enthalpy_from_pt_ice_poly(CTnew); + + LocFrazilHTend(ICell, K) = HTend; + LocFrazilTTend(ICell, K) = TTend; + LocFrazilSTend(ICell, K) = STend; + } + }); } } // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index d739fc9e1546..732bfda52e24 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -8,7 +8,8 @@ //===----------------------------------------------------------------------===// #include "Config.h" -#include "DataTypes.h" +// #include "DataTypes.h" +#include "GlobalConstants.h" #include "HorzMesh.h" #include "OmegaKokkos.h" #include "VertCoord.h" @@ -17,6 +18,8 @@ #include #include +#include + namespace OMEGA { enum class FrazilType { @@ -25,6 +28,85 @@ enum class FrazilType { TeosFrazil ///< Placeholder TEOS frazil option }; +class FrazilMelt { + public: + /// constructor declaration + FrazilMelt(); + + // The functor takes the full arrays of specific volume (inout), + // the indices ICell and KChunk, and the ocean tracers (conservative) + // temperature, and (absolute) salinity as inputs, and outputs the + // specific volume according to the Roquet et al. 2015 75 term expansion. + KOKKOS_FUNCTION void operator()(Real &AccMIce, Real &AccMLiq, Real &AccMSalt, + Real &AccELiq, Real &AccEIce, Real &HTend, + Real &TTend, Real &STend) const { + (void)AccMIce; + (void)AccMLiq; + (void)AccMSalt; + (void)AccELiq; + (void)AccEIce; + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + } +}; + +class FrazilFormation { + public: + /// constructor declaration + FrazilFormation(); + + /// Parameters for FrazilFormation (overwritten by config file if set there) + Real Phi = 0.75_Real; ///< liquid mass fraction of frazil for export + Real MassLimit = + 0.1_Real; ///< layer mass fraction limit for thickness tendency + + // The functor takes the full arrays of specific volume (inout), + // the indices ICell and KChunk, and the ocean tracers (conservative) + // temperature, and (absolute) salinity as inputs, and outputs the + // specific volume according to the Roquet et al. 2015 75 term expansion. + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P, + const Real h, const Real Dt, Real &AccMIce, + Real &AccMLiq, Real &AccMSalt, Real &AccELiq, + Real &AccEIce, Real &HTend, Real &TTend, + Real &STend) const { + + Real CTnew; + Real SAnew; + Real wIh = 0.0_Real; + Real solidMass = 0.0_Real; + Real liquidMass = 0.0_Real; + + // type casting for now because gsw expects double + // and Real can be either float or double depending on build configuration + double SAnew_d = 0.0; + double CTnew_d = 0.0; + double wIh_d = 0.0; + + gsw_frazil_properties_potential_poly( + static_cast(SA), static_cast(Cp0Sw * CT), + static_cast(P), &SAnew_d, &CTnew_d, &wIh_d); + + SAnew = static_cast(SAnew_d); + CTnew = static_cast(CTnew_d); + wIh = static_cast(wIh_d); + + const Real OneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - Phi); + solidMass = h * Kokkos::min(wIh, OneMinusPhi * MassLimit); + liquidMass = (Phi / OneMinusPhi) * solidMass; + + HTend = -(solidMass + liquidMass) / Dt; + TTend = (h - solidMass - liquidMass) * CTnew - h * CT; + STend = (h - solidMass - liquidMass) * SAnew - h * SA; + + AccMIce += solidMass; + AccMLiq += liquidMass; + AccMSalt += liquidMass * SAnew; + AccELiq += liquidMass * Cp0Sw * CTnew; + AccEIce += solidMass * gsw_pot_enthalpy_from_pt_ice_poly(CTnew); + } +}; + class Frazil { public: static void init(); @@ -55,7 +137,8 @@ class Frazil { Array1DReal AccELiq; Array1DReal AccMSalt; - void checkFrazil(); + void computeFrazil(const Array2DReal &CT, const Array2DReal &SA, + const Array2DReal &P, const Array2DReal &H); private: static Frazil *DefaultFrazil; @@ -70,13 +153,16 @@ class Frazil { Frazil &operator=(Frazil &&) = delete; FrazilType frazilChoice; + FrazilFormation computeFrazilFormation; + FrazilMelt computeFrazilMelt; Real massLimit; Real phi; + // Real dt; I4 NCellsAll; I4 NChunks; - void computeFrazilFormation(); - void computeFrazilMelt(); + const HorzMesh *MeshPtr; + const VertCoord *VCoordPtr; }; } // namespace OMEGA From 8302bbf9523ee40c9e99a8b00b63eeaa19522243 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 28 May 2026 15:22:28 -0700 Subject: [PATCH 03/11] added Frazil unit tests on Formation --- components/omega/src/ocn/Frazil.h | 5 +- components/omega/test/CMakeLists.txt | 7 + components/omega/test/ocn/FrazilTest.cpp | 216 +++++++++++++++++++++++ 3 files changed, 227 insertions(+), 1 deletion(-) create mode 100644 components/omega/test/ocn/FrazilTest.cpp diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index 732bfda52e24..eba2d1912eb7 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -99,9 +99,12 @@ class FrazilFormation { TTend = (h - solidMass - liquidMass) * CTnew - h * CT; STend = (h - solidMass - liquidMass) * SAnew - h * SA; + // these are not currently mass or energy, they all need a RhoSw factor AccMIce += solidMass; AccMLiq += liquidMass; - AccMSalt += liquidMass * SAnew; + AccMSalt += + liquidMass * SAnew * PPt2Salt; // the PPt2Salt scaling could be moved + // to the coupler interaction AccELiq += liquidMass * Cp0Sw * CTnew; AccEIce += solidMass * gsw_pot_enthalpy_from_pt_ice_poly(CTnew); } diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 50c42bf0a2f0..a7c949da1f19 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -512,6 +512,13 @@ add_omega_test( "-n;2" ) +add_omega_test( + FRAZIL_TEST + testFrazil.exe + ocn/FrazilTest.cpp + "-n;2" +) + ################ # VertCoord test ################ diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp new file mode 100644 index 000000000000..8c5f882c53bc --- /dev/null +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -0,0 +1,216 @@ +//===-- Test driver for OMEGA FrazilFormation ---------------------------*- C++ +//-*-===// +// +/// \file +/// \brief Minimal test driver for OMEGA frazil formation functor +// +//===-----------------------------------------------------------------------===/ + +#include "Frazil.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Field.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 "VertCoord.h" +#include "mpi.h" + +using namespace OMEGA; + +constexpr int NVertLayers = 60; + +void initFrazilTest(const std::string &mesh) { + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + LOG_INFO("------ Frazil Unit Tests ------"); + + Config("Omega"); + Config::readAll("omega.yml"); + + IO::init(DefComm); + IOStream::init(); + Decomp::init(mesh); + Halo::init(); + HorzMesh::init(); + VertCoord::init(false); + Frazil::init(); +} + +void finalizeFrazilTest() { + Frazil::clear(); + VertCoord::clear(); + HorzMesh::clear(); + Halo::clear(); + Decomp::clear(); + Field::clear(); + Dimension::clear(); + MachEnv::removeAll(); +} + +void testFrazilFormationCold() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = -2.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real Dt = 1800.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + FrazilFormation ComputeFrazilFormation; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, Dt, AccMIce, AccMLiq, AccMSalt, + AccELiq, AccEIce, HTend, TTend, STend); + + if (AccMIce <= 0.0_Real) { + ABORT_ERROR("FrazilTestCold: accumulated ice mass is non-positive: {}", + AccMIce); + } + if (AccMLiq <= 0.0_Real) { + ABORT_ERROR("FrazilTestCold: accumulated liquid mass is non-positive: {}", + AccMLiq); + } + if (AccMSalt <= 0.0_Real) { + ABORT_ERROR("FrazilTestCold: accumulated salt mass is non-positive: {}", + AccMSalt); + } + + if (isApprox(AccELiq, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilTestCold: accumulated liquid energy is effectively zero: {}", + AccELiq); + } + if (isApprox(AccEIce, 0.0_Real, RTol)) { + ABORT_ERROR( + "FrazilTestCold: accumulated ice energy is effectively zero: {}", + AccEIce); + } + if (isApprox(HTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTestCold: HTend is effectively zero: {}", HTend); + } + + if (isApprox(TTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTestCold: TTend is zero: {}", TTend); + } + + if (isApprox(STend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTestCold: STend is effectively zero: {}", STend); + } + LOG_INFO("FrazilTestCold: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +void testFrazilFormationWarm() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + VCoord->NVertLayers = NVertLayers; + + const Real SAIn = 35.0_Real; + const Real CTIn = 10.0_Real; + const Real PIn = 100.0_Real; + const Real h = 10.0_Real; + const Real Dt = 1800.0_Real; + const Real RTol = 1e-10_Real; + + (void)Mesh; + + FrazilFormation ComputeFrazilFormation; + + Real AccMIce = 0.0_Real; + Real AccMLiq = 0.0_Real; + Real AccMSalt = 0.0_Real; + Real AccELiq = 0.0_Real; + Real AccEIce = 0.0_Real; + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + ComputeFrazilFormation(SAIn, CTIn, PIn, h, Dt, AccMIce, AccMLiq, AccMSalt, + AccELiq, AccEIce, HTend, TTend, STend); + + if (!isApprox(AccMIce, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero AccMIce, got {}", AccMIce); + } + + if (!isApprox(AccMSalt, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero AccMSalt, got {}", AccMSalt); + } + if (!isApprox(AccMLiq, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero AccMLiq, got {}", AccMLiq); + } + + if (!isApprox(AccELiq, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero AccELiq, got {}", AccELiq); + } + if (!isApprox(AccEIce, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero AccEIce, got {}", AccEIce); + } + if (!isApprox(HTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero HTend, got {}", HTend); + } + + if (!isApprox(TTend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero TTend, got {}", TTend); + } + + if (!isApprox(STend, 0.0_Real, RTol)) { + ABORT_ERROR("FrazilTest warm: expected zero STend, got {}", STend); + } + LOG_INFO("FrazilTestWarm: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); +} + +void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { + initFrazilTest(MeshFile); + testFrazilFormationCold(); + testFrazilFormationWarm(); + finalizeFrazilTest(); +} + +int main(int argc, char *argv[]) { + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + frazilTest(); + + LOG_INFO("------ Frazil Unit Tests Successful ------"); + + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + return 0; +} From 49488a9fee223e30fe16d2753e68256dbe5cfd7a Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Fri, 29 May 2026 11:49:28 -0700 Subject: [PATCH 04/11] added a multi-layer test (incl kernel logging) and FrazilMelt --- components/omega/src/ocn/Frazil.cpp | 36 +++--- components/omega/src/ocn/Frazil.h | 138 ++++++++++++++++++----- components/omega/test/ocn/FrazilTest.cpp | 106 ++++++++++++++++- 3 files changed, 230 insertions(+), 50 deletions(-) diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 5fde6eca5669..7b253b8d89ce 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -170,30 +170,36 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, Real HTend = 0.0_Real; Real TTend = 0.0_Real; Real STend = 0.0_Real; - Real Dt = 1800.0_Real; // hard-coded for now (30min in s) if (CTIn < Tfrz) { - LocComputeFrazilFormation( - SAIn, CTIn, PIn, H, Dt, LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), - LocAccEIce(ICell), HTend, TTend, STend); + LocComputeFrazilFormation(SAIn, CTIn, PIn, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), + HTend, TTend, STend); } else { - LocComputeFrazilMelt(LocAccMIce(ICell), LocAccMLiq(ICell), - LocAccMSalt(ICell), LocAccELiq(ICell), - LocAccEIce(ICell), HTend, TTend, STend); + LocComputeFrazilMelt(SAIn, CTIn, PIn, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), + HTend, TTend, STend); } - // LocAccMIce(ICell) += solidMass; - // LocAccMLiq(ICell) += liquidMass; - // LocAccMSalt(ICell) += liquidMass * SAnew; - // LocAccELiq(ICell) += liquidMass * Cp0Sw * CTnew; - // LocAccEIce(ICell) += - // solidMass * gsw_pot_enthalpy_from_pt_ice_poly(CTnew); + // temporary log -- TBRemoved + LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " + "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", + ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), + LocAccEIce(ICell)); - LocFrazilHTend(ICell, K) = HTend; + LocFrazilHTend(ICell, K) = HTend; // not scaled by dt LocFrazilTTend(ICell, K) = TTend; LocFrazilSTend(ICell, K) = STend; } + // Convert to coupler units + LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; + LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; + LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; + LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; + LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; }); } diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index eba2d1912eb7..87ae523aebca 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -37,17 +37,84 @@ class FrazilMelt { // the indices ICell and KChunk, and the ocean tracers (conservative) // temperature, and (absolute) salinity as inputs, and outputs the // specific volume according to the Roquet et al. 2015 75 term expansion. - KOKKOS_FUNCTION void operator()(Real &AccMIce, Real &AccMLiq, Real &AccMSalt, - Real &AccELiq, Real &AccEIce, Real &HTend, - Real &TTend, Real &STend) const { - (void)AccMIce; - (void)AccMLiq; - (void)AccMSalt; - (void)AccELiq; - (void)AccEIce; - HTend = 0.0_Real; - TTend = 0.0_Real; - STend = 0.0_Real; + KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P, + const Real h, Real &AccMIce, Real &AccMLiq, + Real &AccMSalt, Real &AccELiq, Real &AccEIce, + Real &HTend, Real &TTend, + Real &STend) const { + + constexpr Real Eps = 1.0e-12_Real; + + if (AccMIce <= Eps) { // potential leak if we dont redistribute + AccMIce = 0.0_Real; + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + return; + } + + const Real safeAccMLiq = Kokkos::max(AccMLiq, Eps); + const Real safeAccMIce = Kokkos::max(AccMIce, Eps); + const Real frazilIceFraction = Kokkos::min( + 1.0_Real, + Kokkos::max(0.0_Real, AccMIce / (safeAccMLiq + safeAccMIce))); + const Real brineSalinity = AccMSalt / safeAccMLiq; + const Real brineEnthalpy = AccELiq / safeAccMLiq; + const Real potEnthalpyIce = AccEIce / safeAccMIce; + + const Real layerMass = h + AccMIce; + const Real safeLayerMass = Kokkos::max(layerMass, Eps); + + const Real layerIceFraction = + Kokkos::min(1.0_Real, Kokkos::max(0.0_Real, AccMIce / safeLayerMass)); + + // typecasting for now but will be simplified + const double SA_d = static_cast(SA); + const double CT_d = static_cast(CT); + const double P_d = static_cast(P); + const double wIhIn_d = static_cast(layerIceFraction); + const double pt0Ice_d = + gsw_pt_from_pot_enthalpy_ice(static_cast(potEnthalpyIce)); + const double tIce_d = gsw_t_from_pt0_ice(pt0Ice_d, P_d); + double SAnew_d = SA_d; + double CTnew_d = CT_d; + double wIhOut_d = wIhIn_d; + + gsw_melting_ice_into_seawater(SA_d, CT_d, P_d, wIhIn_d, tIce_d, &SAnew_d, + &CTnew_d, &wIhOut_d); + + const Real wIhOut = Kokkos::min( + 1.0_Real, Kokkos::max(0.0_Real, static_cast(wIhOut_d))); + + const Real finalSolidMass = + Kokkos::min(AccMIce, Kokkos::max(0.0_Real, wIhOut * safeLayerMass)); + const Real solidMass = Kokkos::max(0.0_Real, AccMIce - finalSolidMass); + + if (solidMass <= Eps) { + return; + } + + const Real liquidMass = + Kokkos::min(AccMLiq, solidMass * (1.0_Real - frazilIceFraction) / + Kokkos::max(frazilIceFraction, Eps)); + const Real solidEnthalpy = solidMass * potEnthalpyIce; + const Real liquidEnthalpy = liquidMass * brineEnthalpy; + + HTend = +(solidMass + liquidMass); + TTend = +(liquidEnthalpy + solidEnthalpy) / Cp0Sw; + STend = +(liquidMass * brineSalinity); + + AccMIce = Kokkos::max(0.0_Real, AccMIce - solidMass); + AccMLiq = Kokkos::max(0.0_Real, AccMLiq - liquidMass); + AccMSalt = Kokkos::max(0.0_Real, AccMSalt - liquidMass * brineSalinity); + AccELiq -= liquidEnthalpy; + AccEIce -= solidEnthalpy; + if (AccMIce <= Eps) { + AccEIce = 0.0_Real; + } + if (AccMLiq <= Eps) { + AccELiq = 0.0_Real; + } } }; @@ -66,26 +133,31 @@ class FrazilFormation { // temperature, and (absolute) salinity as inputs, and outputs the // specific volume according to the Roquet et al. 2015 75 term expansion. KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P, - const Real h, const Real Dt, Real &AccMIce, - Real &AccMLiq, Real &AccMSalt, Real &AccELiq, - Real &AccEIce, Real &HTend, Real &TTend, + const Real h, Real &AccMIce, Real &AccMLiq, + Real &AccMSalt, Real &AccELiq, Real &AccEIce, + Real &HTend, Real &TTend, Real &STend) const { Real CTnew; Real SAnew; - Real wIh = 0.0_Real; - Real solidMass = 0.0_Real; - Real liquidMass = 0.0_Real; + Real wIh = 0.0_Real; + Real solidMass = 0.0_Real; + Real liquidMass = 0.0_Real; + Real solidEnthalpy = 0.0_Real; + Real liquidEnthalpy = 0.0_Real; // type casting for now because gsw expects double // and Real can be either float or double depending on build configuration double SAnew_d = 0.0; double CTnew_d = 0.0; double wIh_d = 0.0; + // double PTnew_d = 0.0; gsw_frazil_properties_potential_poly( static_cast(SA), static_cast(Cp0Sw * CT), static_cast(P), &SAnew_d, &CTnew_d, &wIh_d); + double PTnew_d = + gsw_pt_from_ct(SAnew_d, CTnew_d); // convert to potential temperature SAnew = static_cast(SAnew_d); CTnew = static_cast(CTnew_d); @@ -94,19 +166,24 @@ class FrazilFormation { const Real OneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - Phi); solidMass = h * Kokkos::min(wIh, OneMinusPhi * MassLimit); liquidMass = (Phi / OneMinusPhi) * solidMass; - - HTend = -(solidMass + liquidMass) / Dt; - TTend = (h - solidMass - liquidMass) * CTnew - h * CT; - STend = (h - solidMass - liquidMass) * SAnew - h * SA; - - // these are not currently mass or energy, they all need a RhoSw factor - AccMIce += solidMass; - AccMLiq += liquidMass; - AccMSalt += - liquidMass * SAnew * PPt2Salt; // the PPt2Salt scaling could be moved - // to the coupler interaction - AccELiq += liquidMass * Cp0Sw * CTnew; - AccEIce += solidMass * gsw_pot_enthalpy_from_pt_ice_poly(CTnew); + solidEnthalpy = solidMass * gsw_pot_enthalpy_from_pt_ice_poly(PTnew_d); + liquidEnthalpy = liquidMass * Cp0Sw * CTnew; + + // per timestep (not scaled by dt here) + HTend = -(solidMass + liquidMass); + // TTend = (h - solidMass - liquidMass) * CTnew - h * CT; + // STend = (h - solidMass - liquidMass) * SAnew - h * SA; + TTend = -(liquidEnthalpy + solidEnthalpy) / + Cp0Sw; // convert back to CT tendency + STend = -(liquidMass * SAnew); + + // Local unit of mass is pseudo thickness (m) + // these all need a RhoSw factor before coupling + AccMIce += solidMass; // m + AccMLiq += liquidMass; // m + AccMSalt += liquidMass * SAnew; // (m)(g/kg) + AccELiq += liquidEnthalpy; // (m)(J/kg) + AccEIce += solidEnthalpy; // (m)(J/kg) } }; @@ -160,7 +237,6 @@ class Frazil { FrazilMelt computeFrazilMelt; Real massLimit; Real phi; - // Real dt; I4 NCellsAll; I4 NChunks; diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp index 8c5f882c53bc..98cdab844d1b 100644 --- a/components/omega/test/ocn/FrazilTest.cpp +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -69,7 +69,6 @@ void testFrazilFormationCold() { const Real CTIn = -2.0_Real; const Real PIn = 100.0_Real; const Real h = 10.0_Real; - const Real Dt = 1800.0_Real; const Real RTol = 1e-10_Real; (void)Mesh; @@ -86,7 +85,7 @@ void testFrazilFormationCold() { Real TTend = 0.0_Real; Real STend = 0.0_Real; - ComputeFrazilFormation(SAIn, CTIn, PIn, h, Dt, AccMIce, AccMLiq, AccMSalt, + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); if (AccMIce <= 0.0_Real) { @@ -138,7 +137,6 @@ void testFrazilFormationWarm() { const Real CTIn = 10.0_Real; const Real PIn = 100.0_Real; const Real h = 10.0_Real; - const Real Dt = 1800.0_Real; const Real RTol = 1e-10_Real; (void)Mesh; @@ -155,7 +153,7 @@ void testFrazilFormationWarm() { Real TTend = 0.0_Real; Real STend = 0.0_Real; - ComputeFrazilFormation(SAIn, CTIn, PIn, h, Dt, AccMIce, AccMLiq, AccMSalt, + ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); if (!isApprox(AccMIce, 0.0_Real, RTol)) { @@ -191,10 +189,110 @@ void testFrazilFormationWarm() { AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); } +void testComputeFrazilColumn() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilTestColumn: default frazil object is null"); + } + + const Real RTol = 1e-10_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 10.0_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 4) { + ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + deepCopy(CT, CTH); + + TestFrazil->computeFrazil(CT, SA, P, H); + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + if (HTendH(ICell, KBottom0) >= 0.0_Real || + TTendH(ICell, KBottom0) <= 0.0_Real || + STendH(ICell, KBottom0) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: bottom cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KBottom1) >= 0.0_Real || + TTendH(ICell, KBottom1) <= 0.0_Real || + STendH(ICell, KBottom1) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: second cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + if (HTendH(ICell, KWarm) < 0.0_Real || TTendH(ICell, KWarm) > 0.0_Real || + STendH(ICell, KWarm) < 0.0_Real) { + ABORT_ERROR("FrazilTestColumn: warm layer sign check failed (HTend>=0, " + "TTend<=0, STend>=0 expected)"); + } + + if (HTendH(ICell, KTopCold) >= 0.0_Real || + TTendH(ICell, KTopCold) <= 0.0_Real || + STendH(ICell, KTopCold) >= 0.0_Real) { + ABORT_ERROR( + "FrazilTestColumn: top cold layer sign check failed (HTend<0, " + "TTend>0, STend<0 expected)"); + } + + LOG_INFO( + "FrazilTestColumn: XTend branch-switching checks passed for ICell={}", + ICell); +} + void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { initFrazilTest(MeshFile); testFrazilFormationCold(); testFrazilFormationWarm(); + testComputeFrazilColumn(); finalizeFrazilTest(); } From fe1b1b0c1d5e5d58cd2d74af4a8643f4c0f2ae4b Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 30 Jun 2026 09:48:44 -0700 Subject: [PATCH 05/11] added conservation check, depth limit, and tendencies --- components/omega/configs/Default.yml | 8 +- components/omega/src/ocn/Frazil.cpp | 143 ++++++++++++++++-- components/omega/src/ocn/Frazil.h | 5 + components/omega/src/ocn/OceanFinal.cpp | 3 +- components/omega/src/ocn/Tendencies.cpp | 15 +- components/omega/src/ocn/Tendencies.h | 1 + components/omega/src/ocn/TendencyTerms.cpp | 68 +++++++++ components/omega/src/ocn/TendencyTerms.h | 22 +++ components/omega/test/ocn/FrazilTest.cpp | 160 +++++++++++++++++---- 9 files changed, 385 insertions(+), 40 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 0ae0ecbce9cc..bb49ff8c6bbb 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -63,6 +63,7 @@ Omega: VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true PressureGradTendencyEnable: true + FrazilTendencyEnable: true ManufacturedSolution: WavelengthX: 5.0e6 WavelengthY: 4.33013e6 @@ -77,9 +78,12 @@ Omega: DRhoDS: 0.8 RhoT0S0: 1000.0 Frazil: + Enable: true FrazilType: teos - massLimit: 0.1 - phi: 0.75 + MassLimit: 0.1 + Phi: 0.75 + DepthLimit: -1.0 + ConservationCheck: false VertMix: Background: Diffusivity: 1.0e-5 diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 7b253b8d89ce..d2db57c8031f 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -11,6 +11,21 @@ namespace OMEGA { +namespace { + +KOKKOS_INLINE_FUNCTION bool isApprox(const Real X, const Real Y, + const Real RTol, const Real ATol = 0) { + if (Kokkos::isnan(X) || Kokkos::isnan(Y) || Kokkos::isinf(X) || + Kokkos::isinf(Y)) { + return false; + } + + return Kokkos::abs(X - Y) <= + Kokkos::max(ATol, RTol * Kokkos::max(Kokkos::abs(X), Kokkos::abs(Y))); +} + +} // namespace + Frazil *Frazil::DefaultFrazil = nullptr; std::map> Frazil::AllFrazil; @@ -97,13 +112,29 @@ Frazil *Frazil::create(const std::string &Name) { ABORT_ERROR("Frazil::create: Unknown FrazilType requested"); } - Err += FrazilConfig.get("massLimit", + Err += FrazilConfig.get("MassLimit", NewFrazil->computeFrazilFormation.MassLimit); CHECK_ERROR_ABORT(Err, - "Frazil::create: massLimit not found in Frazil config"); + "Frazil::create: MassLimit not found in Frazil config"); - Err += FrazilConfig.get("phi", NewFrazil->computeFrazilFormation.Phi); - CHECK_ERROR_ABORT(Err, "Frazil::create: phi not found in Frazil config"); + Err += FrazilConfig.get("Phi", NewFrazil->computeFrazilFormation.Phi); + CHECK_ERROR_ABORT(Err, "Frazil::create: Phi not found in Frazil config"); + + Error CheckColumnErr = + FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); + if (!CheckColumnErr.isSuccess()) { + NewFrazil->conservationCheck = false; + } + + Error EnabledErr = FrazilConfig.get("Enable", NewFrazil->Enabled); + if (!EnabledErr.isSuccess()) { + NewFrazil->Enabled = true; + } + + Error DepthLimitErr = FrazilConfig.get("DepthLimit", NewFrazil->depthLimit); + if (!DepthLimitErr.isSuccess()) { + NewFrazil->depthLimit = -1.0_Real; + } if (Name == "Default") { DefaultFrazil = NewFrazil; @@ -137,10 +168,74 @@ void Frazil::clear() { DefaultFrazil = nullptr; } +void Frazil::checkColumnConservation() const { + auto MinLayerCellH = createHostMirrorCopy(VCoordPtr->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoordPtr->MaxLayerCell); + auto FrazilHTendH = createHostMirrorCopy(FrazilHTend); + auto FrazilTTendH = createHostMirrorCopy(FrazilTTend); + auto FrazilSTendH = createHostMirrorCopy(FrazilSTend); + auto AccMIceH = createHostMirrorCopy(AccMIce); + auto AccMLiqH = createHostMirrorCopy(AccMLiq); + auto AccMSaltH = createHostMirrorCopy(AccMSalt); + auto AccELiqH = createHostMirrorCopy(AccELiq); + auto AccEIceH = createHostMirrorCopy(AccEIce); + + constexpr Real RTol = 1.0e-10_Real; + + for (I4 ICell = 0; ICell < NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + Real MassTend = 0.0_Real; + Real EnergyTend = 0.0_Real; + Real SaltTend = 0.0_Real; + + for (I4 K = KMin; K <= KMax; ++K) { + MassTend += FrazilHTendH(ICell, K); + EnergyTend += FrazilTTendH(ICell, K); + SaltTend += FrazilSTendH(ICell, K); + } + + const Real MassTotal = AccMIceH(ICell) + AccMLiqH(ICell); + const Real EnergyTotal = AccELiqH(ICell) + AccEIceH(ICell); + const Real SaltTotal = AccMSaltH(ICell); + if (ICell == 0) { + LOG_INFO("Frazil column conservation check: cell {} MassTend={} " + "MassTotal={} " + "EnergyTend={} EnergyTotal={} SaltTend={} SaltTotal={}", + ICell, MassTend * RhoSw, MassTotal, + EnergyTend * Cp0Sw * RhoSw, EnergyTotal, + SaltTend * RhoSw * PPt2Salt, SaltTotal); + LOG_INFO("Frazil column conservation check: cell {} EpsMass={} " + "EpsE={} EpsS={} ", + ICell, MassTend * RhoSw + MassTotal, + EnergyTend * Cp0Sw * RhoSw + EnergyTotal, + SaltTend * RhoSw * PPt2Salt + SaltTotal); + } + + if (!isApprox(-MassTend * RhoSw, MassTotal, RTol)) { + ABORT_ERROR( + "Frazil column mass check failed: cell {} tendency={} total={}", + ICell, -MassTend * RhoSw, MassTotal); + } + if (!isApprox(-EnergyTend * Cp0Sw * RhoSw, EnergyTotal, RTol)) { + ABORT_ERROR( + "Frazil column energy check failed: cell {} tendency={} total={}", + ICell, -EnergyTend * Cp0Sw * RhoSw, EnergyTotal); + } + if (!isApprox(-SaltTend * RhoSw * PPt2Salt, SaltTotal, RTol)) { + ABORT_ERROR( + "Frazil column salt check failed: cell {} tendency={} total={}", + ICell, -SaltTend * RhoSw * PPt2Salt, SaltTotal); + } + } +} + void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const Array2DReal &P, const Array2DReal &LayerH) { OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell); + OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid); OMEGA_SCOPE(LocComputeFrazilFormation, computeFrazilFormation); OMEGA_SCOPE(LocComputeFrazilMelt, computeFrazilMelt); @@ -158,8 +253,29 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const I4 KMin = MinLayerCell(ICell); const I4 KMax = MaxLayerCell(ICell); + I4 Klim = KMax; + bool HasKlim = true; + const bool Limit = (depthLimit >= 0.0_Real); + if (Limit) { + HasKlim = false; + for (I4 K = KMax; K >= KMin; --K) { + if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) { + Klim = K; + HasKlim = true; + break; + } + } + } + // Explicit accumulation order: bottom layer to top layer. for (I4 K = KMax; K >= KMin; --K) { + if (!HasKlim || K > Klim) { + LocFrazilHTend(ICell, K) = 0.0_Real; + LocFrazilTTend(ICell, K) = 0.0_Real; + LocFrazilSTend(ICell, K) = 0.0_Real; + continue; + } + const Real SAIn = SA(ICell, K); const Real CTIn = CT(ICell, K); const Real PIn = P(ICell, K); @@ -184,11 +300,16 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, } // temporary log -- TBRemoved - LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " - "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", - ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), - LocAccEIce(ICell)); + if (ICell == 0) { + LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " + "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", + ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell)); + LOG_INFO(" HTend= {} " + "TTend= {} STend= {}", + HTend, TTend, STend); + } LocFrazilHTend(ICell, K) = HTend; // not scaled by dt LocFrazilTTend(ICell, K) = TTend; @@ -201,6 +322,10 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; }); + + if (conservationCheck) { + checkColumnConservation(); + } } } // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index 87ae523aebca..d10428baeaa2 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -216,9 +216,12 @@ class Frazil { Array1DReal AccMLiq; Array1DReal AccELiq; Array1DReal AccMSalt; + bool Enabled = true; void computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const Array2DReal &P, const Array2DReal &H); + bool conservationCheck = false; + Real depthLimit = -1.0_Real; private: static Frazil *DefaultFrazil; @@ -242,6 +245,8 @@ class Frazil { const HorzMesh *MeshPtr; const VertCoord *VCoordPtr; + + void checkColumnConservation() const; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index 1ae7061bd79f..edd131dfce5c 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -10,6 +10,7 @@ #include "Eos.h" #include "Field.h" #include "Forcing.h" +#include "Frazil.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" @@ -48,6 +49,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time AuxiliaryState::clear(); Forcing::clear(); OceanState::clear(); + Frazil::clear(); VertAdv::clear(); VertCoord::clear(); Dimension::clear(); @@ -61,7 +63,6 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time // destroy singletons that use raw new/delete Eos::destroyInstance(); VertMix::destroyInstance(); - return RetVal; } // end ocnFinalize diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index ba0df5dd692a..2dd13399a464 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -263,6 +263,10 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); + Err += TendConfig.get("FrazilTendencyEnable", this->FrazilTerm.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: FrazilTendencyEnable not found in TendConfig"); + Err += TendConfig.get("SurfaceTracerRestoringEnable", this->SurfaceTracerRestoring.Enabled); CHECK_ERROR_ABORT( @@ -392,7 +396,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies SfcStressForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), - CustomThicknessTend(InCustomThicknessTend), + FrazilTerm(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { // Tendency arrays @@ -823,6 +827,15 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:surfaceTracerRestoring", 2); } + if (FrazilTerm.Enabled) { + Pacer::start("Tend:frazil", 2); + const auto &PressureMid = VCoord->PressureMid; + Array2DReal PseudoThickness = State->getPseudoThickness(ThickTimeLevel); + FrazilTerm(PseudoThicknessTend, TracerTend, TracerArray, PressureMid, + PseudoThickness); + Pacer::stop("Tend:frazil", 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 5d26ca041458..866f1794b982 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -72,6 +72,7 @@ class Tendencies { TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; SurfaceTracerRestoringOnCell SurfaceTracerRestoring; + FrazilOnCell FrazilTerm; std::string Name; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 4d7ce499941f..c482002cb491 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -11,6 +11,7 @@ #include "TendencyTerms.h" #include "AuxiliaryState.h" #include "DataTypes.h" +#include "Error.h" #include "HorzMesh.h" #include "HorzOperators.h" #include "OceanState.h" @@ -115,6 +116,73 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( const HorzMesh *Mesh) {} +FrazilOnCell::FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) + : NCellsAll(Mesh->NCellsAll), TempTracerIndex(-1), SaltTracerIndex(-1), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) { + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + OMEGA_REQUIRE(TempTracerIndex >= 0, + "FrazilOnCell: Temperature tracer index is undefined"); + OMEGA_REQUIRE(SaltTracerIndex >= 0, + "FrazilOnCell: Salinity tracer index is undefined"); +} + +void FrazilOnCell::operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const { + auto *DefaultFrazil = Frazil::getDefault(); + if (!Enabled || !DefaultFrazil || !DefaultFrazil->Enabled) { + return; + } + + deepCopy(DefaultFrazil->FrazilTTend, 0.0_Real); + deepCopy(DefaultFrazil->FrazilSTend, 0.0_Real); + deepCopy(DefaultFrazil->FrazilHTend, 0.0_Real); + deepCopy(DefaultFrazil->AccMIce, 0.0_Real); + deepCopy(DefaultFrazil->AccEIce, 0.0_Real); + deepCopy(DefaultFrazil->AccMLiq, 0.0_Real); + deepCopy(DefaultFrazil->AccELiq, 0.0_Real); + deepCopy(DefaultFrazil->AccMSalt, 0.0_Real); + + const auto ConservTemp = + Kokkos::subview(TracerArray, TempTracerIndex, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, SaltTracerIndex, Kokkos::ALL, Kokkos::ALL); + + DefaultFrazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid, + PseudoThickness); + + const auto FrazilHTend = DefaultFrazil->FrazilHTend; + const auto FrazilTTend = DefaultFrazil->FrazilTTend; + const auto FrazilSTend = DefaultFrazil->FrazilSTend; + const I4 TempIndex = TempTracerIndex; + const I4 SaltIndex = SaltTracerIndex; + + OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, TracerTend); + OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); + OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); + OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); + OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); + + parallelForOuter( + {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTend(ICell, K) += LocFrazilHTend(ICell, K); + LocTracerTend(TempIndex, ICell, K) += LocFrazilTTend(ICell, K); + LocTracerTend(SaltIndex, ICell, K) += LocFrazilSTend(ICell, K); + }); + }); +} + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const VertCoord *VCoord = this->VerticalCoord; diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 389d6bb1aafb..0d46101ec398 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -11,6 +11,7 @@ //===----------------------------------------------------------------------===// #include "AuxiliaryState.h" +#include "Frazil.h" #include "GlobalConstants.h" #include "HorzMesh.h" #include "MachEnv.h" @@ -599,5 +600,26 @@ class SurfaceTracerRestoringOnCell { } }; +/// Frazil tendency hook-up term +class FrazilOnCell { + public: + bool Enabled = false; + + FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + void operator()(const Array2DReal &PseudoThicknessTend, + const Array3DReal &TracerTend, + const Array3DReal &TracerArray, + const Array2DReal &PressureMid, + const Array2DReal &PseudoThickness) const; + + private: + I4 NCellsAll; + I4 TempTracerIndex; + I4 SaltTracerIndex; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + } // namespace OMEGA #endif diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp index 98cdab844d1b..9620c85c3130 100644 --- a/components/omega/test/ocn/FrazilTest.cpp +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -157,36 +157,45 @@ void testFrazilFormationWarm() { AccELiq, AccEIce, HTend, TTend, STend); if (!isApprox(AccMIce, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero AccMIce, got {}", AccMIce); + ABORT_ERROR("FrazilFormationTest warm: expected zero AccMIce, got {}", + AccMIce); } if (!isApprox(AccMSalt, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero AccMSalt, got {}", AccMSalt); + ABORT_ERROR("FrazilFormationTest warm: expected zero AccMSalt, got {}", + AccMSalt); } if (!isApprox(AccMLiq, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero AccMLiq, got {}", AccMLiq); + ABORT_ERROR("FrazilFormationTest warm: expected zero AccMLiq, got {}", + AccMLiq); } if (!isApprox(AccELiq, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero AccELiq, got {}", AccELiq); + ABORT_ERROR("FrazilFormationTest warm: expected zero AccELiq, got {}", + AccELiq); } if (!isApprox(AccEIce, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero AccEIce, got {}", AccEIce); + ABORT_ERROR("FrazilFormationTest warm: expected zero AccEIce, got {}", + AccEIce); } if (!isApprox(HTend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero HTend, got {}", HTend); + ABORT_ERROR("FrazilFormationTest warm: expected zero HTend, got {}", + HTend); } if (!isApprox(TTend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero TTend, got {}", TTend); + ABORT_ERROR("FrazilFormationTest warm: expected zero TTend, got {}", + TTend); } if (!isApprox(STend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTest warm: expected zero STend, got {}", STend); + ABORT_ERROR("FrazilFormationTest warm: expected zero STend, got {}", + STend); } - LOG_INFO("FrazilTestWarm: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " - "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", - AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); + LOG_INFO( + "FrazilFormationTestWarm: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); } void testComputeFrazilColumn() { @@ -198,12 +207,13 @@ void testComputeFrazilColumn() { ABORT_ERROR("FrazilTestColumn: default frazil object is null"); } - const Real RTol = 1e-10_Real; - const Real SACold = 35.0_Real; - const Real PRef = 100.0_Real; - const Real HRef = 10.0_Real; - const Real CTCold = -2.0_Real; - const Real CTWarm = 10.0_Real; + const Real RTol = 1e-10_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); @@ -235,19 +245,32 @@ void testComputeFrazilColumn() { ICell); } - const I4 KBottom0 = KMax; - const I4 KBottom1 = KMax - 1; - const I4 KWarm = KMax - 2; - const I4 KTopCold = KMax - 3; - - auto CTH = createHostMirrorCopy(CT); - CTH(ICell, KBottom0) = CTCold; - CTH(ICell, KBottom1) = CTCold; - CTH(ICell, KWarm) = CTWarm; - CTH(ICell, KTopCold) = CTCold; + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - .5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; deepCopy(CT, CTH); + const bool SavedConservationCheck = TestFrazil->conservationCheck; + TestFrazil->conservationCheck = true; TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); @@ -288,11 +311,94 @@ void testComputeFrazilColumn() { ICell); } +void testComputeFrazilDepthLimit() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *TestFrazil = Frazil::getDefault(); + + if (!TestFrazil) { + ABORT_ERROR("FrazilTestColumn: default frazil object is null"); + } + + const Real RTol = 1e-10_Real; + const Real SACold = 35.0_Real; + const Real PRef = 100.0_Real; + const Real HRef = 10.0_Real; + const Real CTCold = -2.0_Real; + const Real CTWarm = 0.0_Real; + const Real CTWarm2 = -1.9_Real; + + Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers); + Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers); + Array2DReal P("P", Mesh->NCellsSize, NVertLayers); + Array2DReal H("H", Mesh->NCellsSize, NVertLayers); + + deepCopy(SA, SACold); + deepCopy(CT, CTWarm); + deepCopy(P, PRef); + deepCopy(H, HRef); + + deepCopy(TestFrazil->AccMIce, 0.0_Real); + deepCopy(TestFrazil->AccMLiq, 0.0_Real); + deepCopy(TestFrazil->AccMSalt, 0.0_Real); + deepCopy(TestFrazil->AccELiq, 0.0_Real); + deepCopy(TestFrazil->AccEIce, 0.0_Real); + deepCopy(TestFrazil->FrazilHTend, 0.0_Real); + deepCopy(TestFrazil->FrazilTTend, 0.0_Real); + deepCopy(TestFrazil->FrazilSTend, 0.0_Real); + + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + + const I4 ICell = 0; + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + if ((KMax - KMin + 1) < 4) { + ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers", + ICell); + } + + const I4 KBottom0 = KMax; + const I4 KBottom1 = KMax - 1; + const I4 KWarm = KMax - 2; + const I4 KTopCold = KMax - 3; + const I4 KCold2 = KMin + 3; + const I4 KCold3 = KMin + 2; + const I4 KWarm2 = KMin + 1; + const I4 KTopCold2 = KMin; + + auto CTH = createHostMirrorCopy(CT); + CTH(ICell, KBottom0) = CTCold; + CTH(ICell, KBottom1) = CTCold; + CTH(ICell, KWarm) = CTWarm; + CTH(ICell, KTopCold) = CTCold; + CTH(ICell, KCold2 + 2) = CTCold; + CTH(ICell, KCold2 + 1) = CTCold - .5_Real; + CTH(ICell, KCold2) = CTCold; + CTH(ICell, KCold3) = CTCold; + CTH(ICell, KWarm2) = CTWarm2; + CTH(ICell, KTopCold2) = CTCold; + deepCopy(CT, CTH); + + const bool SavedConservationCheck = TestFrazil->conservationCheck; + const bool SavedDepthLimit = TestFrazil->depthLimit; + TestFrazil->conservationCheck = true; + TestFrazil->depthLimit = 500.0_Real; + TestFrazil->computeFrazil(CT, SA, P, H); + TestFrazil->conservationCheck = SavedConservationCheck; + TestFrazil->depthLimit = SavedDepthLimit; + + auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); + auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); + auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); +} + void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { initFrazilTest(MeshFile); testFrazilFormationCold(); testFrazilFormationWarm(); testComputeFrazilColumn(); + testComputeFrazilDepthLimit(); finalizeFrazilTest(); } From c8811ae9e7c0baeb3859bd6333fae63c7679bd54 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 30 Jun 2026 11:18:30 -0700 Subject: [PATCH 06/11] update to Fraziltest incl. depthlimit --- components/omega/test/ocn/FrazilTest.cpp | 132 +++++++++++++++++------ 1 file changed, 101 insertions(+), 31 deletions(-) diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp index 9620c85c3130..ddc1b26e6b9a 100644 --- a/components/omega/test/ocn/FrazilTest.cpp +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -21,6 +21,7 @@ #include "OceanTestCommon.h" #include "OmegaKokkos.h" #include "Pacer.h" +#include "TimeMgr.h" #include "VertCoord.h" #include "mpi.h" @@ -39,11 +40,18 @@ void initFrazilTest(const std::string &mesh) { Config("Omega"); Config::readAll("omega.yml"); + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + IO::init(DefComm); - IOStream::init(); Decomp::init(mesh); + Field::init(ModelClock); + IOStream::init(ModelClock); Halo::init(); - HorzMesh::init(); + HorzMesh::init(ModelClock); VertCoord::init(false); Frazil::init(); } @@ -56,9 +64,15 @@ void finalizeFrazilTest() { Decomp::clear(); Field::clear(); Dimension::clear(); + IOStream::finalize(); MachEnv::removeAll(); } +// this test only excercises the frazil formation functor (no melt) +// in a cold case: the frazil terms should be +// - strictly positive for ice, liquid, and salt mass +//- strictly negative for ice and liquid energy +// - positive for T tendency and negative for S, H tendencies void testFrazilFormationCold() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -89,44 +103,54 @@ void testFrazilFormationCold() { AccELiq, AccEIce, HTend, TTend, STend); if (AccMIce <= 0.0_Real) { - ABORT_ERROR("FrazilTestCold: accumulated ice mass is non-positive: {}", - AccMIce); + ABORT_ERROR( + "FrazilFormationTestCold: accumulated ice mass is non-positive: {}", + AccMIce); } if (AccMLiq <= 0.0_Real) { - ABORT_ERROR("FrazilTestCold: accumulated liquid mass is non-positive: {}", + ABORT_ERROR("FrazilFormationTestCold: accumulated liquid mass is " + "non-positive: {}", AccMLiq); } if (AccMSalt <= 0.0_Real) { - ABORT_ERROR("FrazilTestCold: accumulated salt mass is non-positive: {}", - AccMSalt); + ABORT_ERROR( + "FrazilFormationTestCold: accumulated salt mass is non-positive: {}", + AccMSalt); } - if (isApprox(AccELiq, 0.0_Real, RTol)) { - ABORT_ERROR( - "FrazilTestCold: accumulated liquid energy is effectively zero: {}", - AccELiq); + if (AccELiq >= 0.0_Real) { + ABORT_ERROR("FrazilFormationTestCold: accumulated liquid energy is " + "positive (exp. negative): {}", + AccELiq); } - if (isApprox(AccEIce, 0.0_Real, RTol)) { - ABORT_ERROR( - "FrazilTestCold: accumulated ice energy is effectively zero: {}", - AccEIce); + if (AccEIce >= 0.0_Real) { + ABORT_ERROR("FrazilFormationTestCold: accumulated ice energy is positive " + "(exp. negative): {}", + AccEIce); } - if (isApprox(HTend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTestCold: HTend is effectively zero: {}", HTend); + if (HTend >= 0.0_Real) { + ABORT_ERROR( + "FrazilFormationTestCold: HTend is positive (exp. negative): {}", + HTend); } - - if (isApprox(TTend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTestCold: TTend is zero: {}", TTend); + if (TTend <= 0.0_Real) { + ABORT_ERROR( + "FrazilFormationTestCold: TTend is negative (exp. positive): {}", + TTend); } - - if (isApprox(STend, 0.0_Real, RTol)) { - ABORT_ERROR("FrazilTestCold: STend is effectively zero: {}", STend); + if (STend >= 0.0_Real) { + ABORT_ERROR( + "FrazilFormationTestCold: STend is positive (exp. negative): {}", + STend); } - LOG_INFO("FrazilTestCold: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " - "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", - AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); + LOG_INFO( + "FrazilFormationTestCold: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, " + "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}", + AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); } +// this test only excercises the frazil formation functor (no melt) +// in a warm case: the frazil FORMATION terms should all be zero void testFrazilFormationWarm() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -198,6 +222,10 @@ void testFrazilFormationWarm() { AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend); } +// this test exercises the frazil formation and melt functors +// in a column of water with both cold and warm layers. +// It turns to frazil column conservation check. +// In dev, there is extra verbose logging in the frazil code (TBRemoved). void testComputeFrazilColumn() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -311,6 +339,9 @@ void testComputeFrazilColumn() { ICell); } +// this test exercises the frazil formation and melt functors +// with a depth limit set. Layers deeper than the depth limit +// should have zero frazil tendencies. void testComputeFrazilDepthLimit() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -320,7 +351,7 @@ void testComputeFrazilDepthLimit() { ABORT_ERROR("FrazilTestColumn: default frazil object is null"); } - const Real RTol = 1e-10_Real; + const Real RTol = 1e-12_Real; const Real SACold = 35.0_Real; const Real PRef = 100.0_Real; const Real HRef = 10.0_Real; @@ -353,8 +384,8 @@ void testComputeFrazilDepthLimit() { const I4 ICell = 0; const I4 KMin = MinLayerCellH(ICell); const I4 KMax = MaxLayerCellH(ICell); - if ((KMax - KMin + 1) < 4) { - ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers", + if ((KMax - KMin + 1) < 10) { + ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 10 active layers", ICell); } @@ -382,8 +413,20 @@ void testComputeFrazilDepthLimit() { const bool SavedConservationCheck = TestFrazil->conservationCheck; const bool SavedDepthLimit = TestFrazil->depthLimit; - TestFrazil->conservationCheck = true; - TestFrazil->depthLimit = 500.0_Real; + const Real TestDepthLimit = 35.0_Real; // this needs to be positive + // if TestDepthLimit is negative, test will fail: + // - the code assume depthlimit < 0 mean no limit (i.e. full depth frazil) + // - the test below will exclude all layers and fail because Tend !=0. + + // Populate GeomZMid explicitly for the test column. + auto GeomZMidH = createHostMirrorCopy(VCoord->GeomZMid); + for (I4 K = KMin; K <= KMax; ++K) { + GeomZMidH(ICell, K) = -10.0_Real * (K - KMin + 1); + } + deepCopy(VCoord->GeomZMid, GeomZMidH); + + TestFrazil->conservationCheck = true; + TestFrazil->depthLimit = TestDepthLimit; TestFrazil->computeFrazil(CT, SA, P, H); TestFrazil->conservationCheck = SavedConservationCheck; TestFrazil->depthLimit = SavedDepthLimit; @@ -391,6 +434,33 @@ void testComputeFrazilDepthLimit() { auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend); auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend); auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend); + + bool FoundExcludedLayer = false; + for (I4 K = KMin; K <= KMax; ++K) { + const Real Depth = GeomZMidH(ICell, K); + const Real AbsDepth = Depth < 0.0_Real ? -Depth : Depth; + + if (AbsDepth > TestDepthLimit) { + FoundExcludedLayer = true; + if (!isApprox(HTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(TTendH(ICell, K), 0.0_Real, RTol) || + !isApprox(STendH(ICell, K), 0.0_Real, RTol)) { + ABORT_ERROR("FrazilDepthLimitTest: excluded layer K={} has " + "non-zero tendencies (HTend={}, TTend={}, STend={})", + K, HTendH(ICell, K), TTendH(ICell, K), + STendH(ICell, K)); + } + } + } + if (!FoundExcludedLayer) { + ABORT_ERROR("FrazilDepthLimitTest: no layers were excluded for ICell={} " + "with depthLimit={}", + ICell, TestDepthLimit); + } + + LOG_INFO("FrazilDepthLimitTest: depthLimit={} exclusion check passed for " + "ICell={}", + TestDepthLimit, ICell); } void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") { From 259fbc229c9d78b9cd805ce9840ce06f9800e4b1 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 30 Jun 2026 14:47:13 -0700 Subject: [PATCH 07/11] clean up of redundant flag. Frazil is based on FrazilTendencyEnable --- components/omega/configs/Default.yml | 1 - components/omega/src/ocn/Frazil.cpp | 42 ++++++++++++++-------- components/omega/src/ocn/Frazil.h | 1 - components/omega/src/ocn/TendencyTerms.cpp | 30 ++++++++-------- 4 files changed, 43 insertions(+), 31 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index bb49ff8c6bbb..86b032a5f77c 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -78,7 +78,6 @@ Omega: DRhoDS: 0.8 RhoT0S0: 1000.0 Frazil: - Enable: true FrazilType: teos MassLimit: 0.1 Phi: 0.75 diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index d2db57c8031f..26d0e0607363 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -42,6 +42,26 @@ void Frazil::init() { } if (!DefaultFrazil) { + Error Err; + bool FrazilTendencyEnable = false; + Config *OmegaConfig = Config::getOmegaConfig(); + Config TendConfig("Tendencies"); + + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, + "Frazil::init: Tendencies group not found in Config"); + + Err += TendConfig.get("FrazilTendencyEnable", FrazilTendencyEnable); + CHECK_ERROR_ABORT( + Err, "Frazil::init: FrazilTendencyEnable not found in Tendencies"); + + if (!FrazilTendencyEnable) { + LOG_INFO("Frazil::init: Frazil tendency disabled; skipping default " + "frazil object creation"); + LOG_INFO("All frazil is off - frazil parameters will be ignored"); + return; + } + DefaultFrazil = create("Default"); } } @@ -101,9 +121,11 @@ Frazil *Frazil::create(const std::string &Name) { if ((FrazilTypeStr == "Basic") or (FrazilTypeStr == "basic") or (FrazilTypeStr == "BasicFrazil")) { NewFrazil->frazilChoice = FrazilType::BasicFrazil; + ABORT_ERROR("Frazil::create: BasicFrazil not supported yet"); } else if ((FrazilTypeStr == "Simple") or (FrazilTypeStr == "simple") or (FrazilTypeStr == "SimpleFrazil")) { NewFrazil->frazilChoice = FrazilType::SimpleFrazil; + ABORT_ERROR("Frazil::create: SimpleFrazil not supported yet"); } else if ((FrazilTypeStr == "teos") or (FrazilTypeStr == "Teos") or (FrazilTypeStr == "TEOS") or (FrazilTypeStr == "Teos10") or (FrazilTypeStr == "teos10") or (FrazilTypeStr == "TEOS10")) { @@ -120,21 +142,13 @@ Frazil *Frazil::create(const std::string &Name) { Err += FrazilConfig.get("Phi", NewFrazil->computeFrazilFormation.Phi); CHECK_ERROR_ABORT(Err, "Frazil::create: Phi not found in Frazil config"); - Error CheckColumnErr = - FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); - if (!CheckColumnErr.isSuccess()) { - NewFrazil->conservationCheck = false; - } + Err += FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); + CHECK_ERROR_ABORT( + Err, "Frazil::create: ConservationCheck not found in Frazil config"); - Error EnabledErr = FrazilConfig.get("Enable", NewFrazil->Enabled); - if (!EnabledErr.isSuccess()) { - NewFrazil->Enabled = true; - } - - Error DepthLimitErr = FrazilConfig.get("DepthLimit", NewFrazil->depthLimit); - if (!DepthLimitErr.isSuccess()) { - NewFrazil->depthLimit = -1.0_Real; - } + Err += FrazilConfig.get("DepthLimit", NewFrazil->depthLimit); + CHECK_ERROR_ABORT(Err, + "Frazil::create: DepthLimit not found in Frazil config"); if (Name == "Default") { DefaultFrazil = NewFrazil; diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index d10428baeaa2..ddbdd904fd3c 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -216,7 +216,6 @@ class Frazil { Array1DReal AccMLiq; Array1DReal AccELiq; Array1DReal AccMSalt; - bool Enabled = true; void computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const Array2DReal &P, const Array2DReal &H); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index c482002cb491..f62f14430dc3 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -133,31 +133,31 @@ void FrazilOnCell::operator()(const Array2DReal &PseudoThicknessTend, const Array3DReal &TracerArray, const Array2DReal &PressureMid, const Array2DReal &PseudoThickness) const { - auto *DefaultFrazil = Frazil::getDefault(); - if (!Enabled || !DefaultFrazil || !DefaultFrazil->Enabled) { + auto *Frazil = Frazil::getDefault(); + if (!Enabled || !Frazil) { return; } - deepCopy(DefaultFrazil->FrazilTTend, 0.0_Real); - deepCopy(DefaultFrazil->FrazilSTend, 0.0_Real); - deepCopy(DefaultFrazil->FrazilHTend, 0.0_Real); - deepCopy(DefaultFrazil->AccMIce, 0.0_Real); - deepCopy(DefaultFrazil->AccEIce, 0.0_Real); - deepCopy(DefaultFrazil->AccMLiq, 0.0_Real); - deepCopy(DefaultFrazil->AccELiq, 0.0_Real); - deepCopy(DefaultFrazil->AccMSalt, 0.0_Real); + deepCopy(Frazil->FrazilTTend, 0.0_Real); + deepCopy(Frazil->FrazilSTend, 0.0_Real); + deepCopy(Frazil->FrazilHTend, 0.0_Real); + deepCopy(Frazil->AccMIce, 0.0_Real); + deepCopy(Frazil->AccEIce, 0.0_Real); + deepCopy(Frazil->AccMLiq, 0.0_Real); + deepCopy(Frazil->AccELiq, 0.0_Real); + deepCopy(Frazil->AccMSalt, 0.0_Real); const auto ConservTemp = Kokkos::subview(TracerArray, TempTracerIndex, Kokkos::ALL, Kokkos::ALL); const auto AbsSalinity = Kokkos::subview(TracerArray, SaltTracerIndex, Kokkos::ALL, Kokkos::ALL); - DefaultFrazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid, - PseudoThickness); + Frazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid, + PseudoThickness); - const auto FrazilHTend = DefaultFrazil->FrazilHTend; - const auto FrazilTTend = DefaultFrazil->FrazilTTend; - const auto FrazilSTend = DefaultFrazil->FrazilSTend; + const auto FrazilHTend = Frazil->FrazilHTend; + const auto FrazilTTend = Frazil->FrazilTTend; + const auto FrazilSTend = Frazil->FrazilSTend; const I4 TempIndex = TempTracerIndex; const I4 SaltIndex = SaltTracerIndex; From 51d02e5b5c6a2d803b3d0cdf347087cf790e8bc8 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 30 Jun 2026 17:55:35 -0700 Subject: [PATCH 08/11] fix to P units + tendencytest --- components/omega/src/ocn/Frazil.cpp | 10 +- components/omega/test/infra/IOStreamTest.cpp | 7 +- components/omega/test/ocn/TendenciesTest.cpp | 148 ++++++++++++++++++- 3 files changed, 156 insertions(+), 9 deletions(-) diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 26d0e0607363..479335a046b2 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -293,21 +293,22 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const Real SAIn = SA(ICell, K); const Real CTIn = CT(ICell, K); const Real PIn = P(ICell, K); + const Real PDb = PIn * Pa2Db; const Real H = LayerH(ICell, K); - const Real Tfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real); + const Real Tfrz = gsw_ct_freezing_poly(SAIn, PDb, 0.0_Real); Real HTend = 0.0_Real; Real TTend = 0.0_Real; Real STend = 0.0_Real; if (CTIn < Tfrz) { - LocComputeFrazilFormation(SAIn, CTIn, PIn, H, LocAccMIce(ICell), + LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell), LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), LocAccEIce(ICell), HTend, TTend, STend); } else { - LocComputeFrazilMelt(SAIn, CTIn, PIn, H, LocAccMIce(ICell), + LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), LocAccMLiq(ICell), LocAccMSalt(ICell), LocAccELiq(ICell), LocAccEIce(ICell), HTend, TTend, STend); @@ -315,6 +316,9 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, // temporary log -- TBRemoved if (ICell == 0) { + LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, PIn = " + "{}, H = {}, Tfrz = {}", + ICell, SAIn, CTIn, PIn, H, Tfrz); LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 1678a08470f6..51620505568c 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -108,12 +108,13 @@ void initIOStreamTest(Clock *&ModelClock // Model clock PressureGrad::init(); + // Initialize Tracers before Tendencies so tracer indices are available + // during FrazilOnCell construction inside Tendencies::init(). + Tracers::init(); + // Intialize Tendencies Tendencies::init(); - // Initialize Tracers - Tracers::init(); - // Initialize Aux State AuxiliaryState::init(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5c8bb8d3d3f0..5df40c343a2c 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -138,6 +138,7 @@ int initTendenciesTest(const std::string &mesh) { VertCoord::init(); Tracers::init(); + Frazil::init(); VertAdv::init(); PressureGrad::init(); Eos::init(); @@ -241,6 +242,7 @@ int testTendencies() { "BaselineNormalVelocityTend", Mesh->NEdgesSize, VCoord->NVertLayers); DefTendencies->SfcStressForcing.Enabled = false; + DefTendencies->FrazilTerm.Enabled = false; // frazil needs to be off for now DefTendencies->computeAllTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, TracerTimeLevel, Time, Interval); @@ -302,6 +304,145 @@ int testTendencies() { DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled; + // Verify frazil tendencies are plumbed into pseudo-thickness and tracer + // tendencies by comparing runs with FrazilTerm disabled/enabled. + I4 TempTracerIndex = -1; + I4 SaltTracerIndex = -1; + Tracers::getIndex(TempTracerIndex, "Temperature"); + Tracers::getIndex(SaltTracerIndex, "Salinity"); + + if (TempTracerIndex < 0 || SaltTracerIndex < 0) { + Err++; + LOG_ERROR("TendenciesTest: missing Temperature/Salinity tracer indices"); + } else { + auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell); + auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell); + auto TracerArrayH = createHostMirrorCopy(TracerArray); + + // Set Salinity in active layers to [31, 33], top-to-bottom ramp. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + + for (I4 K = KMin; K <= KMax; ++K) { + Real frac = 0.0_Real; + if (KMax > KMin) { + frac = + static_cast(K - KMin) / static_cast(KMax - KMin); + } + TracerArrayH(SaltTracerIndex, ICell, K) = + 31.0_Real + 2.0_Real * frac; + } + } + + // Optional: keep Temperature warm everywhere first, then add one cold + // point. + for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) { + const I4 KMin = MinLayerCellH(ICell); + const I4 KMax = MaxLayerCellH(ICell); + for (I4 K = KMin; K <= KMax; ++K) { + TracerArrayH(TempTracerIndex, ICell, K) = 8.0_Real; + } + } + const I4 FrazilTestCell = 0; + const I4 FrazilTestLayer = MinLayerCellH(FrazilTestCell); + TracerArrayH(TempTracerIndex, FrazilTestCell, FrazilTestLayer) = + -2.0_Real; + deepCopy(TracerArray, TracerArrayH); + + LOG_INFO("TendenciesTest: minLayerCell = {}", FrazilTestLayer); + + Array2DReal BaselinePseudoThicknessTend( + "BaselinePseudoThicknessTend", Mesh->NCellsSize, VCoord->NVertLayers); + Array3DReal BaselineTracerTend("BaselineTracerTend", + Tracers::getNumTracers(), Mesh->NCellsSize, + VCoord->NVertLayers); + + const bool OrigFrazilEnabled = DefTendencies->FrazilTerm.Enabled; + + DefTendencies->FrazilTerm.Enabled = false; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + deepCopy(BaselinePseudoThicknessTend, DefTendencies->PseudoThicknessTend); + deepCopy(BaselineTracerTend, DefTendencies->TracerTend); + + deepCopy(DefTendencies->PseudoThicknessTend, 0.0_Real); + deepCopy(DefTendencies->NormalVelocityTend, 0.0_Real); + deepCopy(DefTendencies->TracerTend, 0.0_Real); + + DefTendencies->FrazilTerm.Enabled = true; + DefTendencies->computeAllTendencies(State, AuxState, TracerArray, + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time, Interval); + + Array2DReal PseudoThicknessTendDiff( + "PseudoThicknessTendDiff", Mesh->NCellsSize, VCoord->NVertLayers); + Array2DReal TempTracerTendDiff("TempTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + Array2DReal SaltTracerTendDiff("SaltTracerTendDiff", Mesh->NCellsSize, + VCoord->NVertLayers); + + OMEGA_SCOPE(LocPseudoThicknessTendDiff, PseudoThicknessTendDiff); + OMEGA_SCOPE(LocTempTracerTendDiff, TempTracerTendDiff); + OMEGA_SCOPE(LocSaltTracerTendDiff, SaltTracerTendDiff); + OMEGA_SCOPE(LocPseudoThicknessTend, DefTendencies->PseudoThicknessTend); + OMEGA_SCOPE(LocBaselinePseudoThicknessTend, BaselinePseudoThicknessTend); + OMEGA_SCOPE(LocTracerTend, DefTendencies->TracerTend); + OMEGA_SCOPE(LocBaselineTracerTend, BaselineTracerTend); + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + + parallelForOuter( + "TendenciesTest:FrazilTendDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + + parallelForInner( + Team, Range{KMin, KMax}, INNER_LAMBDA(int K) { + LocPseudoThicknessTendDiff(ICell, K) = + Kokkos::abs(LocPseudoThicknessTend(ICell, K) - + LocBaselinePseudoThicknessTend(ICell, K)); + LocTempTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(TempTracerIndex, ICell, K) - + LocBaselineTracerTend(TempTracerIndex, ICell, K)); + LocSaltTracerTendDiff(ICell, K) = Kokkos::abs( + LocTracerTend(SaltTracerIndex, ICell, K) - + LocBaselineTracerTend(SaltTracerIndex, ICell, K)); + }); + }); + + const Real PseudoDelta = sum(PseudoThicknessTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real TempDelta = sum(TempTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + const Real SaltDelta = sum(SaltTracerTendDiff, Mesh->NCellsOwned, + VCoord->MinLayerCell, VCoord->MaxLayerCell); + + constexpr Real FrazilDeltaATol = 1e-12_Real; + if (!Kokkos::isfinite(PseudoDelta) || + isApprox(PseudoDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change " + "PseudoThicknessTend"); + } + if (!Kokkos::isfinite(TempDelta) || + isApprox(TempDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Temperature tracer " + "tendency"); + } + if (!Kokkos::isfinite(SaltDelta) || + isApprox(SaltDelta, 0._Real, 0._Real, FrazilDeltaATol)) { + Err++; + LOG_ERROR("TendenciesTest: Frazil did not change Salinity tracer " + "tendency"); + } + + DefTendencies->FrazilTerm.Enabled = OrigFrazilEnabled; + } + // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; int NEdgesOwned = Mesh->NEdgesOwned; @@ -323,12 +464,12 @@ int testTendencies() { LOG_ERROR("TendenciesTest: NormVelTendSum FAIL"); } - const Real TraceTendSum = + const Real TracerTendSum = sum(DefTendencies->TracerTend, NTracers, NCellsOwned, VCoord->MinLayerCell, VCoord->MaxLayerCell); - if (!Kokkos::isfinite(TraceTendSum) || TraceTendSum == 0) { + if (!Kokkos::isfinite(TracerTendSum) || TracerTendSum == 0) { Err++; - LOG_ERROR("TendenciesTest: TraceTendSum FAIL"); + LOG_ERROR("TendenciesTest: TracerTendSum FAIL"); } Tendencies::clear(); @@ -339,6 +480,7 @@ int testTendencies() { void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); + Frazil::clear(); PressureGrad::clear(); Eos::destroyInstance(); AuxiliaryState::clear(); From d85d5d7599b228f3f26fd3f4da71bc31db3ba9d1 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 9 Jul 2026 13:42:34 -0700 Subject: [PATCH 09/11] using the Eos for CtFrz calculation --- components/omega/src/ocn/Eos.h | 18 ++- components/omega/src/ocn/Frazil.cpp | 180 ++++++++++++++++------------ components/omega/src/ocn/Frazil.h | 10 +- 3 files changed, 128 insertions(+), 80 deletions(-) diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index a530549d493d..effaefec3ada 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -354,8 +354,8 @@ class Teos10Eos { /// Calculates freezing Conservative Temperature using TEOS-10 polynomial /// (polynomial error in [-5e-4, 6e-4] K, from GSW package) - KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, - const Real SaturationFract) const { + KOKKOS_FUNCTION static Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract) { constexpr Real Sso = 35.16504; constexpr Real C0 = 0.017947064327968736; constexpr Real C1 = -6.076099099929818; @@ -756,6 +756,20 @@ class Eos { /// Convert potential temperature to Conservative Temperature Real calcCtFromPt(const Real &Sa, const Real &Pt) const; + /// Calculate Conservative Temperature at freezing point. + /// P is expected in dbar to match TEOS polynomial convention. + KOKKOS_FUNCTION static Real calcCtFreezing(const Real Sa, const Real P, + const Real SaturationFract, + const EosType Choice) { + if (Choice == EosType::Teos10Eos) { + return Teos10Eos::calcCtFreezing(Sa, P, SaturationFract); + } + + Kokkos::abort("Eos::calcCtFreezing: CtFreezing not implemented for " + "non-TEOS-10 EOS"); + return 0.0_Real; + } + /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 479335a046b2..1d05d2193b53 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -6,6 +6,7 @@ //===----------------------------------------------------------------------===// #include "Frazil.h" +#include "Eos.h" #include "Error.h" #include "Logging.h" @@ -41,6 +42,9 @@ void Frazil::init() { ABORT_ERROR("Frazil::init: HorzMesh and VertCoord must be initialized"); } + // Frazil freezing-temperature calculations depend on EOS configuration. + Eos::init(); + if (!DefaultFrazil) { Error Err; bool FrazilTendencyEnable = false; @@ -228,17 +232,17 @@ void Frazil::checkColumnConservation() const { } if (!isApprox(-MassTend * RhoSw, MassTotal, RTol)) { - ABORT_ERROR( + LOG_INFO( "Frazil column mass check failed: cell {} tendency={} total={}", ICell, -MassTend * RhoSw, MassTotal); } if (!isApprox(-EnergyTend * Cp0Sw * RhoSw, EnergyTotal, RTol)) { - ABORT_ERROR( + LOG_INFO( "Frazil column energy check failed: cell {} tendency={} total={}", ICell, -EnergyTend * Cp0Sw * RhoSw, EnergyTotal); } if (!isApprox(-SaltTend * RhoSw * PPt2Salt, SaltTotal, RTol)) { - ABORT_ERROR( + LOG_INFO( "Frazil column salt check failed: cell {} tendency={} total={}", ICell, -SaltTend * RhoSw * PPt2Salt, SaltTotal); } @@ -247,12 +251,26 @@ void Frazil::checkColumnConservation() const { void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, const Array2DReal &P, const Array2DReal &LayerH) { + Eos *DefEos = Eos::getInstance(); + if (!DefEos) { + ABORT_ERROR("Frazil::computeFrazil: Eos must be initialized before " + "computeFrazil"); + } + + const EosType LocEosChoice = DefEos->EosChoice; + if (LocEosChoice != EosType::Teos10Eos) { + ABORT_ERROR("Frazil::computeFrazil: CtFreezing not implemented for " + "non-TEOS-10 EOS"); + } + OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell); OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid); OMEGA_SCOPE(LocComputeFrazilFormation, computeFrazilFormation); OMEGA_SCOPE(LocComputeFrazilMelt, computeFrazilMelt); + OMEGA_SCOPE(LocComputeBasicFrazilFormation, computeBasicFrazilFormation); + OMEGA_SCOPE(LocComputeBasicFrazilMelt, computeBasicFrazilMelt); OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); @@ -264,86 +282,98 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, parallelFor( {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { - const I4 KMin = MinLayerCell(ICell); - const I4 KMax = MaxLayerCell(ICell); - - I4 Klim = KMax; - bool HasKlim = true; - const bool Limit = (depthLimit >= 0.0_Real); - if (Limit) { - HasKlim = false; - for (I4 K = KMax; K >= KMin; --K) { - if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) { - Klim = K; - HasKlim = true; - break; - } - } - } - - // Explicit accumulation order: bottom layer to top layer. - for (I4 K = KMax; K >= KMin; --K) { - if (!HasKlim || K > Klim) { - LocFrazilHTend(ICell, K) = 0.0_Real; - LocFrazilTTend(ICell, K) = 0.0_Real; - LocFrazilSTend(ICell, K) = 0.0_Real; - continue; - } + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + I4 Klim = KMax; + bool HasKlim = true; + const bool Limit = (depthLimit >= 0.0_Real); + + if (Limit) { + HasKlim = false; + for (I4 K = KMax; K >= KMin; --K) { + if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) { + Klim = K; + HasKlim = true; + break; + } + } + } - const Real SAIn = SA(ICell, K); - const Real CTIn = CT(ICell, K); - const Real PIn = P(ICell, K); - const Real PDb = PIn * Pa2Db; - const Real H = LayerH(ICell, K); - - const Real Tfrz = gsw_ct_freezing_poly(SAIn, PDb, 0.0_Real); - - Real HTend = 0.0_Real; - Real TTend = 0.0_Real; - Real STend = 0.0_Real; - - if (CTIn < Tfrz) { - LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), - LocAccELiq(ICell), LocAccEIce(ICell), - HTend, TTend, STend); - } else { - LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), - LocAccELiq(ICell), LocAccEIce(ICell), - HTend, TTend, STend); + // Explicit accumulation order: bottom layer to top layer. + for (I4 K = KMax; K >= KMin; --K) { + if (!HasKlim || K > Klim) { + LocFrazilHTend(ICell, K) = 0.0_Real; + LocFrazilTTend(ICell, K) = 0.0_Real; + LocFrazilSTend(ICell, K) = 0.0_Real; + continue; + } + + const Real SAIn = SA(ICell, K); + const Real CTIn = CT(ICell, K); + const Real PIn = P(ICell, K); + const Real PDb = PIn * Pa2Db; + const Real H = LayerH(ICell, K); + + const Real Tfrz = + Eos::calcCtFreezing(SAIn, PDb, 0.0_Real, LocEosChoice); + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + if (CTIn < Tfrz) { + LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), + HTend, TTend, STend); + } + } + else { + if (LocAccMIce(ICell) > 0.0_Real) { + LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), HTend, + TTend, STend); + } + } + // else { + // temporary kernel logging - TBRemoved + // LOG_INFO( + // "warm layer but no ice to melt"); + // } } // temporary log -- TBRemoved - if (ICell == 0) { - LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, PIn = " - "{}, H = {}, Tfrz = {}", - ICell, SAIn, CTIn, PIn, H, Tfrz); - LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " - "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", - ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), - LocAccELiq(ICell), LocAccEIce(ICell)); - LOG_INFO(" HTend= {} " - "TTend= {} STend= {}", - HTend, TTend, STend); - } + // if (ICell == 0) { + // LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, PIn = " + // "{}, H = {}, Tfrz = {}", + // ICell, SAIn, CTIn, PIn, H, Tfrz); + // LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " + // "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", + // ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), + // LocAccMLiq(ICell), LocAccMSalt(ICell), + // LocAccELiq(ICell), LocAccEIce(ICell)); + // LOG_INFO(" HTend= {} " + // "TTend= {} STend= {}", + // HTend, TTend, STend); + // } LocFrazilHTend(ICell, K) = HTend; // not scaled by dt LocFrazilTTend(ICell, K) = TTend; LocFrazilSTend(ICell, K) = STend; - } - // Convert to coupler units - LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; - LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; - LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; - LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; - LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; - }); - - if (conservationCheck) { - checkColumnConservation(); - } +} +// Convert to coupler units +LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; +LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; +LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; +LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; +LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; +}); + +if (conservationCheck) { + checkColumnConservation(); +} } } // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index ddbdd904fd3c..1121338ef284 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -23,7 +23,7 @@ namespace OMEGA { enum class FrazilType { - BasicFrazil, ///< Placeholder basic frazil option + BasicFrazil, ///< MPAS-O style basic frazil option SimpleFrazil, ///< Placeholder simple frazil option TeosFrazil ///< Placeholder TEOS frazil option }; @@ -164,8 +164,10 @@ class FrazilFormation { wIh = static_cast(wIh_d); const Real OneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - Phi); - solidMass = h * Kokkos::min(wIh, OneMinusPhi * MassLimit); - liquidMass = (Phi / OneMinusPhi) * solidMass; + // anything called mass below is in pseudo-thickness units (m) and needs + // to be scaled by RhoSw for coupling + solidMass = h * Kokkos::min(wIh, OneMinusPhi * MassLimit); + liquidMass = (Phi / OneMinusPhi) * solidMass; solidEnthalpy = solidMass * gsw_pot_enthalpy_from_pt_ice_poly(PTnew_d); liquidEnthalpy = liquidMass * Cp0Sw * CTnew; @@ -235,6 +237,8 @@ class Frazil { Frazil &operator=(Frazil &&) = delete; FrazilType frazilChoice; + BasicFrazilFormation computeBasicFrazilFormation; + BasicFrazilMelt computeBasicFrazilMelt; FrazilFormation computeFrazilFormation; FrazilMelt computeFrazilMelt; Real massLimit; From a7c3be012a01cdcca9e4aa2d354f61ba15cdb8e2 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 9 Jul 2026 18:02:19 -0700 Subject: [PATCH 10/11] added limit to melt, refactored melt terms, and clean-up --- components/omega/src/ocn/Frazil.cpp | 189 ++++++++++--------- components/omega/src/ocn/Frazil.h | 162 ++++++++-------- components/omega/test/ocn/FrazilTest.cpp | 4 + components/omega/test/ocn/TendenciesTest.cpp | 2 +- 4 files changed, 183 insertions(+), 174 deletions(-) diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 1d05d2193b53..9d0f1ee0b386 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -1,7 +1,7 @@ //===-- ocn/Frazil.cpp - Frazil Ice Formation -----------------*- C++ -*-===// // // The Frazil class manages frazil tendencies and accumulators. -// This initial scaffold wires allocation and configuration only. +// This initial implementation only has a teos-10 configuration. // //===----------------------------------------------------------------------===// @@ -71,8 +71,7 @@ void Frazil::init() { } Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) - : frazilChoice(FrazilType::TeosFrazil), massLimit(0.1_Real), phi(0.75_Real), - NCellsAll(Mesh->NCellsAll), + : frazilChoice(FrazilType::TeosFrazil), NCellsAll(Mesh->NCellsAll), NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh), VCoordPtr(VCoord), computeFrazilFormation(), computeFrazilMelt() { @@ -139,11 +138,12 @@ Frazil *Frazil::create(const std::string &Name) { } Err += FrazilConfig.get("MassLimit", - NewFrazil->computeFrazilFormation.MassLimit); + NewFrazil->computeFrazilFormation.massLimit); + Err += FrazilConfig.get("MassLimit", NewFrazil->computeFrazilMelt.massLimit); CHECK_ERROR_ABORT(Err, "Frazil::create: MassLimit not found in Frazil config"); - Err += FrazilConfig.get("Phi", NewFrazil->computeFrazilFormation.Phi); + Err += FrazilConfig.get("Phi", NewFrazil->computeFrazilFormation.phi); CHECK_ERROR_ABORT(Err, "Frazil::create: Phi not found in Frazil config"); Err += FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck); @@ -269,8 +269,6 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, OMEGA_SCOPE(LocComputeFrazilFormation, computeFrazilFormation); OMEGA_SCOPE(LocComputeFrazilMelt, computeFrazilMelt); - OMEGA_SCOPE(LocComputeBasicFrazilFormation, computeBasicFrazilFormation); - OMEGA_SCOPE(LocComputeBasicFrazilMelt, computeBasicFrazilMelt); OMEGA_SCOPE(LocFrazilTTend, FrazilTTend); OMEGA_SCOPE(LocFrazilSTend, FrazilSTend); OMEGA_SCOPE(LocFrazilHTend, FrazilHTend); @@ -282,98 +280,105 @@ void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA, parallelFor( {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { - const I4 KMin = MinLayerCell(ICell); - const I4 KMax = MaxLayerCell(ICell); - - I4 Klim = KMax; - bool HasKlim = true; - const bool Limit = (depthLimit >= 0.0_Real); - - if (Limit) { - HasKlim = false; - for (I4 K = KMax; K >= KMin; --K) { - if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) { - Klim = K; - HasKlim = true; - break; - } - } - } + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + I4 Klim = KMax; + bool HasKlim = true; + const bool Limit = (depthLimit >= 0.0_Real); + + if (Limit) { + HasKlim = false; + for (I4 K = KMax; K >= KMin; --K) { + if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) { + Klim = K; + HasKlim = true; + break; + } + } + } + + // Explicit accumulation order: bottom layer to top layer. + for (I4 K = KMax; K >= KMin; --K) { + if (!HasKlim || K > Klim) { + LocFrazilHTend(ICell, K) = 0.0_Real; + LocFrazilTTend(ICell, K) = 0.0_Real; + LocFrazilSTend(ICell, K) = 0.0_Real; + continue; + } - // Explicit accumulation order: bottom layer to top layer. - for (I4 K = KMax; K >= KMin; --K) { - if (!HasKlim || K > Klim) { - LocFrazilHTend(ICell, K) = 0.0_Real; - LocFrazilTTend(ICell, K) = 0.0_Real; - LocFrazilSTend(ICell, K) = 0.0_Real; - continue; - } - - const Real SAIn = SA(ICell, K); - const Real CTIn = CT(ICell, K); - const Real PIn = P(ICell, K); - const Real PDb = PIn * Pa2Db; - const Real H = LayerH(ICell, K); - - const Real Tfrz = - Eos::calcCtFreezing(SAIn, PDb, 0.0_Real, LocEosChoice); - - Real HTend = 0.0_Real; - Real TTend = 0.0_Real; - Real STend = 0.0_Real; - - if (CTIn < Tfrz) { - LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), - LocAccELiq(ICell), LocAccEIce(ICell), - HTend, TTend, STend); - } - } - else { - if (LocAccMIce(ICell) > 0.0_Real) { - LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), - LocAccMLiq(ICell), LocAccMSalt(ICell), - LocAccELiq(ICell), LocAccEIce(ICell), HTend, - TTend, STend); - } - } - // else { - // temporary kernel logging - TBRemoved - // LOG_INFO( - // "warm layer but no ice to melt"); - // } + const Real SAIn = SA(ICell, K); + const Real CTIn = CT(ICell, K); + const Real PIn = P(ICell, K); + const Real PDb = PIn * Pa2Db; + const Real H = LayerH(ICell, K); + + const Real Tfrz = + Eos::calcCtFreezing(SAIn, PDb, 0.0_Real, LocEosChoice); + + Real HTend = 0.0_Real; + Real TTend = 0.0_Real; + Real STend = 0.0_Real; + + if (CTIn < Tfrz) { + LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), + HTend, TTend, STend); + } + + else { + if (LocAccMIce(ICell) > 0.0_Real) { + LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell), + LocAccMLiq(ICell), LocAccMSalt(ICell), + LocAccELiq(ICell), LocAccEIce(ICell), + HTend, TTend, STend); + } + // else { + // temporary kernel logging - TBRemoved + // LOG_INFO("warm layer but no ice to melt"); + // } } - // temporary log -- TBRemoved - // if (ICell == 0) { - // LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, PIn = " - // "{}, H = {}, Tfrz = {}", - // ICell, SAIn, CTIn, PIn, H, Tfrz); - // LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " - // "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", - // ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), - // LocAccMLiq(ICell), LocAccMSalt(ICell), - // LocAccELiq(ICell), LocAccEIce(ICell)); - // LOG_INFO(" HTend= {} " - // "TTend= {} STend= {}", - // HTend, TTend, STend); + // // temporary kernel logging -- TBRemoved + // if (ICell == 0) { + // const Real Hf = H + HTend; + // const Real SAf = (H * SAIn + STend) / Hf; + // const Real Tf = (H * CTIn + TTend) / Hf; + // // LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, + // PIn + // // = " + // // "{}, H = {}, Tfrz = {}", + // // ICell, SAIn, CTIn, PIn, H, Tfrz); + // LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} " + // "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}", + // ICell, K, (CTIn < Tfrz), LocAccMIce(ICell), + // LocAccMLiq(ICell), LocAccMSalt(ICell), + // LocAccELiq(ICell), LocAccEIce(ICell)); + // LOG_INFO(" HTend= {} " + // "TTend= {} STend= {}", + // HTend, TTend, STend); + // LOG_INFO(" " + // "H= {}-->{}, T= {}-->{}, S= {}-->{}", + // H, Hf, CTIn, Tf, SAIn, SAf); // } LocFrazilHTend(ICell, K) = HTend; // not scaled by dt LocFrazilTTend(ICell, K) = TTend; LocFrazilSTend(ICell, K) = STend; -} -// Convert to coupler units -LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; -LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; -LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; -LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; -LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; -}); - -if (conservationCheck) { - checkColumnConservation(); -} -} + } + + // Convert to coupler units + LocAccMIce(ICell) = LocAccMIce(ICell) * RhoSw; + LocAccMLiq(ICell) = LocAccMLiq(ICell) * RhoSw; + LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt; + LocAccELiq(ICell) = LocAccELiq(ICell) * RhoSw; + LocAccEIce(ICell) = LocAccEIce(ICell) * RhoSw; + }); // end parallelFor + + if (conservationCheck) { + checkColumnConservation(); + } +} // end of computeFrazil } // namespace OMEGA diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h index 1121338ef284..bf99f3caedd5 100644 --- a/components/omega/src/ocn/Frazil.h +++ b/components/omega/src/ocn/Frazil.h @@ -2,8 +2,9 @@ #define OMEGA_FRAZIL_H //===-- ocn/Frazil.h - Frazil Ice Formation -------------------*- C++ -*-===// // -// This header defines a scaffold for frazil-related tendencies and -// accumulators. Physics implementations are intentionally left empty. +// The Frazil class manages frazil tendencies and accumulators. +// This initial implementation only has a teos-10 configuration. +// but carries scaffolding for other implementations. // //===----------------------------------------------------------------------===// @@ -33,10 +34,14 @@ class FrazilMelt { /// constructor declaration FrazilMelt(); - // The functor takes the full arrays of specific volume (inout), - // the indices ICell and KChunk, and the ocean tracers (conservative) - // temperature, and (absolute) salinity as inputs, and outputs the - // specific volume according to the Roquet et al. 2015 75 term expansion. + // masslimit parameter (set in config) + Real massLimit; + + // The functor for FrazilMelt takes as inputs: + // the local ocean layer state (SA, CT, P, h), + // the accumulated frazil solid and liquid mass, energy, and salt + // and outputs the frazil tendencies (HTend, TTend, STend) and updated + // accumulators. KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P, const Real h, Real &AccMIce, Real &AccMLiq, Real &AccMSalt, Real &AccELiq, Real &AccEIce, @@ -45,36 +50,45 @@ class FrazilMelt { constexpr Real Eps = 1.0e-12_Real; - if (AccMIce <= Eps) { // potential leak if we dont redistribute + // this check on AccMIce and E should be done in the calling function, but + // is here for safety we can do a better implementation of the checks + if (AccMIce <= Eps || + AccEIce >= Eps) { // potential leak if we dont redistribute AccMIce = 0.0_Real; + AccEIce = 0.0_Real; HTend = 0.0_Real; TTend = 0.0_Real; STend = 0.0_Real; return; } - const Real safeAccMLiq = Kokkos::max(AccMLiq, Eps); - const Real safeAccMIce = Kokkos::max(AccMIce, Eps); - const Real frazilIceFraction = Kokkos::min( - 1.0_Real, - Kokkos::max(0.0_Real, AccMIce / (safeAccMLiq + safeAccMIce))); - const Real brineSalinity = AccMSalt / safeAccMLiq; - const Real brineEnthalpy = AccELiq / safeAccMLiq; - const Real potEnthalpyIce = AccEIce / safeAccMIce; - - const Real layerMass = h + AccMIce; - const Real safeLayerMass = Kokkos::max(layerMass, Eps); - - const Real layerIceFraction = - Kokkos::min(1.0_Real, Kokkos::max(0.0_Real, AccMIce / safeLayerMass)); - - // typecasting for now but will be simplified - const double SA_d = static_cast(SA); - const double CT_d = static_cast(CT); - const double P_d = static_cast(P); - const double wIhIn_d = static_cast(layerIceFraction); - const double pt0Ice_d = - gsw_pt_from_pot_enthalpy_ice(static_cast(potEnthalpyIce)); + if (AccMLiq <= Eps || + AccELiq >= Eps) { // potential leak if we dont redistribute + AccMLiq = 0.0_Real; + AccELiq = 0.0_Real; + HTend = 0.0_Real; + TTend = 0.0_Real; + STend = 0.0_Real; + return; + } + + // 1. we start by adding the solid ice to the ocean layer (no brine yet) + const Real potEnthalpyIce = AccEIce / AccMIce; + const Real newLayerMass = + Kokkos::max(h + AccMIce, Eps); // max unnecessary but for safety + const Real newLayerIceFraction = AccMIce / newLayerMass; + + // 2. we calculate the (mass- and energy-conserving) ocean layer evolution + // ... how much (pure) ice does this layer melt? + + // typecasting for now but will be simplified once gsw functions are + // ported + const double SA_d = static_cast(SA); + const double CT_d = static_cast(CT); + const double P_d = static_cast(P); + const double wIhIn_d = static_cast(newLayerIceFraction); + const double pt0Ice_d = gsw_pt_from_pot_enthalpy_ice_poly( + static_cast(potEnthalpyIce)); const double tIce_d = gsw_t_from_pt0_ice(pt0Ice_d, P_d); double SAnew_d = SA_d; double CTnew_d = CT_d; @@ -83,38 +97,32 @@ class FrazilMelt { gsw_melting_ice_into_seawater(SA_d, CT_d, P_d, wIhIn_d, tIce_d, &SAnew_d, &CTnew_d, &wIhOut_d); - const Real wIhOut = Kokkos::min( - 1.0_Real, Kokkos::max(0.0_Real, static_cast(wIhOut_d))); - - const Real finalSolidMass = - Kokkos::min(AccMIce, Kokkos::max(0.0_Real, wIhOut * safeLayerMass)); - const Real solidMass = Kokkos::max(0.0_Real, AccMIce - finalSolidMass); - - if (solidMass <= Eps) { - return; - } - - const Real liquidMass = - Kokkos::min(AccMLiq, solidMass * (1.0_Real - frazilIceFraction) / - Kokkos::max(frazilIceFraction, Eps)); - const Real solidEnthalpy = solidMass * potEnthalpyIce; - const Real liquidEnthalpy = liquidMass * brineEnthalpy; - - HTend = +(solidMass + liquidMass); - TTend = +(liquidEnthalpy + solidEnthalpy) / Cp0Sw; - STend = +(liquidMass * brineSalinity); - - AccMIce = Kokkos::max(0.0_Real, AccMIce - solidMass); - AccMLiq = Kokkos::max(0.0_Real, AccMLiq - liquidMass); - AccMSalt = Kokkos::max(0.0_Real, AccMSalt - liquidMass * brineSalinity); - AccELiq -= liquidEnthalpy; - AccEIce -= solidEnthalpy; - if (AccMIce <= Eps) { - AccEIce = 0.0_Real; - } - if (AccMLiq <= Eps) { - AccELiq = 0.0_Real; - } + const Real wIhOut = static_cast( + wIhOut_d); // by def 0<= wIhOut <= 1 ; To-do: check function behavior + + // 3. we calculate the mass fraction of the frazil (pure) ice that was + // melted - limited by a total mass limit of 0.1h + const Real solidMassMelted = Kokkos::max( + 0.0_Real, + AccMIce - wIhOut * newLayerMass); // original - left-over solid ice, + const Real frazilFractionMelted = Kokkos::min( + solidMassMelted / AccMIce, + h * massLimit / (AccMIce + AccMLiq)); // total added mass < 0.1h + + // the frazil fraction based on the solid ice also sets the (proportional) + // contributions from the frazil brine + HTend = +(frazilFractionMelted * (AccMLiq + AccMIce)); + TTend = +(frazilFractionMelted * (AccELiq + AccEIce)) / Cp0Sw; + STend = +(frazilFractionMelted * AccMSalt); + + const Real frazilFractionLeft = + Kokkos::max(0.0_Real, 1.0_Real - frazilFractionMelted); + + AccMIce = frazilFractionLeft * AccMIce; + AccMLiq = frazilFractionLeft * AccMLiq; + AccMSalt = frazilFractionLeft * AccMSalt; + AccELiq = frazilFractionLeft * AccELiq; + AccEIce = frazilFractionLeft * AccEIce; } }; @@ -123,15 +131,15 @@ class FrazilFormation { /// constructor declaration FrazilFormation(); - /// Parameters for FrazilFormation (overwritten by config file if set there) - Real Phi = 0.75_Real; ///< liquid mass fraction of frazil for export - Real MassLimit = - 0.1_Real; ///< layer mass fraction limit for thickness tendency + /// Parameters for FrazilFormation (set by yaml file) + Real phi; ///< liquid mass fraction of new frazil (0 < Phi < 1) + Real massLimit; ///< layer mass fraction limit for thickness tendency - // The functor takes the full arrays of specific volume (inout), - // the indices ICell and KChunk, and the ocean tracers (conservative) - // temperature, and (absolute) salinity as inputs, and outputs the - // specific volume according to the Roquet et al. 2015 75 term expansion. + // The functor for FrazilFormation takes as inputs: + // the local ocean layer state (SA, CT, P, h), + // the accumulated frazil solid and liquid mass, energy, and salt + // and outputs the frazil tendencies (HTend, TTend, STend) and updated + // accumulators. KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P, const Real h, Real &AccMIce, Real &AccMLiq, Real &AccMSalt, Real &AccELiq, Real &AccEIce, @@ -163,20 +171,16 @@ class FrazilFormation { CTnew = static_cast(CTnew_d); wIh = static_cast(wIh_d); - const Real OneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - Phi); + const Real oneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - phi); // anything called mass below is in pseudo-thickness units (m) and needs // to be scaled by RhoSw for coupling - solidMass = h * Kokkos::min(wIh, OneMinusPhi * MassLimit); - liquidMass = (Phi / OneMinusPhi) * solidMass; + solidMass = h * Kokkos::min(wIh, oneMinusPhi * massLimit); + liquidMass = (phi / oneMinusPhi) * solidMass; solidEnthalpy = solidMass * gsw_pot_enthalpy_from_pt_ice_poly(PTnew_d); liquidEnthalpy = liquidMass * Cp0Sw * CTnew; - // per timestep (not scaled by dt here) HTend = -(solidMass + liquidMass); - // TTend = (h - solidMass - liquidMass) * CTnew - h * CT; - // STend = (h - solidMass - liquidMass) * SAnew - h * SA; - TTend = -(liquidEnthalpy + solidEnthalpy) / - Cp0Sw; // convert back to CT tendency + TTend = -(liquidEnthalpy + solidEnthalpy) / Cp0Sw; STend = -(liquidMass * SAnew); // Local unit of mass is pseudo thickness (m) @@ -237,12 +241,8 @@ class Frazil { Frazil &operator=(Frazil &&) = delete; FrazilType frazilChoice; - BasicFrazilFormation computeBasicFrazilFormation; - BasicFrazilMelt computeBasicFrazilMelt; FrazilFormation computeFrazilFormation; FrazilMelt computeFrazilMelt; - Real massLimit; - Real phi; I4 NCellsAll; I4 NChunks; diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp index ddc1b26e6b9a..b118b0f12d02 100644 --- a/components/omega/test/ocn/FrazilTest.cpp +++ b/components/omega/test/ocn/FrazilTest.cpp @@ -88,6 +88,8 @@ void testFrazilFormationCold() { (void)Mesh; FrazilFormation ComputeFrazilFormation; + ComputeFrazilFormation.phi = 0.75_Real; + ComputeFrazilFormation.massLimit = 0.1_Real; Real AccMIce = 0.0_Real; Real AccMLiq = 0.0_Real; @@ -166,6 +168,8 @@ void testFrazilFormationWarm() { (void)Mesh; FrazilFormation ComputeFrazilFormation; + ComputeFrazilFormation.phi = 0.75_Real; + ComputeFrazilFormation.massLimit = 0.1_Real; Real AccMIce = 0.0_Real; Real AccMLiq = 0.0_Real; diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5df40c343a2c..e0b151cf4b88 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -439,7 +439,7 @@ int testTendencies() { LOG_ERROR("TendenciesTest: Frazil did not change Salinity tracer " "tendency"); } - + // add a log info if no errors? DefTendencies->FrazilTerm.Enabled = OrigFrazilEnabled; } From 090bc04517a4bf8449f0b21a05295b355dbab156 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Thu, 9 Jul 2026 18:16:55 -0700 Subject: [PATCH 11/11] clean up initializer order to remove build warning --- components/omega/src/ocn/Frazil.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp index 9d0f1ee0b386..9a323f51b76f 100644 --- a/components/omega/src/ocn/Frazil.cpp +++ b/components/omega/src/ocn/Frazil.cpp @@ -71,9 +71,10 @@ void Frazil::init() { } Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord) - : frazilChoice(FrazilType::TeosFrazil), NCellsAll(Mesh->NCellsAll), + : frazilChoice(FrazilType::TeosFrazil), computeFrazilFormation(), + computeFrazilMelt(), NCellsAll(Mesh->NCellsAll), NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh), - VCoordPtr(VCoord), computeFrazilFormation(), computeFrazilMelt() { + VCoordPtr(VCoord) { FrazilTTend = Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers);