From 25f32a572e743f12d75555fd9ee3e755d9375a19 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 16:06:23 -0400 Subject: [PATCH 01/24] Add array type structs --- components/omega/src/infra/OmegaKokkos.h | 47 ++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/components/omega/src/infra/OmegaKokkos.h b/components/omega/src/infra/OmegaKokkos.h index 9c5466574e19..aff22787dbda 100644 --- a/components/omega/src/infra/OmegaKokkos.h +++ b/components/omega/src/infra/OmegaKokkos.h @@ -75,6 +75,53 @@ template struct ArrayRank { static constexpr bool Is5D = T::rank == 5; }; +// Array type traits: Map scalar types to concrete OMEGA array types +template struct Array1D; +template<> struct Array1D { using type = Array1DI4; }; +template<> struct Array1D { using type = Array1DI8; }; +template<> struct Array1D { using type = Array1DR4; }; +template<> struct Array1D { using type = Array1DR8; }; + +template struct Array2D; +template<> struct Array2D { using type = Array2DI4; }; +template<> struct Array2D { using type = Array2DI8; }; +template<> struct Array2D { using type = Array2DR4; }; +template<> struct Array2D { using type = Array2DR8; }; + +template struct Array3D; +template<> struct Array3D { using type = Array3DI4; }; +template<> struct Array3D { using type = Array3DI8; }; +template<> struct Array3D { using type = Array3DR4; }; +template<> struct Array3D { using type = Array3DR8; }; + +template struct Array4D; +template<> struct Array4D { using type = Array4DI4; }; +template<> struct Array4D { using type = Array4DI8; }; +template<> struct Array4D { using type = Array4DR4; }; +template<> struct Array4D { using type = Array4DR8; }; + +template struct Array5D; +template<> struct Array5D { using type = Array5DI4; }; +template<> struct Array5D { using type = Array5DI8; }; +template<> struct Array5D { using type = Array5DR4; }; +template<> struct Array5D { using type = Array5DR8; }; + +// Helper aliases for cleaner syntax +template +using Array1D_t = typename Array1D::type; + +template +using Array2D_t = typename Array2D::type; + +template +using Array3D_t = typename Array3D::type; + +template +using Array4D_t = typename Array4D::type; + +template +using Array5D_t = typename Array5D::type; + template auto createHostMirrorCopy(const V &View) -> Kokkos::View { From 13fa0fcd6694f68d24287b885c3fb6f428fec535 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 16:12:52 -0400 Subject: [PATCH 02/24] Make IOStream::create public --- components/omega/src/infra/IOStream.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/components/omega/src/infra/IOStream.h b/components/omega/src/infra/IOStream.h index ba2556867154..3f6208c9a3b5 100644 --- a/components/omega/src/infra/IOStream.h +++ b/components/omega/src/infra/IOStream.h @@ -151,15 +151,6 @@ class IOStream { bool Validated; //---- Private utility functions to support public interfaces - /// Creates a new stream and adds to the list of all streams, based on - /// options in the input model configuration. This routine is called by - /// the IOStreams initialize function. It requires an initialized model - /// clock so that stream alarm can be attached to this clock during creation. - static void create(const std::string &StreamName, ///< [in] name of stream - Config &StreamConfig, ///< [in] stream configuration - Clock *&ModelClock ///< [inout] Omega model clock - ); - /// Read all dimensions from an input file and determine the dimension ID. /// The file must be in data mode. void readAllDims( @@ -277,6 +268,15 @@ class IOStream { /// Overloaded init with no args, helpful for tests where no Clock exists static void init(void); + /// Creates a new stream and adds to the list of all streams, based on + /// options in the input model configuration. This routine is called by + /// the IOStreams initialize function. It requires an initialized model + /// clock so that stream alarm can be attached to this clock during creation. + static void create(const std::string &StreamName, ///< [in] name of stream + Config &StreamConfig, ///< [in] stream configuration + Clock *&ModelClock ///< [inout] Omega model clock + ); + //--------------------------------------------------------------------------- /// Performs a final write of any streams that have the OnShutdown option and /// then removes all streams to clean up. From 233b3c3b31ff069f238266e1bb59253c3a858ea3 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 16:16:19 -0400 Subject: [PATCH 03/24] Add IOStream::getAlarm method --- components/omega/src/infra/IOStream.cpp | 9 +++++++++ components/omega/src/infra/IOStream.h | 5 +++++ 2 files changed, 14 insertions(+) diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index 95e082480ae9..7dfe07df2257 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/components/omega/src/infra/IOStream.cpp @@ -149,6 +149,15 @@ IOStream::getFilename(const std::string &StreamName ///< [in] name of stream return StreamPtr->Filename; } // End getFilename +//------------------------------------------------------------------------------ +// Retrieves a pointer to the Read/Write alarm for the stream with the input +// name +Alarm *IOStream::getAlarm(const std::string &StreamName ///< [in] name of stream +) { + auto StreamPtr = get(StreamName); + return &StreamPtr->MyAlarm; +} + //------------------------------------------------------------------------------ // Changes the filename associated with a stream (eg during unit testing) void IOStream::changeFilename( diff --git a/components/omega/src/infra/IOStream.h b/components/omega/src/infra/IOStream.h index 3f6208c9a3b5..7f5fe7d82c71 100644 --- a/components/omega/src/infra/IOStream.h +++ b/components/omega/src/infra/IOStream.h @@ -300,6 +300,11 @@ class IOStream { getFilename(const std::string &StreamName ///< [in] name of stream ); + //--------------------------------------------------------------------------- + // Retrieve the Alarm for a given stream + static Alarm *getAlarm(const std::string &StreamName ///< [in] name of stream + ); + //--------------------------------------------------------------------------- /// Changes the filename in a stream (typically during unit testing) static void changeFilename( From 2d691e6c415b4f0217d1910cea627d8324e78a7d Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 16:31:31 -0400 Subject: [PATCH 04/24] Add TimeInterval::isDivisibleBy method Adds the TimeInterval::isDivisible by method to the TimeMgr module. Checks if one TimeInterval is evenly divisible by another. Has specific logic to handle calendar-based intervals and variable-length months --- components/omega/src/infra/TimeMgr.cpp | 172 +++++++++++++++++++++++++ components/omega/src/infra/TimeMgr.h | 11 ++ 2 files changed, 183 insertions(+) diff --git a/components/omega/src/infra/TimeMgr.cpp b/components/omega/src/infra/TimeMgr.cpp index 7653ff44b118..952b5f097e65 100644 --- a/components/omega/src/infra/TimeMgr.cpp +++ b/components/omega/src/infra/TimeMgr.cpp @@ -3159,6 +3159,178 @@ bool TimeInterval::isPositive(void) { } // end TimeInterval::isPositive +//------------------------------------------------------------------------------ +// Check if this TimeInterval is evenly divisible by another TimeInterval +bool TimeInterval::isDivisibleBy( + const TimeInterval &Divisor ///< [in] divisor interval +) const { + + if (!Calendar::isDefined()) { + ABORT_ERROR("TimeInterval::isDivisibleBy: No calendar has been defined. " + "Call Calendar::init() before using time intervals."); + } + + // First check for calendar-based relationships that can be handled directly + // without converting to seconds (e.g., years and months) + + // Years and Months: Use calendar relationship (1 year = 12 months) + if (IsCalendar && Divisor.IsCalendar) { + if ((Units == TimeUnits::Years || Units == TimeUnits::Months) && + (Divisor.Units == TimeUnits::Years || + Divisor.Units == TimeUnits::Months)) { + + // Convert both to months for comparison + I8 ThisInMonths = 0; + I8 DivisorInMonths = 0; + + if (Units == TimeUnits::Years) { + ThisInMonths = CalInterval * MONTHS_PER_YEAR; + } else { // Months + ThisInMonths = CalInterval; + } + + if (Divisor.Units == TimeUnits::Years) { + DivisorInMonths = Divisor.CalInterval * MONTHS_PER_YEAR; + } else { // Months + DivisorInMonths = Divisor.CalInterval; + } + + return (ThisInMonths % DivisorInMonths) == 0; + } + } + + // Special handling for months with day-based intervals + // For calendars with variable month lengths (Gregorian, NoLeap, Julian), + // months are only evenly divisible by 1-day intervals or intervals + // evenly divisible into 1 day. + if ((IsCalendar && Units == TimeUnits::Months) || + (Divisor.IsCalendar && Divisor.Units == TimeUnits::Months)) { + + CalendarKind CalKind = Calendar::OmegaCal->getKind(); + bool HasVariableMonths = + (CalKind == CalendarGregorian || CalKind == CalendarNoLeap || + CalKind == CalendarJulian); + + if (HasVariableMonths) { + // For variable-month calendars, check if comparing month with + // day-based or shorter units + bool ComparingMonthWithDayOrShorter = + (IsCalendar && Units == TimeUnits::Months && + (Divisor.IsCalendar && (Divisor.Units == TimeUnits::Days || + Divisor.Units == TimeUnits::Hours || + Divisor.Units == TimeUnits::Minutes || + Divisor.Units == TimeUnits::Seconds) || + !Divisor.IsCalendar)) || + (Divisor.IsCalendar && Divisor.Units == TimeUnits::Months && + (IsCalendar && + (Units == TimeUnits::Days || Units == TimeUnits::Hours || + Units == TimeUnits::Minutes || + Units == TimeUnits::Seconds) || + !IsCalendar)); + + if (ComparingMonthWithDayOrShorter) { + I8 SecsPerDay = Calendar::OmegaCal->SecondsPerDay; + + if (Divisor.IsCalendar && Divisor.Units == TimeUnits::Months) { + // Divisor is months, dividend is day-based or shorter + // For variable-month calendars, day-based intervals cannot be + // evenly divided by months (since months have variable length) + return false; + } else if (Divisor.IsCalendar && Divisor.Units == TimeUnits::Days) { + // Divisor is in days - only 1 day divides all variable-length + // months + return (Divisor.CalInterval == 1); + } else if (Divisor.IsCalendar && + (Divisor.Units == TimeUnits::Hours || + Divisor.Units == TimeUnits::Minutes || + Divisor.Units == TimeUnits::Seconds)) { + // Divisor is calendar-based hours/minutes/seconds - convert to + // seconds + I8 DivisorSec = + Divisor.CalInterval * Calendar::OmegaCal->SecondsPerDay; + if (Divisor.Units == TimeUnits::Hours) { + DivisorSec = Divisor.CalInterval * SECONDS_PER_HOUR; + } else if (Divisor.Units == TimeUnits::Minutes) { + DivisorSec = Divisor.CalInterval * SECONDS_PER_MINUTE; + } else { // Seconds + DivisorSec = Divisor.CalInterval; + } + + // Check if it divides evenly into a day + if (SecsPerDay % DivisorSec == 0) { + return true; + } else { + return false; + } + } else if (!Divisor.IsCalendar) { + // Non-calendar interval in seconds (hours, minutes, seconds) + I8 DivisorSec = Divisor.Interval.getWhole(); + + // Check if this interval divides evenly into a day + if (SecsPerDay % DivisorSec == 0) { + return true; + } else { + return false; + } + } + } + } + } + + // For all other cases, convert both intervals to seconds for comparison + I8 ThisInSeconds = 0; + I8 DivisorInSeconds = 0; + + // Convert this interval to seconds + if (IsCalendar) { + switch (Units) { + case TimeUnits::Years: + ThisInSeconds = CalInterval * Calendar::OmegaCal->SecondsPerYear; + break; + case TimeUnits::Months: + ThisInSeconds = CalInterval * + (Calendar::OmegaCal->SecondsPerYear / MONTHS_PER_YEAR); + break; + case TimeUnits::Days: + ThisInSeconds = CalInterval * Calendar::OmegaCal->SecondsPerDay; + break; + default: + ABORT_ERROR("TimeInterval::isDivisibleBy: Unknown calendar units"); + } + } else { + // Non-calendar interval - get from TimeFrac as seconds + ThisInSeconds = Interval.getWhole(); + } + + // Convert divisor interval to seconds + if (Divisor.IsCalendar) { + switch (Divisor.Units) { + case TimeUnits::Years: + DivisorInSeconds = + Divisor.CalInterval * Calendar::OmegaCal->SecondsPerYear; + break; + case TimeUnits::Months: + DivisorInSeconds = + Divisor.CalInterval * + (Calendar::OmegaCal->SecondsPerYear / MONTHS_PER_YEAR); + break; + case TimeUnits::Days: + DivisorInSeconds = + Divisor.CalInterval * Calendar::OmegaCal->SecondsPerDay; + break; + default: + ABORT_ERROR("TimeInterval::isDivisibleBy: Unknown calendar units"); + } + } else { + // Non-calendar interval - get from TimeFrac as seconds + DivisorInSeconds = Divisor.Interval.getWhole(); + } + + // Check divisibility using integer modulo + return (ThisInSeconds % DivisorInSeconds) == 0; + +} // end TimeInterval::isDivisibleBy + //------------------------------------------------------------------------------ // TimeInstant definitions //------------------------------------------------------------------------------ diff --git a/components/omega/src/infra/TimeMgr.h b/components/omega/src/infra/TimeMgr.h index 1237fee89d23..f17aa3fe4e4c 100644 --- a/components/omega/src/infra/TimeMgr.h +++ b/components/omega/src/infra/TimeMgr.h @@ -404,6 +404,8 @@ class Calendar { I8 &Month, ///< [in,out] calendar month for time to be advanced I8 &Day ///< [in,out] calendar day for time to be advanced ); + + friend class TimeInterval; }; // end class Calendar /// The TimeInterval class represents an interval of time -- the amount of time @@ -571,6 +573,15 @@ class TimeInterval { /// Check whether a time interval is positive bool isPositive(void); + /// Check if this TimeInterval is evenly divisible by an integer multiple + /// of another TimeInterval. Uses calendar relationships when possible + /// (e.g., 1 year = 12 months exactly). For calendars with variable month + /// lengths, months are only divisible by 1-day intervals or intervals + /// divisible into a 1-day interval. Return true if evenly divisible, + /// false otherwise. + bool isDivisibleBy(const TimeInterval &Divisor ///< [in] divisor interval + ) const; + /// commutative multiplication operators need to be defined as /// free functions, and therefore need to be given acces to /// private members of TimeInterval From 50bc708262576d013b05f1039236beb1101e254a Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 16:57:32 -0400 Subject: [PATCH 05/24] Add inline specifiers to Reductions.h functions Mark non-template functions as inline to prevent multiple definition linker errors when Reductions.h is included in multiple files --- components/omega/src/base/Reductions.h | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/components/omega/src/base/Reductions.h b/components/omega/src/base/Reductions.h index 64b6aab4682f..73bfeb591c69 100644 --- a/components/omega/src/base/Reductions.h +++ b/components/omega/src/base/Reductions.h @@ -46,7 +46,7 @@ static MPI_Op MPI_SUMDD; // special MPI operator for reproducible R8 sum // Special sum function to use in custom double precision MPI reduction. // Implements the DD algorithm for reproducible sums as implemented in // -void sumDD(void *InBuffer, void *OutBuffer, int *Len, MPI_Datatype *DataType) { +inline void sumDD(void *InBuffer, void *OutBuffer, int *Len, MPI_Datatype *DataType) { complex *DDa = (complex *)InBuffer; complex *DDb = (complex *)OutBuffer; double E, T1, T2; @@ -60,7 +60,7 @@ void sumDD(void *InBuffer, void *OutBuffer, int *Len, MPI_Datatype *DataType) { } // Local iteration of a DD sum for use in local sums -void sumDDLocal(complex &DDb, //< [inout] local sum and residual +inline void sumDDLocal(complex &DDb, //< [inout] local sum and residual const double &DDa //< [in] local double to add to sum ) { double T1 = DDa + real(DDb); @@ -138,7 +138,7 @@ void getReduceArrayInfo( } /// Copies some array index info to device -void copyReduceInfoToDevice( +inline void copyReduceInfoToDevice( Array1DI4 &DevRange, ///< [out] index range for reduction Array1DI8 &DevStrides, ///< [out] array stride for each dim std::vector &IRange, ///< [in] index range for reduction @@ -157,7 +157,7 @@ void copyReduceInfoToDevice( //------------------------------------------------------------------------------ // Initialize the special DD sum operator for double precision reproducible sums -void globalSumInit() { +inline void globalSumInit() { int Err = MPI_Op_create(&sumDD, 1, &MPI_SUMDD); if (Err == 0) R8SumNotInitialized = false; @@ -168,7 +168,7 @@ void globalSumInit() { /// with the communicator. Return value is sum across all tasks and sums are /// bit reproducible. /// I4 specific interface -I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summed +inline I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { I4 Result; @@ -179,7 +179,7 @@ I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summed } /// I8 specific interface -I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summed +inline I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { I8 Result; @@ -190,7 +190,7 @@ I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summed } /// R4 specific interface -R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summed +inline R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { R8 LocalTmp, GlobalTmp; @@ -203,7 +203,7 @@ R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summed } /// R8 specific interface -R8 globalSum(const R8 &Val, ///< [in] local scalar value to be summed +inline R8 globalSum(const R8 &Val, ///< [in] local scalar value to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { // initialize reproducible MPI_SUMDD operator @@ -228,7 +228,7 @@ R8 globalSum(const R8 &Val, ///< [in] local scalar value to be summed /// Multifield scalar sums - sums a set of scalar values across MPI tasks /// Return value is a vector containg the sum of each scalar across all tasks /// I4 specific interface -std::vector globalSum( +inline std::vector globalSum( const std::vector &Scalars, ///< [in] vector of scalars to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { @@ -242,7 +242,7 @@ std::vector globalSum( } /// I8 specific interface -std::vector globalSum( +inline std::vector globalSum( const std::vector &Scalars, ///< [in] vector of scalars to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { @@ -256,7 +256,7 @@ std::vector globalSum( } /// R4 specific interface -std::vector globalSum( +inline std::vector globalSum( const std::vector &Scalars, ///< [in] vector of scalars to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { @@ -279,7 +279,7 @@ std::vector globalSum( } /// R8 specific interface -std::vector globalSum( +inline std::vector globalSum( const std::vector &Scalars, ///< [in] vector of scalars to be summed const MPI_Comm Comm ///< [in] MPI communicator ) { From a3c761974d0c8d5ff5cc1e98b5915fca1168e788 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 17:45:20 -0400 Subject: [PATCH 06/24] Add globalMaskedSum functions to Reductions.h Implement globalMaskedSum functions that compute global sums of array products where a higher-dimensional value array is multiplied by a lower-dimensional mask/weight array (1D or 2D). --- components/omega/src/base/Reductions.h | 528 ++++++++++++++++++++++++- 1 file changed, 518 insertions(+), 10 deletions(-) diff --git a/components/omega/src/base/Reductions.h b/components/omega/src/base/Reductions.h index 73bfeb591c69..a5e0c1915743 100644 --- a/components/omega/src/base/Reductions.h +++ b/components/omega/src/base/Reductions.h @@ -46,7 +46,8 @@ static MPI_Op MPI_SUMDD; // special MPI operator for reproducible R8 sum // Special sum function to use in custom double precision MPI reduction. // Implements the DD algorithm for reproducible sums as implemented in // -inline void sumDD(void *InBuffer, void *OutBuffer, int *Len, MPI_Datatype *DataType) { +inline void sumDD(void *InBuffer, void *OutBuffer, int *Len, + MPI_Datatype *DataType) { complex *DDa = (complex *)InBuffer; complex *DDb = (complex *)OutBuffer; double E, T1, T2; @@ -61,7 +62,7 @@ inline void sumDD(void *InBuffer, void *OutBuffer, int *Len, MPI_Datatype *DataT // Local iteration of a DD sum for use in local sums inline void sumDDLocal(complex &DDb, //< [inout] local sum and residual - const double &DDa //< [in] local double to add to sum + const double &DDa //< [in] local double to add to sum ) { double T1 = DDa + real(DDb); double E = T1 - DDa; @@ -168,8 +169,8 @@ inline void globalSumInit() { /// with the communicator. Return value is sum across all tasks and sums are /// bit reproducible. /// I4 specific interface -inline I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summed - const MPI_Comm Comm ///< [in] MPI communicator +inline I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summed + const MPI_Comm Comm ///< [in] MPI communicator ) { I4 Result; int Err = MPI_Allreduce(&Val, &Result, 1, MPI_INT32_T, MPI_SUM, Comm); @@ -179,8 +180,8 @@ inline I4 globalSum(const I4 &Val, ///< [in] local scalar value to be summe } /// I8 specific interface -inline I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summed - const MPI_Comm Comm ///< [in] MPI communicator +inline I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summed + const MPI_Comm Comm ///< [in] MPI communicator ) { I8 Result; int Err = MPI_Allreduce(&Val, &Result, 1, MPI_INT64_T, MPI_SUM, Comm); @@ -190,8 +191,8 @@ inline I8 globalSum(const I8 &Val, ///< [in] local scalar value to be summe } /// R4 specific interface -inline R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summed - const MPI_Comm Comm ///< [in] MPI communicator +inline R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summed + const MPI_Comm Comm ///< [in] MPI communicator ) { R8 LocalTmp, GlobalTmp; LocalTmp = Val; // convert to double for reproducibility @@ -203,8 +204,8 @@ inline R4 globalSum(const R4 &Val, ///< [in] local scalar value to be summe } /// R8 specific interface -inline R8 globalSum(const R8 &Val, ///< [in] local scalar value to be summed - const MPI_Comm Comm ///< [in] MPI communicator +inline R8 globalSum(const R8 &Val, ///< [in] local scalar value to be summed + const MPI_Comm Comm ///< [in] MPI communicator ) { // initialize reproducible MPI_SUMDD operator if (R8SumNotInitialized) @@ -1607,6 +1608,513 @@ globalMaxVal(const Kokkos::View Arr1, ///< [in] 1st array in product return GlobalMax; } +///----------------------------------------------------------------------------- +/// Sums a product of two arrays (eg a mask), where the two arrays need not +/// be the same size and type. The Mask array (Arr2) is either a 1D array (based +/// on the horizontal index) or a 2D array (1st index: horiz, 2nd index: +/// vertical) An optional index range can be specified through a vector of min, +/// max indices for each dimension of the arrays. +/// I4 specific interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I4> +globalMaskedSum(const Kokkos::View Arr1, + const Kokkos::View Arr2, const MPI_Comm Comm, + const std::vector *IndxRange = nullptr) { + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "globalMaskedSum: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "globalMaskedSum: Arr1 rank must be >= Arr2 rank"); + // Verify dimension matching + if (Arr2.rank == 1) { + int HorizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(HorizDim1) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE( + Arr1.rank >= 2, + "globalMaskedSum: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "globalMaskedSum: Vertical dimensions must match"); + } + + // Get array and layout information + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, + IndxRange); ///< [out] true if host array + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + // Compute local sum on host or device + I4 LocalSum = 0; + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + LocalSum += Arr1.data()[Addr1] * Arr2.data()[Addr2]; + } + } + } + } + } + } else { // on device + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, I4 &Lsum) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + Lsum += LocArr1.data()[Addr1] * LocArr2.data()[Addr2]; + } + }, + LocalSum); + } // end if onHost + // Compute final sum by adding local sums from each MPI task + I4 GlobalSum; + int Err = + MPI_Allreduce(&LocalSum, &GlobalSum, 1, MPI_INT32_T, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedSum: Error in MPI_Allreduce"); + return GlobalSum; +} + +/// I8 specific interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I8> +globalMaskedSum(const Kokkos::View Arr1, + const Kokkos::View Arr2, const MPI_Comm Comm, + const std::vector *IndxRange = nullptr) { + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "globalMaskedSum: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "globalMaskedSum: Arr1 rank must be >= Arr2 rank"); + // Verify dimension matching + if (Arr2.rank == 1) { + int HorizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(HorizDim1) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE( + Arr1.rank >= 2, + "globalMaskedSum: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "globalMaskedSum: Vertical dimensions must match"); + } + + // Get array and layout information + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, + IndxRange); ///< [out] true if host array + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + // Compute local sum on host or device + I8 LocalSum = 0; + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + LocalSum += Arr1.data()[Addr1] * Arr2.data()[Addr2]; + } + } + } + } + } + } else { // on device + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, I8 &Lsum) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + Lsum += LocArr1.data()[Addr1] * LocArr2.data()[Addr2]; + } + }, + LocalSum); + } // end if onHost + // Compute final sum by adding local sums from each MPI task + I8 GlobalSum; + int Err = + MPI_Allreduce(&LocalSum, &GlobalSum, 1, MPI_INT64_T, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedSum: Error in MPI_Allreduce"); + return GlobalSum; +} + +/// R4 specific interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R4> +globalMaskedSum(const Kokkos::View Arr1, + const Kokkos::View Arr2, const MPI_Comm Comm, + const std::vector *IndxRange = nullptr) { + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "globalMaskedSum: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "globalMaskedSum: Arr1 rank must be >= Arr2 rank"); + // Verify dimension matching + if (Arr2.rank == 1) { + int HorizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(HorizDim1) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE( + Arr1.rank >= 2, + "globalMaskedSum: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "globalMaskedSum: Vertical dimensions must match"); + } + + // Get array and layout information + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, + IndxRange); ///< [out] true if host array + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + // Compute local sum on host or device + R8 LocalSum = 0; + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + // convert each to R8 to be sure prod is computed in R8 + R8 DTmp1 = Arr1.data()[Addr1]; + R8 DTmp2 = Arr2.data()[Addr2]; + LocalSum += DTmp1 * DTmp2; + } + } + } + } + } + } else { // on device + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, R8 &Lsum) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + // convert each to R8 to be sure prod is computed in R8 + R8 DTmp1 = LocArr1.data()[Addr1]; + R8 DTmp2 = static_cast(LocArr2.data()[Addr2]); + Lsum += DTmp1 * DTmp2; + } + }, + LocalSum); + } // end if onHost + // Compute final sum by adding local sums from each MPI task + R8 GlobalTmp; + int Err = MPI_Allreduce(&LocalSum, &GlobalTmp, 1, MPI_DOUBLE, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedSum: Error in MPI_Allreduce"); + R4 GlobalSum = GlobalTmp; + return GlobalSum; +} + +/// R8 specific interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R8> +globalMaskedSum(const Kokkos::View Arr1, + const Kokkos::View Arr2, const MPI_Comm Comm, + const std::vector *IndxRange = nullptr) { + + // initialize reproducible MPI_SUMDD operator + if (R8SumNotInitialized) + globalSumInit(); + + using Scalar2 = typename Kokkos::View::value_type; + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "globalMaskedSum: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "globalMaskedSum: Arr1 rank must be >= Arr2 rank"); + + // Verify dimension matching + if (Arr2.rank == 1) { + int HorizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(HorizDim1) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE( + Arr1.rank >= 2, + "globalMaskedSum: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "globalMaskedSum: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "globalMaskedSum: Vertical dimensions must match"); + } + + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, IndxRange); + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + complex DDTmp(0.0, 0.0); + + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + R8 ProdTmp = Arr1.data()[Addr1] * + static_cast(Arr2.data()[Addr2]); + sumDDLocal(DDTmp, ProdTmp); + } + } + } + } + } + } else { + R8 LocalSum = 0.0; + + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, R8 &Lsum) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + Lsum += LocArr1.data()[Addr1] * + static_cast(LocArr2.data()[Addr2]); + } + }, + LocalSum); + + DDTmp = complex(LocalSum, 0.0); + } + + // Compute final sum by adding local sums from each MPI task + complex GlobalTmp(0.0, 0.0); + int Err = MPI_Allreduce(&DDTmp, &GlobalTmp, 1, MPI_C_DOUBLE_COMPLEX, + MPI_SUMDD, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedSum: Error in MPI_Allreduce"); + + R8 GlobalSum = real(GlobalTmp); + return GlobalSum; +} //------------------------------------------------------------------------------ } // end namespace OMEGA From c83f9af03a96495578726e6e2b225a8e3d6fbeaf Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 18:05:20 -0400 Subject: [PATCH 07/24] Add globalMaskedMin/Max functions to Reductions.h Implement globalMaskedMin and globalMaskedMax funcitons that compute global minimum and maximum values of array products where a higher-dimensional value array is multiplied by a lower-dimensional mask/weight array (1D or 2D). --- components/omega/src/base/Reductions.h | 478 ++++++++++++++++++++++++- 1 file changed, 476 insertions(+), 2 deletions(-) diff --git a/components/omega/src/base/Reductions.h b/components/omega/src/base/Reductions.h index a5e0c1915743..65c31a71bb9c 100644 --- a/components/omega/src/base/Reductions.h +++ b/components/omega/src/base/Reductions.h @@ -1610,8 +1610,8 @@ globalMaxVal(const Kokkos::View Arr1, ///< [in] 1st array in product } ///----------------------------------------------------------------------------- /// Sums a product of two arrays (eg a mask), where the two arrays need not -/// be the same size and type. The Mask array (Arr2) is either a 1D array (based -/// on the horizontal index) or a 2D array (1st index: horiz, 2nd index: +/// be the same shape and type. The Mask array (Arr2) is either a 1D array +/// (based on the horizontal index) or a 2D array (1st index: horiz, 2nd index: /// vertical) An optional index range can be specified through a vector of min, /// max indices for each dimension of the arrays. /// I4 specific interface @@ -2116,6 +2116,480 @@ globalMaskedSum(const Kokkos::View Arr1, return GlobalSum; } +///----------------------------------------------------------------------------- +/// Finds the local minimum of a Kokkos array product where the arrays need not +/// be the same shape and type. +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + void> +localMaskedMin( + const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask/weight (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + IT &LocalMinVal, ///< [out] local min + const std::vector *IndxRange = nullptr ///< [in] index range +) { + using Scalar2 = typename Kokkos::View::value_type; + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "localMaskedMin: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "localMaskedMin: Arr1 rank must be >= Arr2 rank"); + + // Verify dimension matching + if (Arr2.rank == 1) { + int horizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(horizDim1) == Arr2.extent(0), + "localMaskedMin: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE(Arr1.rank >= 2, + "localMaskedMin: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "localMaskedMin: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "localMaskedMin: Vertical dimensions must match"); + } + + // Get array and layout information + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, IndxRange); + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + // Compute local minimum on host or device + IT LocalMin = std::numeric_limits::max(); + + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + R8 MaskVal = static_cast(Arr2.data()[Addr2]); + if (MaskVal != 0.0) { + IT TestVal = Arr1.data()[Addr1] * MaskVal; + if (TestVal < LocalMin) + LocalMin = TestVal; + } + } + } + } + } + } + } else { // on device + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, IT &Lmin) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + R8 MaskVal = static_cast(LocArr2.data()[Addr2]); + if (MaskVal != 0.0) { + IT TestVal = LocArr1.data()[Addr1] * MaskVal; + Lmin = Kokkos::min(TestVal, Lmin); + } + } + }, + Kokkos::Min(LocalMin)); + } + + LocalMinVal = LocalMin; + return; +} + +///----------------------------------------------------------------------------- +/// Global minimum of Kokkos array product with broadcasting (R8 interface) +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R8> +globalMaskedMin( + const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask/weight (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local min (also performs error checks) + R8 LocalMin; + localMaskedMin(Arr1, Arr2, Comm, LocalMin, IndxRange); + + // Compute final minimum across MPI tasks + R8 GlobalMin; + int Err = MPI_Allreduce(&LocalMin, &GlobalMin, 1, MPI_DOUBLE, MPI_MIN, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMin: Error in MPI_Allreduce"); + + return GlobalMin; +} + +/// Specific R4 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R4> +globalMaskedMin( + const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask/weight (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local min (also performs error checks) + R4 LocalMin; + localMaskedMin(Arr1, Arr2, Comm, LocalMin, IndxRange); + + // Compute final minimum across MPI tasks + R4 GlobalMin; + int Err = MPI_Allreduce(&LocalMin, &GlobalMin, 1, MPI_FLOAT, MPI_MIN, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMin: Error in MPI_Allreduce"); + + return GlobalMin; +} +/// Specific I4 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I4> +globalMaskedMin( + const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask/weight (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local min (also performs error checks) + I4 LocalMin; + localMaskedMin(Arr1, Arr2, Comm, LocalMin, IndxRange); + + // Compute final minimum across MPI tasks + I4 GlobalMin; + int Err = + MPI_Allreduce(&LocalMin, &GlobalMin, 1, MPI_INT32_T, MPI_MIN, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMin: Error in MPI_Allreduce"); + + return GlobalMin; +} + +/// Specific I8 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I8> +globalMaskedMin( + const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask/weight (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local min (also performs error checks) + I8 LocalMin; + localMaskedMin(Arr1, Arr2, Comm, LocalMin, IndxRange); + + // Compute final minimum across MPI tasks + I8 GlobalMin; + int Err = + MPI_Allreduce(&LocalMin, &GlobalMin, 1, MPI_INT64_T, MPI_MIN, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMin: Error in MPI_Allreduce"); + + return GlobalMin; +} + +///----------------------------------------------------------------------------- +/// Finds the local maximum of a Kokkos array product where the arrays need not +/// be the same shape and type. +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + void> +localMaskedMax(const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + IT &LocalMaxVal, ///< [out] local max + const std::vector *IndxRange = nullptr ///< [in] index range +) { + using Scalar2 = typename Kokkos::View::value_type; + + // Error checks + OMEGA_REQUIRE(Arr2.rank == 1 || Arr2.rank == 2, + "localMaskedMin: Arr2 must be 1D or 2D"); + OMEGA_REQUIRE(Arr1.rank >= Arr2.rank, + "localMaskedMin: Arr1 rank must be >= Arr2 rank"); + + // Verify dimension matching + if (Arr2.rank == 1) { + int horizDim1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + OMEGA_REQUIRE(Arr1.extent(horizDim1) == Arr2.extent(0), + "localMaskedMin: Horizontal dimensions must match"); + } else { // Arr2.rank == 2 + OMEGA_REQUIRE(Arr1.rank >= 2, + "localMaskedMin: Arr1 must be at least 2D when Arr2 is 2D"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 2) == Arr2.extent(0), + "localMaskedMin: Horizontal dimensions must match"); + OMEGA_REQUIRE(Arr1.extent(Arr1.rank - 1) == Arr2.extent(1), + "localMaskedMin: Vertical dimensions must match"); + } + + // Get array and layout information + bool IsHost = isReduceArrayOnHost(Arr1); + std::vector Strides1(5, 0); + std::vector IRange(10, 0); + getReduceArrayInfo(IRange, Strides1, Arr1, IndxRange); + + // Compute Arr2 strides + I8 Stride2H = (Arr2.rank == 1) ? 1 : Arr2.extent(1); + I8 Stride2V = 1; + + // Compute local maximum on host or device + IT LocalMax = std::numeric_limits::lowest(); + + if (IsHost) { + for (I4 I = IRange[0]; I <= IRange[1]; ++I) { + for (I4 J = IRange[2]; J <= IRange[3]; ++J) { + for (I4 K = IRange[4]; K <= IRange[5]; ++K) { + for (I4 L = IRange[6]; L <= IRange[7]; ++L) { + for (I4 M = IRange[8]; M <= IRange[9]; ++M) { + + size_t Addr1 = I * Strides1[0] + J * Strides1[1] + + K * Strides1[2] + L * Strides1[3] + + M * Strides1[4]; + + size_t Addr2 = 0; + std::array Indices = {I, J, K, L, M}; + + if (Arr2.rank == 1) { + int HorizIdx1 = (Arr1.rank == 1) ? 0 : Arr1.rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Arr1.rank - 2]; + int VertIdx = Indices[Arr1.rank - 1]; + Addr2 = HorizIdx * Stride2H + VertIdx * Stride2V; + } + + R8 MaskVal = static_cast(Arr2.data()[Addr2]); + if (MaskVal != 0.0) { + IT TestVal = Arr1.data()[Addr1] * MaskVal; + if (TestVal > LocalMax) + LocalMax = TestVal; + } + } + } + } + } + } + } else { // on device + const int Arr1Rank = Arr1.rank; + const int Arr2Rank = Arr2.rank; + + Array1DI4 DevRange("IRange", 10); + Array1DI8 DevStrides("Strides", 5); + copyReduceInfoToDevice(DevRange, DevStrides, IRange, Strides1); + OMEGA_SCOPE(LocStrides, DevStrides); + OMEGA_SCOPE(LocRange, DevRange); + OMEGA_SCOPE(LocArr1, Arr1); + OMEGA_SCOPE(LocArr2, Arr2); + + parallelReduce( + {IRange[1] + 1, IRange[3] + 1, IRange[5] + 1, IRange[7] + 1, + IRange[9] + 1}, + KOKKOS_LAMBDA(int I, int J, int K, int L, int M, IT &Lmax) { + if (I >= LocRange(0) and J >= LocRange(2) and K >= LocRange(4) and + L >= LocRange(6) and M >= LocRange(8)) { + + size_t Addr1 = I * LocStrides(0) + J * LocStrides(1) + + K * LocStrides(2) + L * LocStrides(3) + + M * LocStrides(4); + + size_t Addr2 = 0; + int Indices[5] = {I, J, K, L, M}; + if (Arr2Rank == 1) { + int HorizIdx1 = (Arr1Rank == 1) ? 0 : Arr1Rank - 2; + Addr2 = Indices[HorizIdx1]; + } else { + int HorizIdx = Indices[Kokkos::max(0, Arr1Rank - 2)]; + int VertIdx = Indices[Arr1Rank - 1]; + Addr2 = HorizIdx * LocArr2.extent(1) + VertIdx; + } + + R8 MaskVal = static_cast(LocArr2.data()[Addr2]); + if (MaskVal != 0.0) { + IT TestVal = LocArr1.data()[Addr1] * MaskVal; + Lmax = Kokkos::max(TestVal, Lmax); + } + } + }, + Kokkos::Max(LocalMax)); + } + + LocalMaxVal = LocalMax; + return; +} + +///----------------------------------------------------------------------------- +/// Global maximum of masked Kokkos array product (R8 interface) +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R8> +globalMaskedMax(const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local max (also performs error checks) + R8 LocalMax; + localMaskedMax(Arr1, Arr2, Comm, LocalMax, IndxRange); + + // Compute final maximum across MPI tasks + R8 GlobalMax; + int Err = MPI_Allreduce(&LocalMax, &GlobalMax, 1, MPI_DOUBLE, MPI_MAX, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMax: Error in MPI_Allreduce"); + + return GlobalMax; +} + +/// Specific R4 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + R4> +globalMaskedMax(const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local max (also performs error checks) + R4 LocalMax; + localMaskedMax(Arr1, Arr2, Comm, LocalMax, IndxRange); + + // Compute final maximum across MPI tasks + R4 GlobalMax; + int Err = MPI_Allreduce(&LocalMax, &GlobalMax, 1, MPI_FLOAT, MPI_MAX, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMax: Error in MPI_Allreduce"); + + return GlobalMax; +} + +/// Specific I4 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I4> +globalMaskedMax(const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local max (also performs error checks) + I4 LocalMax; + localMaskedMax(Arr1, Arr2, Comm, LocalMax, IndxRange); + + // Compute final maximum across MPI tasks + I4 GlobalMax; + int Err = + MPI_Allreduce(&LocalMax, &GlobalMax, 1, MPI_INT32_T, MPI_MAX, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMax: Error in MPI_Allreduce"); + + return GlobalMax; +} + +/// Specific I8 interface +template +std::enable_if_t< + std::is_same_v::value_type> && + std::is_arithmetic_v::value_type>, + I8> +globalMaskedMax(const Kokkos::View Arr1, ///< [in] 1st array + const Kokkos::View Arr2, ///< [in] mask (1D or 2D) + const MPI_Comm Comm, ///< [in] MPI Communicator + const std::vector *IndxRange = nullptr ///< [in] index range +) { + // Call routine to find local max (also performs error checks) + I8 LocalMax; + localMaskedMax(Arr1, Arr2, Comm, LocalMax, IndxRange); + + // Compute final maximum across MPI tasks + I8 GlobalMax; + int Err = + MPI_Allreduce(&LocalMax, &GlobalMax, 1, MPI_INT64_T, MPI_MAX, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("globalMaskedMax: Error in MPI_Allreduce"); + + return GlobalMax; +} + //------------------------------------------------------------------------------ } // end namespace OMEGA From 3280b142e02ccfc7c98323e2df85be5bc457884c Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Jun 2026 20:10:41 -0400 Subject: [PATCH 08/24] Add Analysis operator factory --- .../omega/src/analysis/AnalysisOpFactory.cpp | 115 +++++++++++ .../omega/src/analysis/AnalysisOpFactory.h | 181 ++++++++++++++++++ 2 files changed, 296 insertions(+) create mode 100644 components/omega/src/analysis/AnalysisOpFactory.cpp create mode 100644 components/omega/src/analysis/AnalysisOpFactory.h diff --git a/components/omega/src/analysis/AnalysisOpFactory.cpp b/components/omega/src/analysis/AnalysisOpFactory.cpp new file mode 100644 index 000000000000..a3675997fbd7 --- /dev/null +++ b/components/omega/src/analysis/AnalysisOpFactory.cpp @@ -0,0 +1,115 @@ +//===-- analysis/AnalysisOpFactory.cpp - Factory implementation --*- C++-*-===// +// +// Implementation of AnalysisOpFactory static methods for operator registration +// and creation. The factory maintains a runtime registry (Meyer's singleton) +// that maps fully-qualified operator type keys to constructor functions. +// Registration occurs during program initialization via static initializers. +// Creation performs type-safe dispatch by inspecting Field metadata to select +// the appropriate templated operator specialization. +// +//===----------------------------------------------------------------------===// + +#include "AnalysisOpFactory.h" +#include + +namespace OMEGA { + +//------------------------------------------------------------------------------ +// Registers an operator variant by adding it to the factory registry. Checks +// for duplicate registration and aborts if the key already exists, which +// indicates a programming error (likely duplicate registration in static +// initializers). +void AnalysisOpFactory::registerOperator(const std::string &Label, + CreatorFunc Creator) { + + auto &Reg = registry(); + + // Check for duplicate registration (programming error) + if (Reg.find(Label) != Reg.end()) { + ABORT_ERROR("AnalysisOpFactory: Operator type {} is already registered", + Label); + } + + // Add constructor function to registry + Reg[Label] = Creator; +} // end registerOperator + +//------------------------------------------------------------------------------ +// Creates an operator instance by inspecting the primary upstream Field's +// metadata to construct a fully-qualified type key, looking up the constructor +// in the registry, and invoking it. The primary Field (first in UpstreamNames) +// determines the array type template parameter for the operator. +std::unique_ptr +AnalysisOpFactory::createOp(const std::string &OpType, + const std::vector &UpstreamNames, + Config Options) { + + // Validate that all upstream Fields exist before attempting creation + for (const auto &FieldName : UpstreamNames) { + auto FieldPtr = Field::get(FieldName); + if (!FieldPtr) { + ABORT_ERROR("Field '{}' not found for operator creation", FieldName); + } + } + + // Extract metadata from primary upstream Field (determines array type) + auto FieldPtr = Field::get(UpstreamNames[0]); + ArrayDataType DType = FieldPtr->getType(); + int Rank = FieldPtr->getNumDims(); + ArrayMemLoc MemLoc = FieldPtr->getMemoryLocation(); + + // Map metadata to array type name (e.g., "Array2DR8_Device") + std::string arrayTypeName = getArrayTypeName(DType, Rank, MemLoc); + + // Build fully-qualified operator key (e.g., "SpatialMean_Array2DR8_Device") + std::string FullOpType = OpType + "_" + arrayTypeName; + + // Look up constructor in registry + auto &Reg = registry(); + auto it = Reg.find(FullOpType); + + if (it == Reg.end()) { + ABORT_ERROR("Operator type {} not found", FullOpType); + } + + // Invoke the registered constructor function + return it->second(UpstreamNames, Options); +} // end createOp + +//------------------------------------------------------------------------------ +// Checks whether an operator variant with the given key exists in the registry. +// Returns true if found, false otherwise. +bool AnalysisOpFactory::hasOperator(const std::string &Type) { + auto &Reg = registry(); + return Reg.find(Type) != Reg.end(); +} // end hasOperator + +//------------------------------------------------------------------------------ +// Constructs an array type name string from Field metadata by expanding the +// OMEGA_ANALYSIS_ARRAY_TYPES macro and checking each combination. Returns a +// string like "Array2DR8_Device" that is used as part of the operator variant +// key. Aborts if the metadata combination is not supported. +std::string AnalysisOpFactory::getArrayTypeName(ArrayDataType DType, int Rank, + ArrayMemLoc MemLoc) { +// Define a macro to check if metadata matches this array type +#define TRY_ARRAY_TYPE(dt, r, ml, ArrayT) \ + if (DType == dt && Rank == r && MemLoc == ml) { \ + return std::string(#ArrayT) + "_" + #ml; \ + } + + // Expand macro over all supported array type combinations + OMEGA_ANALYSIS_ARRAY_TYPES(TRY_ARRAY_TYPE) + +#undef TRY_ARRAY_TYPE + + // If we reach here, the array type is not supported + ABORT_ERROR( + "Unsupported array type/Rank/location: DType={}, Rank={}, MemLoc={}", + static_cast(DType), Rank, static_cast(MemLoc)); + + return {}; +} // end getArrayTypeName + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/analysis/AnalysisOpFactory.h b/components/omega/src/analysis/AnalysisOpFactory.h new file mode 100644 index 000000000000..f394f8118e23 --- /dev/null +++ b/components/omega/src/analysis/AnalysisOpFactory.h @@ -0,0 +1,181 @@ +#ifndef OMEGA_ANALYSISOPFACTORY_H +#define OMEGA_ANALYSISOPFACTORY_H + +//===-- analysis/AnalysisOpFactory.h - Operator factory ---------*- C++ -*-===// +// +/// \file +/// \brief Defines the AnalysisOpFactory for runtime operator creation +/// +/// The AnalysisOpFactory provides a runtime registry that maps operator type +/// names to constructor functions, enabling dynamic instantiation of templated +/// operators without hard-coded switch statements. The factory: +/// - Supports decentralized registration via template helpers +/// - Performs type-safe dispatch based on Field metadata (scalar type, rank, +/// memory location) +/// - Enables extensibility by allowing new operators to integrate without +/// modifying orchestration code +/// +/// Each operator template (e.g., SpatialMaxOp) is registered for all +/// supported array type combinations via registerAllArrayVariants(). At +/// operator creation time, the factory inspects the primary upstream Field's +/// metadata and selects the matching templated specialization. +/// +//===----------------------------------------------------------------------===// + +#include "AnalysisOperator.h" +#include "Config.h" + +#include + +namespace OMEGA { + +/// Macro defining all supported array type combinations for analysis operators. +/// Used by registerAllArrayVariants() to register operator variants for all +/// combinations of scalar type, rank, and memory location. The macro takes a +/// single argument X which should be a macro that processes each +/// (DType, Rank, MemLoc, ArrayT) tuple. +#define OMEGA_ANALYSIS_ARRAY_TYPES(X) \ + /* 1D Arrays */ \ + X(ArrayDataType::I4, 1, ArrayMemLoc::Both, Array1DI4) \ + X(ArrayDataType::I4, 1, ArrayMemLoc::Device, Array1DI4) \ + X(ArrayDataType::I8, 1, ArrayMemLoc::Both, Array1DI8) \ + X(ArrayDataType::I8, 1, ArrayMemLoc::Device, Array1DI8) \ + X(ArrayDataType::R4, 1, ArrayMemLoc::Both, Array1DR4) \ + X(ArrayDataType::R4, 1, ArrayMemLoc::Device, Array1DR4) \ + X(ArrayDataType::R8, 1, ArrayMemLoc::Both, Array1DR8) \ + X(ArrayDataType::R8, 1, ArrayMemLoc::Device, Array1DR8) \ + /* 2D Arrays */ \ + X(ArrayDataType::I4, 2, ArrayMemLoc::Both, Array2DI4) \ + X(ArrayDataType::I4, 2, ArrayMemLoc::Device, Array2DI4) \ + X(ArrayDataType::I8, 2, ArrayMemLoc::Both, Array2DI8) \ + X(ArrayDataType::I8, 2, ArrayMemLoc::Device, Array2DI8) \ + X(ArrayDataType::R4, 2, ArrayMemLoc::Both, Array2DR4) \ + X(ArrayDataType::R4, 2, ArrayMemLoc::Device, Array2DR4) \ + X(ArrayDataType::R8, 2, ArrayMemLoc::Both, Array2DR8) \ + X(ArrayDataType::R8, 2, ArrayMemLoc::Device, Array2DR8) \ + /* 3D Arrays */ \ + X(ArrayDataType::I4, 3, ArrayMemLoc::Both, Array3DI4) \ + X(ArrayDataType::I4, 3, ArrayMemLoc::Device, Array3DI4) \ + X(ArrayDataType::I8, 3, ArrayMemLoc::Both, Array3DI8) \ + X(ArrayDataType::I8, 3, ArrayMemLoc::Device, Array3DI8) \ + X(ArrayDataType::R4, 3, ArrayMemLoc::Both, Array3DR4) \ + X(ArrayDataType::R4, 3, ArrayMemLoc::Device, Array3DR4) \ + X(ArrayDataType::R8, 3, ArrayMemLoc::Both, Array3DR8) \ + X(ArrayDataType::R8, 3, ArrayMemLoc::Device, Array3DR8) // \ +// /* 4D Arrays */ \ +// X(ArrayDataType::I4, 4, ArrayMemLoc::Both, Array4DI4) \ +// X(ArrayDataType::I4, 4, ArrayMemLoc::Device, Array4DI4) \ +// X(ArrayDataType::I8, 4, ArrayMemLoc::Both, Array4DI8) \ +// X(ArrayDataType::I8, 4, ArrayMemLoc::Device, Array4DI8) \ +// X(ArrayDataType::R4, 4, ArrayMemLoc::Both, Array4DR4) \ +// X(ArrayDataType::R4, 4, ArrayMemLoc::Device, Array4DR4) \ +// X(ArrayDataType::R8, 4, ArrayMemLoc::Both, Array4DR8) \ +// X(ArrayDataType::R8, 4, ArrayMemLoc::Device, Array4DR8) \ +// /* 5D Arrays */ \ +// X(ArrayDataType::I4, 5, ArrayMemLoc::Both, Array5DI4) \ +// X(ArrayDataType::I4, 5, ArrayMemLoc::Device, Array5DI4) \ +// X(ArrayDataType::I8, 5, ArrayMemLoc::Both, Array5DI8) \ +// X(ArrayDataType::I8, 5, ArrayMemLoc::Device, Array5DI8) \ +// X(ArrayDataType::R4, 5, ArrayMemLoc::Both, Array5DR4) \ +// X(ArrayDataType::R4, 5, ArrayMemLoc::Device, Array5DR4) \ +// X(ArrayDataType::R8, 5, ArrayMemLoc::Both, Array5DR8) \ +// X(ArrayDataType::R8, 5, ArrayMemLoc::Device, Array5DR8) \ + +/// AnalysisOpFactory is a singleton factory class that manages runtime +/// registration and creation of analysis operators. The factory maintains a +/// registry mapping operator type keys to constructor functions. Keys are +/// formed from base operator name, array type, and memory location (e.g., +/// "SpatialMax_Array2DR8_Device"). This enables type-safe dispatch: the +/// factory inspects the primary upstream Field's metadata at creation time +/// and selects the correct templated specialization. +/// +/// All methods are static; the underlying registry is a Meyer's singleton +/// guaranteed to be initialized before first use. Operators self-register +/// during program initialization via registerAllArrayVariants(). +class AnalysisOpFactory { + public: + /// Function signature for operator constructors. Takes upstream field names + /// and configuration options, returns a unique_ptr to a new operator + /// instance. + using CreatorFunc = std::function( + const std::vector &UpstreamNames, Config Options)>; + + /// Registers a single operator variant in the factory by associating a + /// string label (key) with a constructor function. The label typically + /// includes the operator type, array type, and memory location. Aborts if + /// the label is already registered (duplicate key). + static void registerOperator( + const std::string &Label, ///< [in] unique key for this operator variant + CreatorFunc Creator ///< [in] constructor function for this variant + ); + + /// Creates an operator instance by inspecting the primary upstream Field's + /// metadata (scalar type, rank, memory location) to build a fully-qualified + /// type key, looking up the constructor in the registry, and invoking it + /// with the provided arguments. Aborts if the operator type or array variant + /// is not registered. + static std::unique_ptr + createOp(const std::string &OpType, ///< [in] operator type name + const std::vector + &UpstreamNames, ///< [in] upstream field names + Config Options ///< [in] operator configuration + ); + + /// Returns a list of all registered operator variant keys. Useful for + /// validation and generating informative error messages when an operator + /// type is not found. + static std::vector availableOperators(); + + /// Checks whether an operator variant with the given key is registered. + /// Returns true if found, false otherwise. + static bool + hasOperator(const std::string &Type ///< [in] operator variant key + ); + + /// Registers all array type variants of a templated operator class by + /// expanding OMEGA_ANALYSIS_ARRAY_TYPES over all supported (DType, Rank, + /// MemLoc, ArrayT) combinations. For each combination, constructs a key + /// from baseName + "_" + ArrayT + "_" + MemLoc and registers a lambda + /// that instantiates OperatorTemplate. This enables a single call + /// to register dozens of operator variants. + template