From 8e412a393c6ef22154130a2362ec69e76f20e8ff Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 13 Sep 2026 02:19:18 -0500 Subject: [PATCH 1/5] Add CFUnits, a units algebra for CF units strings Fields derived from other fields need units derived from theirs, and an area-weighted sum or a product cannot simply copy them. CFUnits parses a units string in the CF plain form (m s-1, kg m-3, 1), or the udunits caret and slash spellings still used by some field definitions, into unit symbols with integer exponents; supports the product, quotient and integer power; and formats the result in the plain form with merged exponents, positive exponents first. An empty string is the unknown units of a field with no units attribute, and unknown units propagate through every operation so a derived field never claims units its inputs did not have. A malformed string is a Fail error from parse(), and the string interfaces multiply(), divide() and power() abort on it, since a bad units string is a programming error in a field definition that should stop the run at initialization. Includes a CTest covering parsing, formatting, the algebra, unknown propagation and rejected strings, and a developer-guide page. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/doc/devGuide/CFUnits.md | 46 ++++ components/omega/doc/index.md | 1 + components/omega/src/infra/CFUnits.cpp | 268 ++++++++++++++++++++ components/omega/src/infra/CFUnits.h | 124 +++++++++ components/omega/test/CMakeLists.txt | 11 + components/omega/test/infra/CFUnitsTest.cpp | 204 +++++++++++++++ 6 files changed, 654 insertions(+) create mode 100644 components/omega/doc/devGuide/CFUnits.md create mode 100644 components/omega/src/infra/CFUnits.cpp create mode 100644 components/omega/src/infra/CFUnits.h create mode 100644 components/omega/test/infra/CFUnitsTest.cpp diff --git a/components/omega/doc/devGuide/CFUnits.md b/components/omega/doc/devGuide/CFUnits.md new file mode 100644 index 000000000000..f550619cf7a5 --- /dev/null +++ b/components/omega/doc/devGuide/CFUnits.md @@ -0,0 +1,46 @@ +(omega-dev-cfunits)= + +## CFUnits + +Omega fields carry their units as a string attribute in the plain form used +by the CF standard name table: unit symbols separated by single spaces, each +with an optional signed integer exponent (`m s-1`, `kg m-3`, `W m-2`), and `1` +for a dimensionless quantity. When a field is derived from others, its units +must be derived too. The `CFUnits` class parses such a string, supports the +product, quotient and integer power of units expressions, and formats the +result back in the plain form. + +The string interfaces cover the common cases: +```c++ +#include "CFUnits.h" + +std::string Flux = CFUnits::multiply("m s-1", "m2"); // "m3 s-1" +std::string Speed = CFUnits::divide("m3 s-1", "m2"); // "m s-1" +std::string Energy = CFUnits::power("m s-1", 2); // "m2 s-2" +``` +Exponents of the same symbol are merged (`m s-1` times `m` is `m2 s-1`) and +a symbol whose exponent reaches zero is dropped, so `kg m-3` times `m3` is +`kg`. In the result, symbols with positive exponents come first, then those +with negative exponents, each in the order they first appeared. + +A field with no units attribute has *unknown* units, represented by the empty +string. Unknown units propagate through every operation, so a field derived +from one without units has no units either, rather than wrong ones. + +The parser also accepts the udunits caret and slash spellings (`m^2`, `m/s`), +which some older field definitions still use, but always emits the plain +form. Any other string is an error: the string interfaces abort, and the +value interface reports it so a caller can decide: +```c++ +CFUnits Units; +Error Err = CFUnits::parse("m s-1", Units); // Fail for a malformed string +CFUnits Area; +CFUnits::parse("m2", Area); +std::string Flux = (Units * Area).str(); // "m3 s-1" +``` +`CFUnits` also provides `operator/`, `pow(int)`, `isUnknown()`, +`isDimensionless()` and an order-insensitive `operator==`. + +A malformed units string is a programming error in a field definition, and +aborting when the derived field is created stops the run at initialization +rather than writing output with a wrong or missing attribute. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index b16840d3bc5c..a7f78ede86c3 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -84,6 +84,7 @@ devGuide/Decomp devGuide/Dimension devGuide/Error devGuide/Field +devGuide/CFUnits devGuide/IO devGuide/IOStreams devGuide/Halo diff --git a/components/omega/src/infra/CFUnits.cpp b/components/omega/src/infra/CFUnits.cpp new file mode 100644 index 000000000000..25d8246cac7d --- /dev/null +++ b/components/omega/src/infra/CFUnits.cpp @@ -0,0 +1,268 @@ +//===-- infra/CFUnits.cpp - units algebra for CF units strings --*- C++ -*-===// +// +// Implementation of the CFUnits class. A units string is split on slashes +// into a numerator and denominators, each of those on whitespace into +// factors, and each factor into a unit symbol and an optional integer +// exponent. Factors are kept in order of first appearance so that the +// formatted result reads like the conventional spelling (kg m-3, not +// m-3 kg). +// +//===----------------------------------------------------------------------===// + +#include "CFUnits.h" +#include "Logging.h" + +#include +#include +#include + +namespace OMEGA { + +//------------------------------------------------------------------------------ +// Default constructor creates unknown units +CFUnits::CFUnits() : Known(false) {} + +//------------------------------------------------------------------------------ +// Parses a units string into its factors +Error CFUnits::parse(const std::string &UnitsStr, // [in] units string + CFUnits &Result // [out] parsed units expression +) { + + Error Err; + + Result = CFUnits(); + + // Trim surrounding whitespace; an empty string is the unknown units of a + // field with no units attribute + auto First = UnitsStr.find_first_not_of(" \t"); + if (First == std::string::npos) + return Err; + auto Last = UnitsStr.find_last_not_of(" \t"); + std::string Trimmed = UnitsStr.substr(First, Last - First + 1); + + Result.Known = true; + + // A trailing slash has no denominator after it; getline would skip it + if (Trimmed.back() == '/') { + Result = CFUnits(); + RETURN_ERROR(Err, ErrorCode::Fail, + "CFUnits: cannot parse units string '{}': " + "empty denominator", + UnitsStr); + } + + // Each slash starts a denominator: every factor after it divides + int Sign = 1; + std::string Segment; + std::istringstream Segments(Trimmed); + bool FirstSegment = true; + while (std::getline(Segments, Segment, '/')) { + + // A slash with nothing before or after it is malformed + std::istringstream Tokens(Segment); + std::string Token; + bool AnyToken = false; + while (Tokens >> Token) { + AnyToken = true; + if (!Result.parseFactor(Token, Sign)) { + Result = CFUnits(); + RETURN_ERROR(Err, ErrorCode::Fail, + "CFUnits: cannot parse units string '{}': " + "unexpected factor '{}'", + UnitsStr, Token); + } + } + if (!AnyToken) { + Result = CFUnits(); + RETURN_ERROR(Err, ErrorCode::Fail, + "CFUnits: cannot parse units string '{}': " + "empty {}", + UnitsStr, FirstSegment ? "numerator" : "denominator"); + } + + FirstSegment = false; + Sign = -1; + } + + return Err; + +} // end parse + +//------------------------------------------------------------------------------ +// Parses one factor: a unit symbol of letters and underscores followed by an +// optional integer exponent, with or without a caret (m2, s-1, m^2, s^-1), or +// the digit 1 which contributes nothing +bool CFUnits::parseFactor(const std::string &Token, // [in] one factor + int Sign // [in] +1 or -1 +) { + + if (Token == "1") + return true; + + // Symbol: one or more letters or underscores + std::size_t Pos = 0; + while (Pos < Token.size() && + (std::isalpha(static_cast(Token[Pos])) || + Token[Pos] == '_')) + ++Pos; + if (Pos == 0) + return false; + std::string Symbol = Token.substr(0, Pos); + + // Exponent: optional caret, optional minus sign, then digits to the end + int Exponent = 1; + if (Pos < Token.size()) { + if (Token[Pos] == '^') + ++Pos; + int ExpSign = 1; + if (Pos < Token.size() && Token[Pos] == '-') { + ExpSign = -1; + ++Pos; + } + if (Pos >= Token.size()) + return false; + int Magnitude = 0; + for (; Pos < Token.size(); ++Pos) { + if (!std::isdigit(static_cast(Token[Pos]))) + return false; + Magnitude = 10 * Magnitude + (Token[Pos] - '0'); + } + Exponent = ExpSign * Magnitude; + } + + multiplyFactor(Symbol, Sign * Exponent); + return true; + +} // end parseFactor + +//------------------------------------------------------------------------------ +// Multiplies in one symbol raised to a power +void CFUnits::multiplyFactor(const std::string &Symbol, // [in] unit symbol + int Exponent // [in] its power +) { + + if (Exponent == 0) + return; + + for (auto Iter = Factors.begin(); Iter != Factors.end(); ++Iter) { + if (Iter->first == Symbol) { + Iter->second += Exponent; + if (Iter->second == 0) + Factors.erase(Iter); + return; + } + } + Factors.emplace_back(Symbol, Exponent); + +} // end multiplyFactor + +//------------------------------------------------------------------------------ +// Parses a units string, aborting on failure +CFUnits CFUnits::parseOrAbort(const std::string &UnitsStr // [in] units string +) { + CFUnits Result; + Error Err = parse(UnitsStr, Result); + CHECK_ERROR_ABORT(Err, "CFUnits: invalid units string"); + return Result; +} + +//------------------------------------------------------------------------------ +// String interfaces for the product, quotient and power +std::string CFUnits::multiply(const std::string &UnitsA, // [in] first units + const std::string &UnitsB // [in] second units +) { + return (parseOrAbort(UnitsA) * parseOrAbort(UnitsB)).str(); +} + +std::string CFUnits::divide(const std::string &UnitsA, // [in] numerator + const std::string &UnitsB // [in] denominator +) { + return (parseOrAbort(UnitsA) / parseOrAbort(UnitsB)).str(); +} + +std::string CFUnits::power(const std::string &UnitsA, // [in] units + int Exponent // [in] power +) { + return parseOrAbort(UnitsA).pow(Exponent).str(); +} + +//------------------------------------------------------------------------------ +// Queries +bool CFUnits::isUnknown() const { return !Known; } + +bool CFUnits::isDimensionless() const { return Known && Factors.empty(); } + +//------------------------------------------------------------------------------ +// Formats the units in the CF plain form +std::string CFUnits::str() const { + + if (!Known) + return ""; + if (Factors.empty()) + return "1"; + + std::string Result; + for (int Pass = 0; Pass < 2; ++Pass) { + for (const auto &Factor : Factors) { + bool Positive = Factor.second > 0; + if (Positive != (Pass == 0)) + continue; + if (!Result.empty()) + Result += " "; + Result += Factor.first; + if (Factor.second != 1) + Result += std::to_string(Factor.second); + } + } + return Result; + +} // end str + +//------------------------------------------------------------------------------ +// Algebra +CFUnits CFUnits::operator*(const CFUnits &Other) const { + + if (!Known || !Other.Known) + return CFUnits(); + + CFUnits Result = *this; + for (const auto &Factor : Other.Factors) + Result.multiplyFactor(Factor.first, Factor.second); + return Result; +} + +CFUnits CFUnits::operator/(const CFUnits &Other) const { + return *this * Other.pow(-1); +} + +CFUnits CFUnits::pow(int Exponent) const { + + if (!Known) + return CFUnits(); + + CFUnits Result; + Result.Known = true; + for (const auto &Factor : Factors) + Result.multiplyFactor(Factor.first, Factor.second * Exponent); + return Result; +} + +//------------------------------------------------------------------------------ +// Comparison, ignoring the order of factors +bool CFUnits::operator==(const CFUnits &Other) const { + + if (Known != Other.Known) + return false; + std::map Mine(Factors.begin(), Factors.end()); + std::map Theirs(Other.Factors.begin(), + Other.Factors.end()); + return Mine == Theirs; +} + +bool CFUnits::operator!=(const CFUnits &Other) const { + return !(*this == Other); +} + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/infra/CFUnits.h b/components/omega/src/infra/CFUnits.h new file mode 100644 index 000000000000..cc5b498484a2 --- /dev/null +++ b/components/omega/src/infra/CFUnits.h @@ -0,0 +1,124 @@ +#ifndef OMEGA_CFUNITS_H +#define OMEGA_CFUNITS_H +//===-- infra/CFUnits.h - units algebra for CF units strings ----*- C++ -*-===// +// +/// \file +/// \brief Defines the CFUnits class for deriving units strings +/// +/// Omega fields carry their units as a string attribute in the plain form +/// used by the CF standard name table: unit symbols separated by single +/// spaces, each followed by an optional signed integer exponent (m s-1, +/// kg m-3, W m-2) and the dimensionless quantity written as 1. When a field +/// is derived from others, for example by an area-weighted sum or a product, +/// its units must be derived too. CFUnits parses a units string into its +/// factors, supports the product, quotient and integer power of two such +/// expressions, and formats the result back in the plain form. +/// +/// A field with no units attribute has unknown units, represented by an empty +/// string. Unknown units propagate: any expression involving them is unknown, +/// so a derived field never claims units its inputs did not have. Parsing +/// also accepts the udunits spellings with a caret before the exponent (m^2) +/// and a slash for division (m/s), but always emits the plain form. Anything +/// else is a parse error, since a malformed units string is a programming +/// error in the field definition and should stop the run at initialization. +/// +//===----------------------------------------------------------------------===// + +#include "Error.h" + +#include +#include +#include + +namespace OMEGA { + +/// A units expression as a product of unit symbols raised to integer powers, +/// or the unknown units of a field that has no units attribute +class CFUnits { + + public: + //--------------------------------------------------------------------------- + /// Default constructor creates unknown units + CFUnits(); + + //--------------------------------------------------------------------------- + /// Parses a units string. An empty string gives unknown units and "1" gives + /// dimensionless units. Returns a Fail error if the string is not a product + /// of unit symbols with optional integer exponents, in the CF plain form or + /// the udunits caret and slash spellings. + static Error parse(const std::string &UnitsStr, ///< [in] units string + CFUnits &Result ///< [out] parsed units expression + ); + + //--------------------------------------------------------------------------- + /// Returns the product of two units strings in the CF plain form, or an + /// empty string if either is unknown. Aborts if either string is invalid. + static std::string multiply(const std::string &UnitsA, ///< [in] first units + const std::string &UnitsB ///< [in] second units + ); + + /// Returns the quotient of two units strings in the CF plain form, or an + /// empty string if either is unknown. Aborts if either string is invalid. + static std::string divide(const std::string &UnitsA, ///< [in] numerator + const std::string &UnitsB ///< [in] denominator + ); + + /// Returns a units string raised to an integer power in the CF plain form, + /// or an empty string if it is unknown. Aborts if the string is invalid. + static std::string power(const std::string &UnitsA, ///< [in] units + int Exponent ///< [in] power + ); + + //--------------------------------------------------------------------------- + /// True if the units are unknown (the field has no units attribute) + bool isUnknown() const; + + /// True if the units are known and dimensionless + bool isDimensionless() const; + + /// Formats the units in the CF plain form: symbols with positive exponents + /// first, then those with negative exponents, each in order of first + /// appearance; "1" if dimensionless; an empty string if unknown + std::string str() const; + + //--------------------------------------------------------------------------- + /// Product of two units expressions; unknown if either is unknown + CFUnits operator*(const CFUnits &Other) const; + + /// Quotient of two units expressions; unknown if either is unknown + CFUnits operator/(const CFUnits &Other) const; + + /// Units expression raised to an integer power; unknown if unknown + CFUnits pow(int Exponent) const; + + /// True if both are unknown, or both are known with the same factors in + /// any order + bool operator==(const CFUnits &Other) const; + bool operator!=(const CFUnits &Other) const; + + private: + /// False for the unknown units of a field with no units attribute + bool Known; + + /// Unit symbols and their non-zero exponents, in order of first appearance + std::vector> Factors; + + /// Parses a units string, aborting on failure, for the string interfaces + static CFUnits parseOrAbort(const std::string &UnitsStr); + + /// Multiplies in one symbol raised to a power, merging with an existing + /// factor of the same symbol and dropping the factor if its exponent + /// becomes zero + void multiplyFactor(const std::string &Symbol, int Exponent); + + /// Parses one whitespace-delimited token (a symbol with an optional + /// exponent, or the digit 1) and multiplies it in with the given sign on + /// the exponent. Returns false if the token is malformed. + bool parseFactor(const std::string &Token, int Sign); + +}; // end class CFUnits + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // OMEGA_CFUNITS_H diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 944df2b42102..1c11cd9f8e7d 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -444,6 +444,17 @@ add_omega_test( "-n;1" ) +########################## +# CFUnits units algebra test +########################## + +add_omega_test( + CFUNITS_TEST + testCFUnits.exe + infra/CFUnitsTest.cpp + "-n;1" +) + add_omega_test( TIMEINTERVAL_PARSE_EXTENDED_FORMATS_TEST testTimeIntervalParseExtendedFormats.exe diff --git a/components/omega/test/infra/CFUnitsTest.cpp b/components/omega/test/infra/CFUnitsTest.cpp new file mode 100644 index 000000000000..3c9092679e0f --- /dev/null +++ b/components/omega/test/infra/CFUnitsTest.cpp @@ -0,0 +1,204 @@ +//===-- Test driver for OMEGA CFUnits ---------------------------*- C++ -*-===// +// +/// \file +/// \brief Test driver for the CFUnits units algebra +/// +/// Checks that units strings in the CF plain form and the udunits caret and +/// slash spellings parse and format as expected, that products, quotients +/// and powers give the plain form with merged exponents, that unknown +/// (empty) units propagate, and that malformed strings are rejected. +// +//===----------------------------------------------------------------------===// + +#include "CFUnits.h" +#include "Error.h" +#include "Logging.h" +#include "MachEnv.h" +#include "mpi.h" + +#include + +using namespace OMEGA; + +namespace { + +int NumFailed = 0; + +// Checks that a units string parses and formats to the expected plain form +void checkParse(const std::string &Input, const std::string &Expected) { + CFUnits Units; + Error Err = CFUnits::parse(Input, Units); + if (!Err.isSuccess()) { + LOG_ERROR("CFUnitsTest: '{}' failed to parse: {}", Input, Err.Msg); + ++NumFailed; + } else if (Units.str() != Expected) { + LOG_ERROR("CFUnitsTest: '{}' formatted as '{}', expected '{}'", Input, + Units.str(), Expected); + ++NumFailed; + } +} + +// Checks that a malformed units string is rejected +void checkInvalid(const std::string &Input) { + CFUnits Units; + Error Err = CFUnits::parse(Input, Units); + if (!Err.isFail()) { + LOG_ERROR("CFUnitsTest: '{}' parsed as '{}' but should be rejected", + Input, Units.str()); + ++NumFailed; + } else if (!Units.isUnknown()) { + LOG_ERROR("CFUnitsTest: '{}' was rejected but left units '{}'", Input, + Units.str()); + ++NumFailed; + } +} + +// Checks a string result against its expected value +void checkResult(const std::string &Label, const std::string &Result, + const std::string &Expected) { + if (Result != Expected) { + LOG_ERROR("CFUnitsTest: {} gave '{}', expected '{}'", Label, Result, + Expected); + ++NumFailed; + } +} + +// Checks a boolean condition +void checkTrue(const std::string &Label, bool Condition) { + if (!Condition) { + LOG_ERROR("CFUnitsTest: {} is false", Label); + ++NumFailed; + } +} + +} // namespace + +int main(int argc, char **argv) { + + MPI_Init(&argc, &argv); + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + initLogging(DefEnv); + + LOG_INFO("----- CFUnitsTest -----"); + + // Parsing and formatting in the CF plain form + checkParse("", ""); + checkParse(" ", ""); + checkParse("1", "1"); + checkParse("m", "m"); + checkParse(" m ", "m"); + checkParse("m s-1", "m s-1"); + checkParse("kg m-3", "kg m-3"); + checkParse("m3 kg-1 K-1", "m3 kg-1 K-1"); + checkParse("degree_C", "degree_C"); + checkParse("m s-1 m", "m2 s-1"); + checkParse("m m-1", "1"); + checkParse("m1", "m"); + checkParse("m0 s", "s"); + + // Positive exponents are written before negative ones + checkParse("s-1 m", "m s-1"); + checkParse("m-3 kg", "kg m-3"); + + // The udunits caret and slash spellings are accepted + checkParse("m^2", "m2"); + checkParse("m^-1 s^-1", "m-1 s-1"); + checkParse("m/s", "m s-1"); + checkParse("m/s/s", "m s-2"); + checkParse("1/s", "s-1"); + checkParse("kg m-3/s", "kg m-3 s-1"); + checkParse("m^2/s", "m2 s-1"); + + // Malformed strings are rejected and leave unknown units + checkInvalid("N m^{-2}"); + checkInvalid("m*s"); + checkInvalid("m.s"); + checkInvalid("m^"); + checkInvalid("m^+2"); + checkInvalid("m--1"); + checkInvalid("m2s"); + checkInvalid("2 m"); + checkInvalid("1e6 m3 s-1"); + checkInvalid("m/"); + checkInvalid("/s"); + checkInvalid("m//s"); + checkInvalid("(m s)"); + checkInvalid("m,s"); + + // Queries + { + CFUnits Unknown, Dimless, Meters; + CFUnits::parse("", Unknown); + CFUnits::parse("1", Dimless); + CFUnits::parse("m", Meters); + checkTrue("empty is unknown", Unknown.isUnknown()); + checkTrue("empty is not dimensionless", !Unknown.isDimensionless()); + checkTrue("1 is dimensionless", Dimless.isDimensionless()); + checkTrue("1 is not unknown", !Dimless.isUnknown()); + checkTrue("m is not unknown", !Meters.isUnknown()); + checkTrue("m is not dimensionless", !Meters.isDimensionless()); + checkTrue("default is unknown", CFUnits().isUnknown()); + } + + // Equality ignores the order of factors and distinguishes unknown + { + CFUnits A, B, C, Unknown; + CFUnits::parse("m s-1", A); + CFUnits::parse("s-1 m", B); + CFUnits::parse("m s-2", C); + checkTrue("m s-1 == s-1 m", A == B); + checkTrue("m s-1 != m s-2", A != C); + checkTrue("m s-1 != unknown", A != Unknown); + checkTrue("unknown == unknown", Unknown == CFUnits()); + } + + // Products, quotients and powers of strings + checkResult("m s-1 * m2", CFUnits::multiply("m s-1", "m2"), "m3 s-1"); + checkResult("kg m-3 * m3", CFUnits::multiply("kg m-3", "m3"), "kg"); + checkResult("m s-1 * kg m-3", CFUnits::multiply("m s-1", "kg m-3"), + "kg m-2 s-1"); + checkResult("m * m-1", CFUnits::multiply("m", "m-1"), "1"); + checkResult("1 * m", CFUnits::multiply("1", "m"), "m"); + checkResult("degree_C * m2", CFUnits::multiply("degree_C", "m2"), + "degree_C m2"); + checkResult("W m-2 * m2", CFUnits::multiply("W m-2", "m2"), "W"); + checkResult("m3 s-1 / m2", CFUnits::divide("m3 s-1", "m2"), "m s-1"); + checkResult("m / m", CFUnits::divide("m", "m"), "1"); + checkResult("1 / s", CFUnits::divide("1", "s"), "s-1"); + checkResult("(m s-1)^2", CFUnits::power("m s-1", 2), "m2 s-2"); + checkResult("(m s-1)^-1", CFUnits::power("m s-1", -1), "s m-1"); + checkResult("m^0", CFUnits::power("m", 0), "1"); + + // Unknown units propagate through every operation + checkResult("unknown * m2", CFUnits::multiply("", "m2"), ""); + checkResult("m2 * unknown", CFUnits::multiply("m2", ""), ""); + checkResult("unknown / m2", CFUnits::divide("", "m2"), ""); + checkResult("m2 / unknown", CFUnits::divide("m2", ""), ""); + checkResult("unknown^3", CFUnits::power("", 3), ""); + checkResult("unknown^0", CFUnits::power("", 0), ""); + + // The same algebra on parsed values + { + CFUnits Vel, Area, Dens; + CFUnits::parse("m s-1", Vel); + CFUnits::parse("m2", Area); + CFUnits::parse("kg m-3", Dens); + checkResult("value product", (Vel * Area).str(), "m3 s-1"); + checkResult("value quotient", (Vel / Area).str(), "m-1 s-1"); + checkResult("value chain", (Vel * Area * Dens).str(), "kg s-1"); + checkResult("value power", Dens.pow(-1).str(), "m3 kg-1"); + checkResult("value unknown", (Vel * CFUnits()).str(), ""); + } + + if (NumFailed > 0) { + ABORT_ERROR("CFUnitsTest: {} checks failed", NumFailed); + } + + LOG_INFO("----- CFUnitsTest Successful -----"); + + MPI_Barrier(MPI_COMM_WORLD); + MPI_Finalize(); + return 0; +} From 37d2f5e108590ad4a03147f9a2793407398617e3 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 13 Sep 2026 02:20:29 -0500 Subject: [PATCH 2/5] Derive analysis output metadata from the input field The analysis operators registered their output fields with empty units and standard name, so every global statistic Omega wrote arrived with no units (#545). A mean, minimum, maximum, standard deviation or time mean has the units and standard name of the field it reduces, and CF describes the reduction itself with a cell_methods attribute rather than a different standard name. AnalysisOperator gains helpers for this: inheritMetadata() collects an input's units, standard_name and cell_methods, each empty if the input has none, and appends the cell method of the new reduction so the attribute lists reductions in the order applied; spatialCellMethod() gives "area: " for a horizontal field and "area: depth: " for a layered one; createOutputField() creates the output Field with them. The five operators use these, so a time mean of a spatial mean now carries "area: depth: mean time: mean". An operator whose result has different units derives them with CFUnits before creating its output. Only the CF attribute names are read and written, and a missing attribute is treated as empty, so this works the same whether or not Field::create stores empty attributes. The operator CTest checks the attributes every operator writes, a chained time mean of a spatial mean, and that a field without units yields outputs without units. The Analysis developer's and user's guides describe the output metadata and how a new operator derives it. Fixes #545 Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/doc/devGuide/Analysis.md | 59 +++++++- components/omega/doc/userGuide/Analysis.md | 6 +- .../omega/src/analysis/AnalysisOperator.cpp | 80 ++++++++++ .../omega/src/analysis/AnalysisOperator.h | 60 ++++++++ .../src/analysis/operators/SpatialMaxOp.h | 25 ++-- .../src/analysis/operators/SpatialMeanOp.h | 21 +-- .../src/analysis/operators/SpatialMinOp.h | 25 ++-- .../src/analysis/operators/SpatialStdDevOp.h | 26 ++-- .../omega/src/analysis/operators/TimeMeanOp.h | 30 ++-- .../test/analysis/AnalysisOperatorTest.cpp | 137 ++++++++++++++++-- 10 files changed, 397 insertions(+), 72 deletions(-) diff --git a/components/omega/doc/devGuide/Analysis.md b/components/omega/doc/devGuide/Analysis.md index 6199fb3b7bfb..0b4aec0d3b15 100644 --- a/components/omega/doc/devGuide/Analysis.md +++ b/components/omega/doc/devGuide/Analysis.md @@ -144,7 +144,9 @@ A concrete derived operator is responsible for three things: record the input Field names, construct the output Field name(s), allocate the output data array, and register the output [Field](#omega-dev-field) with the appropriate metadata. Output Field names follow the convention - `_` (for example, `Temperature_SpatialMax`). + `_` (for example, `Temperature_SpatialMax`). The + output's CF metadata is derived from the input's, as described under + [Output metadata](#omega-dev-analysis-output-metadata). 2. **Initialization**: `initialize()` is called by the orchestrator after all Fields exist. The base implementation stores the mesh, vertical coordinate, and communicator; operators may override to perform additional @@ -245,6 +247,53 @@ normalization. The sums accumulate each owned entity's terms in a fixed order and reduce them with `globalSum`, so they are as reproducible as Omega's other global sums. +(omega-dev-analysis-output-metadata)= +### Output metadata + +An output Field inherits its CF metadata from the input it is derived from. +The base class provides the helpers that do this: + +- `inheritMetadata(InputName, CellMethod)` returns an `InheritedMetadata` + struct holding the input's `units`, `standard_name` and `cell_methods`, + each empty if the input has none, with `CellMethod` appended to the cell + methods. CF requires `cell_methods` to list the reductions in the order + they were applied, so a time mean of a spatial mean carries + `area: depth: mean time: mean`. +- `spatialCellMethod(InputName, Method)` returns the cell method for a + reduction over all owned mesh entities and layers of the input: + `area: ` for a horizontal (1D) field and `area: depth: ` + for a field with a vertical dimension. `Method` is a CF cell method name + (`mean`, `minimum`, `maximum`, `standard_deviation`, `sum`). +- `createOutputField(...)` calls `Field::create()` with the inherited units + and standard name and adds the `cell_methods` attribute when it is not + empty. +- `getFieldAttribute(FieldName, AttName)` returns any string attribute of a + Field, or an empty string if the Field has none. + +The built-in reductions (mean, minimum, maximum, standard deviation, time +mean) all have the units and standard name of the field they reduce, so +they pass `inheritMetadata()` straight to `createOutputField()`: +```c++ +auto CellMethod = spatialCellMethod(InputNames[0], "mean"); +auto Meta = inheritMetadata(InputNames[0], CellMethod); +auto OutputField = createOutputField(OutputNames[0], + "Spatial mean of " + InputNames[0], + Meta, ValidMin, ValidMax, NDims, DimNames); +``` +An operator whose result has different units derives them with +[CFUnits](#omega-dev-cfunits) before creating the output, and clears the +standard name if the result is no longer the same physical quantity. For +example, an area-weighted sum multiplies the input's units by those of the +cell area: +```c++ +auto Meta = inheritMetadata(InputNames[0], + spatialCellMethod(InputNames[0], "sum")); +Meta.Units = CFUnits::multiply(Meta.Units, "m2"); +Meta.StdName.clear(); +``` +A field without units yields outputs without units, never with invented ones, +and a malformed units string aborts the run when the operator is constructed. + ## Operator factory and type dispatch Fields carry runtime metadata describing their scalar type, rank, and memory @@ -493,9 +542,11 @@ To add a new operator: operator type name to the base constructor, set `InputNames`, build the output name as `_`, set `OutputNames` and `InstanceName`, allocate the output data array, and register the output Field with - `Field::create()` plus `attachData()`. If the operator depends on another - operator's output (as `SpatialStdDev` depends on `SpatialMean`), append - that dependency to `InputNames`. + `createOutputField()` plus `attachData()`, deriving its units, standard + name and cell methods from the input as described under + [Output metadata](#omega-dev-analysis-output-metadata). If the operator + depends on another operator's output (as `SpatialStdDev` depends on + `SpatialMean`), append that dependency to `InputNames`. 2. **Implement `compute()`** to retrieve the input array(s) from the Field registry, perform the transformation, write the output array(s), and update `LastComputed`/`FieldComputed`. diff --git a/components/omega/doc/userGuide/Analysis.md b/components/omega/doc/userGuide/Analysis.md index e760f15b3e5f..d08540aeee5b 100644 --- a/components/omega/doc/userGuide/Analysis.md +++ b/components/omega/doc/userGuide/Analysis.md @@ -101,7 +101,11 @@ user-specified fields. The mean and standard deviation are weighted so that they do not depend on the mesh resolution: by the area of each cell, edge or vertex for a horizontal field, and by mass (area times pseudo-thickness) for a field with a vertical dimension, with an interface field weighted by half -the mass of each adjacent layer. Only active layers take part. +the mass of each adjacent layer. Only active layers take part. Each statistic +is written with the `units` and `standard_name` of the field it reduces and a +CF `cell_methods` attribute describing the reduction, for example +`area: depth: mean` for the mean of a layered field and +`area: depth: mean time: mean` for its time mean. **Example:** diff --git a/components/omega/src/analysis/AnalysisOperator.cpp b/components/omega/src/analysis/AnalysisOperator.cpp index 59a0a44d4b29..3425dc592530 100644 --- a/components/omega/src/analysis/AnalysisOperator.cpp +++ b/components/omega/src/analysis/AnalysisOperator.cpp @@ -87,6 +87,86 @@ const std::vector AnalysisOperator::getOutputFieldNames() { return OutputNames; } // end getOutputFieldNames +//------------------------------------------------------------------------------ +// Returns a string attribute of a Field, or an empty string if the Field has +// no such attribute. Analysis output is derived from fields that may lack +// units or a standard name, so absence is the empty string, not an error. +std::string AnalysisOperator::getFieldAttribute(const std::string &FieldName, + const std::string &AttName) { + + auto FieldPtr = Field::get(FieldName); + OMEGA_REQUIRE(FieldPtr, "AnalysisOperator: field {} not found", FieldName); + + std::string Value; + if (FieldPtr->hasMetadata(AttName)) { + Error Err = FieldPtr->getMetadata(AttName, Value); + CHECK_ERROR_ABORT(Err, + "AnalysisOperator: attribute {} of field {} is " + "not a string", + AttName, FieldName); + } + return Value; +} // end getFieldAttribute + +//------------------------------------------------------------------------------ +// Collects the units, standard name and cell methods an output inherits from +// its input. A new cell method is appended to the input's cell_methods so the +// attribute lists the reductions in the order they were applied, as CF +// requires (e.g. a time mean of a spatial mean is "area: mean time: mean"). +AnalysisOperator::InheritedMetadata +AnalysisOperator::inheritMetadata(const std::string &InputName, + const std::string &CellMethod) { + + InheritedMetadata Meta; + Meta.Units = getFieldAttribute(InputName, "units"); + Meta.StdName = getFieldAttribute(InputName, "standard_name"); + Meta.CellMethods = getFieldAttribute(InputName, "cell_methods"); + + if (!CellMethod.empty()) { + if (!Meta.CellMethods.empty()) + Meta.CellMethods += " "; + Meta.CellMethods += CellMethod; + } + return Meta; +} // end inheritMetadata + +//------------------------------------------------------------------------------ +// Returns the cell method for a reduction over all owned mesh entities and +// layers of an input Field. The spatial operators treat a 1D field as +// horizontal only and the last dimension of any higher-rank field as +// vertical, so the method covers "area" alone or "area" and "depth". Both +// are names CF allows in cell_methods without a matching coordinate. +std::string AnalysisOperator::spatialCellMethod(const std::string &InputName, + const std::string &Method) { + + auto FieldPtr = Field::get(InputName); + OMEGA_REQUIRE(FieldPtr, "AnalysisOperator: field {} not found", InputName); + + if (FieldPtr->getNumDims() >= 2) + return "area: depth: " + Method; + return "area: " + Method; +} // end spatialCellMethod + +//------------------------------------------------------------------------------ +// Creates an output Field carrying the inherited CF metadata. The cell +// methods are stored as a std::string so the metadata can be retrieved and +// written as a string attribute. +std::shared_ptr AnalysisOperator::createOutputField( + const std::string &OutputName, const std::string &Description, + const InheritedMetadata &Meta, const std::any ValidMin, + const std::any ValidMax, const int NumDims, + const std::vector &DimNames) { + + auto OutputField = + Field::create(OutputName, Description, Meta.Units, Meta.StdName, + ValidMin, ValidMax, NumDims, DimNames); + + if (!Meta.CellMethods.empty()) + OutputField->addMetadata("cell_methods", Meta.CellMethods); + + return OutputField; +} // end createOutputField + //------------------------------------------------------------------------------ // Checks whether the operator's output is valid for the given timestamp. // Returns true if the operator has been computed and the LastComputed diff --git a/components/omega/src/analysis/AnalysisOperator.h b/components/omega/src/analysis/AnalysisOperator.h index aa3800c7090a..ee708255e3bb 100644 --- a/components/omega/src/analysis/AnalysisOperator.h +++ b/components/omega/src/analysis/AnalysisOperator.h @@ -20,6 +20,13 @@ /// data from the Field registry and writes results to operator-owned output /// arrays. /// +/// An output Field inherits its CF metadata from the input it is derived +/// from: the base class provides helpers that collect an input's units, +/// standard name and cell methods, append the cell method describing the new +/// reduction, and create the output Field with them. Operators whose result +/// has different units (an area-weighted sum, a product) derive them with +/// CFUnits before creating the output Field. +/// /// This file also provides helper functions (opParam, makeOpConfig) for /// constructing Config objects inline when instantiating operators, /// providing a uniform parameter-passing mechanism whether operators are @@ -27,6 +34,7 @@ /// //===----------------------------------------------------------------------===// +#include "CFUnits.h" #include "Config.h" #include "DataTypes.h" #include "Dimension.h" @@ -165,6 +173,58 @@ class AnalysisOperator { ); protected: + /// CF metadata an output Field inherits from the input it is derived + /// from. Units and standard name are copied and are empty when the input + /// has none. CellMethods records every reduction applied so far, in + /// order, as the CF cell_methods attribute (e.g. "area: mean time: mean"). + struct InheritedMetadata { + std::string Units; ///< units of the input, empty if unknown + std::string StdName; ///< CF standard name, empty if none + std::string CellMethods; ///< CF cell_methods, empty if none + }; + + /// Collects the CF metadata an output Field inherits from an input Field, + /// with the cell method of the new reduction (if any) appended to those + /// the input already records. The units and standard name are copied + /// unchanged; an operator that changes the units or invalidates the + /// standard name edits the result before creating its output Field. + static InheritedMetadata + inheritMetadata(const std::string &InputName, ///< [in] input field name + const std::string &CellMethod = "" ///< [in] new reduction + ); + + /// Returns the CF cell method for a spatial reduction of an input Field + /// over all of its owned mesh entities and layers: "area: " for a + /// horizontal field and "area: depth: " for a field with a + /// vertical dimension. Method is a CF cell method name (mean, minimum, + /// maximum, standard_deviation, sum). + static std::string + spatialCellMethod(const std::string &InputName, ///< [in] input field name + const std::string &Method ///< [in] CF cell method + ); + + /// Returns the value of a string attribute (units, standard_name, + /// cell_methods, long_name) of a Field, or an empty string if the Field + /// has no such attribute + static std::string + getFieldAttribute(const std::string &FieldName, ///< [in] field name + const std::string &AttName ///< [in] attribute name + ); + + /// Creates an output Field with the given description, valid range and + /// dimensions, carrying the inherited units, standard name and cell + /// methods. Empty units or standard name are passed through to + /// Field::create; an empty cell_methods is not added. + static std::shared_ptr createOutputField( + const std::string &OutputName, ///< [in] output field name + const std::string &Description, ///< [in] long name of the output + const InheritedMetadata &Meta, ///< [in] inherited CF metadata + const std::any ValidMin, ///< [in] min valid value + const std::any ValidMax, ///< [in] max valid value + const int NumDims, ///< [in] number of dimensions + const std::vector &DimNames ///< [in] dimension names + ); + const HorzMesh *Mesh; ///< Horizontal mesh for spatial operations const VertCoord *VCoord; ///< Vertical coordinate for vertical ops MPI_Comm Comm; ///< MPI communicator for parallel reductions diff --git a/components/omega/src/analysis/operators/SpatialMaxOp.h b/components/omega/src/analysis/operators/SpatialMaxOp.h index ab22b1ba5140..b40bd9bea66b 100644 --- a/components/omega/src/analysis/operators/SpatialMaxOp.h +++ b/components/omega/src/analysis/operators/SpatialMaxOp.h @@ -65,17 +65,20 @@ template class SpatialMaxOp : public AnalysisOperator { DimNames[0] = "Scalar"; auto ScalarDim = Dimension::create(DimNames[0], 1); - // Register output Field with metadata - auto OutputField = - Field::create(OutputNames[0], - "Spatial maximum of " + InputNames[0], // Description - "", // Units - "", // Standard name - -std::numeric_limits::max(), // Min valid value - std::numeric_limits::max(), // Max valid value - NDims, // Rank - DimNames // Dimension names - ); + // Register output Field with metadata. A maximum has the units and + // standard name of the field it reduces; its cell_methods record the + // reduction after any the input already carries. + auto CellMethod = spatialCellMethod(InputNames[0], "maximum"); + auto Meta = inheritMetadata(InputNames[0], CellMethod); + auto OutputField = createOutputField( + OutputNames[0], + "Spatial maximum of " + InputNames[0], // Description + Meta, // CF metadata + -std::numeric_limits::max(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Rank + DimNames // Dimension names + ); // Attach output data array to Field OutputField->template attachData>(OutputData); diff --git a/components/omega/src/analysis/operators/SpatialMeanOp.h b/components/omega/src/analysis/operators/SpatialMeanOp.h index ebff0862b86a..fc6b39fe5bc1 100644 --- a/components/omega/src/analysis/operators/SpatialMeanOp.h +++ b/components/omega/src/analysis/operators/SpatialMeanOp.h @@ -66,16 +66,19 @@ template class SpatialMeanOp : public AnalysisOperator { DimNames[0] = "Scalar"; auto ScalarDim = Dimension::create(DimNames[0], 1); - // Register output Field with metadata + // Register output Field with metadata. A mean has the units and + // standard name of the field it reduces; its cell_methods record the + // reduction after any the input already carries. + auto Meta = inheritMetadata(InputNames[0], + spatialCellMethod(InputNames[0], "mean")); auto OutputField = - Field::create(OutputNames[0], - "Spatial mean of " + InputNames[0], // Description - "", // Units - "", // Standard name - -std::numeric_limits::max(), // Min valid value - std::numeric_limits::max(), // Max valid value - NDims, // Rank - DimNames // Dimension names + createOutputField(OutputNames[0], + "Spatial mean of " + InputNames[0], // Description + Meta, // CF metadata + -std::numeric_limits::max(), // Min valid + std::numeric_limits::max(), // Max valid + NDims, // Rank + DimNames // Dim names ); // Attach output data array to Field diff --git a/components/omega/src/analysis/operators/SpatialMinOp.h b/components/omega/src/analysis/operators/SpatialMinOp.h index dcd2c98f3d59..c0b57051119e 100644 --- a/components/omega/src/analysis/operators/SpatialMinOp.h +++ b/components/omega/src/analysis/operators/SpatialMinOp.h @@ -65,17 +65,20 @@ template class SpatialMinOp : public AnalysisOperator { DimNames[0] = "Scalar"; auto ScalarDim = Dimension::create(DimNames[0], 1); - // Register output Field with metadata - auto OutputField = - Field::create(OutputNames[0], - "Spatial minimum of " + InputNames[0], // Description - "", // Units - "", // Standard name - -std::numeric_limits::max(), // Min valid value - std::numeric_limits::max(), // Max valid value - NDims, // Rank - DimNames // Dimension names - ); + // Register output Field with metadata. A minimum has the units and + // standard name of the field it reduces; its cell_methods record the + // reduction after any the input already carries. + auto CellMethod = spatialCellMethod(InputNames[0], "minimum"); + auto Meta = inheritMetadata(InputNames[0], CellMethod); + auto OutputField = createOutputField( + OutputNames[0], + "Spatial minimum of " + InputNames[0], // Description + Meta, // CF metadata + -std::numeric_limits::max(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Rank + DimNames // Dimension names + ); // Attach output data array to Field OutputField->template attachData>(OutputData); diff --git a/components/omega/src/analysis/operators/SpatialStdDevOp.h b/components/omega/src/analysis/operators/SpatialStdDevOp.h index dff76cbb712f..492048f5afe4 100644 --- a/components/omega/src/analysis/operators/SpatialStdDevOp.h +++ b/components/omega/src/analysis/operators/SpatialStdDevOp.h @@ -67,17 +67,21 @@ template class SpatialStdDevOp : public AnalysisOperator { DimNames[0] = "Scalar"; auto ScalarDim = Dimension::create(DimNames[0], 1); - // Register output Field with metadata - auto OutputField = - Field::create(OutputNames[0], - "Standard deviation of " + InputNames[0], // Description - "", // Units - "", // Standard name - static_cast(0), // Min valid value - std::numeric_limits::max(), // Max valid value - NDims, // Rank - DimNames // Dimension names - ); + // Register output Field with metadata. A standard deviation has the + // units and standard name of the field it reduces; its cell_methods + // record the reduction after any the input already carries. + auto Meta = inheritMetadata( + InputNames[0], + spatialCellMethod(InputNames[0], "standard_deviation")); + auto OutputField = createOutputField( + OutputNames[0], + "Standard deviation of " + InputNames[0], // Description + Meta, // CF metadata + static_cast(0), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Rank + DimNames // Dimension names + ); // Attach output data array to Field OutputField->template attachData(OutputData); diff --git a/components/omega/src/analysis/operators/TimeMeanOp.h b/components/omega/src/analysis/operators/TimeMeanOp.h index bddb93a5534d..075709c57153 100644 --- a/components/omega/src/analysis/operators/TimeMeanOp.h +++ b/components/omega/src/analysis/operators/TimeMeanOp.h @@ -83,30 +83,28 @@ template class TimeMeanOp : public AnalysisOperator { std::vector DimNames; InputField->getDimNames(DimNames); - // Fetch input metadata - std::string InputDescr, InputUnits, InputStdName; + // Fetch the input description and valid range + std::string InputDescr; ScalarT InputValidMin, InputValidMax; InputField->getMetadata("Description", InputDescr); - InputField->getMetadata("Units", InputUnits); - InputField->getMetadata("StdName", InputStdName); InputField->getMetadata("ValidMin", InputValidMin); InputField->getMetadata("ValidMax", InputValidMax); - // Register output Field with same dimensions as input but Real type + // Register output Field with same dimensions as input but Real type. + // A time mean has the units and standard name of the field it + // averages; its cell_methods record the reduction after any the input + // already carries (e.g. "area: mean time: mean"). + auto Meta = inheritMetadata(InputNames[0], "time: mean"); auto OutputField = - Field::create(OutputNames[0], - "Time average of " + InputDescr, // Description - InputUnits, // Units - InputStdName, // Standard name - static_cast(InputValidMin), // Min valid value - static_cast(InputValidMax), // Max valid value - NDims, // Rank - DimNames // Dimension names + createOutputField(OutputNames[0], + "Time average of " + InputDescr, // Description + Meta, // CF metadata + static_cast(InputValidMin), // Min valid + static_cast(InputValidMax), // Max valid + NDims, // Rank + DimNames // Dim names ); - // Stamp CF-compliant cell_methods for time-mean output. - OutputField->addMetadata("cell_methods", std::string("time: mean")); - // Store array size for parallel iteration ArraySize = static_cast(InputData.size()); diff --git a/components/omega/test/analysis/AnalysisOperatorTest.cpp b/components/omega/test/analysis/AnalysisOperatorTest.cpp index f324147d61d0..f9593219fb84 100644 --- a/components/omega/test/analysis/AnalysisOperatorTest.cpp +++ b/components/omega/test/analysis/AnalysisOperatorTest.cpp @@ -34,6 +34,11 @@ int NumTests = 0; int NumPassed = 0; int NumFailed = 0; +// Units and standard name given to every test field, which the operators +// must pass through to their outputs +const std::string TestUnits = "m"; +const std::string TestStdName = "test_field_standard_name"; + //===----------------------------------------------------------------------===// // Generic Helper Template Struct //===----------------------------------------------------------------------===// @@ -112,9 +117,9 @@ template struct TestHelper { ValidMax = 1.0e30; } - auto TestField = - Field::create(FieldName, "Test field for multi-type validation", "m", - "", ValidMin, ValidMax, 1, DimNames); + auto TestField = Field::create( + FieldName, "Test field for multi-type validation", TestUnits, + TestStdName, ValidMin, ValidMax, 1, DimNames); ArrayType TestData(FieldName + "_data", Dims[0]); TestField->template attachData(TestData, false); @@ -147,9 +152,9 @@ template struct TestHelper { ValidMax = 1.0e30; } - auto TestField = - Field::create(FieldName, "Test field for multi-type validation", "m", - "", ValidMin, ValidMax, 2, DimNames); + auto TestField = Field::create( + FieldName, "Test field for multi-type validation", TestUnits, + TestStdName, ValidMin, ValidMax, 2, DimNames); ArrayType TestData(FieldName + "_data", Dims[0], Dims[1]); TestField->template attachData(TestData, false); @@ -184,9 +189,9 @@ template struct TestHelper { ValidMax = 1.0e30; } - auto TestField = - Field::create(FieldName, "Test field for multi-type validation", "m", - "", ValidMin, ValidMax, 3, DimNames); + auto TestField = Field::create( + FieldName, "Test field for multi-type validation", TestUnits, + TestStdName, ValidMin, ValidMax, 3, DimNames); ArrayType TestData(FieldName + "_data", Dims[0], Dims[1], Dims[2]); TestField->template attachData(TestData, false); @@ -324,6 +329,54 @@ bool checkClose(const std::string &TestName, Real Computed, Real Expected, return Passed; } +//------------------------------------------------------------------------------ +// Returns a string attribute of a field, or an empty string if the field has +// no such attribute. Field::create may store an empty attribute or omit it; +// both mean the field has none. +std::string getAttribute(const std::string &FieldName, + const std::string &AttName) { + auto FieldPtr = Field::get(FieldName); + std::string Value; + if (FieldPtr && FieldPtr->hasMetadata(AttName)) + FieldPtr->getMetadata(AttName, Value); + return Value; +} + +//------------------------------------------------------------------------------ +// Checks that an operator's output field carries the expected units, standard +// name and cell methods +void checkOutputMetadata(const std::string &TestName, + const std::string &OutputName, + const std::string &ExpectedUnits, + const std::string &ExpectedStdName, + const std::string &ExpectedCellMethods) { + + std::string Units = getAttribute(OutputName, "units"); + std::string StdName = getAttribute(OutputName, "standard_name"); + std::string CellMethods = getAttribute(OutputName, "cell_methods"); + + bool Passed = (Units == ExpectedUnits) && (StdName == ExpectedStdName) && + (CellMethods == ExpectedCellMethods); + reportTest(TestName + " metadata", Passed); + + if (!Passed) { + LOG_ERROR(" units: '{}', expected '{}'", Units, ExpectedUnits); + LOG_ERROR(" standard_name: '{}', expected '{}'", StdName, + ExpectedStdName); + LOG_ERROR(" cell_methods: '{}', expected '{}'", CellMethods, + ExpectedCellMethods); + } +} + +//------------------------------------------------------------------------------ +// Returns the cell method the spatial operators record for a test field of +// the given rank: 1D fields are horizontal only, higher ranks have depth +std::string expectedSpatialCellMethod(int Rank, const std::string &Method) { + if (Rank >= 2) + return "area: depth: " + Method; + return "area: " + Method; +} + //===----------------------------------------------------------------------===// // Operator Test Templates //===----------------------------------------------------------------------===// @@ -416,6 +469,9 @@ void testSpatialMaxOpType(const std::string &TypeName, const MachEnv *Env, bool Passed = (std::abs(ComputedMax - ExpectedMaxReal) <= static_cast(Helper::getTolerance())); reportTest("SpatialMaxOp: " + TypeName, Passed); + checkOutputMetadata("SpatialMaxOp: " + TypeName, FieldName + "_SpatialMax", + TestUnits, TestStdName, + expectedSpatialCellMethod(Rank, "maximum")); if (!Passed) { LOG_ERROR(" Expected: {}, Got: {}", ExpectedMaxReal, ComputedMax); @@ -492,6 +548,9 @@ void testSpatialMinOpType(const std::string &TypeName, const MachEnv *Env, bool Passed = (std::abs(ComputedMin - ExpectedMinReal) <= static_cast(Helper::getTolerance())); reportTest("SpatialMinOp: " + TypeName, Passed); + checkOutputMetadata("SpatialMinOp: " + TypeName, FieldName + "_SpatialMin", + TestUnits, TestStdName, + expectedSpatialCellMethod(Rank, "minimum")); if (!Passed) { LOG_ERROR(" Expected: {}, Got: {}", ExpectedMinReal, ComputedMin); @@ -576,6 +635,9 @@ void testSpatialMeanOpType(const std::string &TypeName, const MachEnv *Env, // Verify checkClose("SpatialMeanOp: " + TypeName, ComputedMean, ExpectedMean, Helper::getRelTolerance()); + checkOutputMetadata("SpatialMeanOp: " + TypeName, FieldName + "_SpatialMean", + TestUnits, TestStdName, + expectedSpatialCellMethod(Rank, "mean")); } //------------------------------------------------------------------------------ @@ -662,6 +724,9 @@ void testSpatialStdDevOpType(const std::string &TypeName, const MachEnv *Env, // Verify checkClose("SpatialStdDevOp: " + TypeName, ComputedStdDev, ExpectedStdDev, Helper::getRelTolerance()); + checkOutputMetadata("SpatialStdDevOp: " + TypeName, + FieldName + "_SpatialStdDev", TestUnits, TestStdName, + expectedSpatialCellMethod(Rank, "standard_deviation")); } //------------------------------------------------------------------------------ @@ -861,6 +926,9 @@ void testTimeMeanOpType(const std::string &TypeName, const MachEnv *Env, } reportTest("TimeMeanOp: " + TypeName, Passed); + checkOutputMetadata("TimeMeanOp: " + TypeName, + FieldName + "_TimeMean" + PeriodLabel, TestUnits, + TestStdName, "time: mean"); if (!Passed) { LOG_ERROR(" Expected mean: {}", ExpectedMean); @@ -1114,6 +1182,55 @@ void testWeightedStats(const MachEnv *Env, const HorzMesh *Mesh, } } +//------------------------------------------------------------------------------ +// Tests that metadata is inherited along an operator chain: a time mean of a +// spatial mean keeps the units and standard name and lists both reductions in +// order in cell_methods. Also tests that a field without units or a standard +// name yields outputs without them rather than with invented ones. +void testInheritedMetadata(const MachEnv *Env, const HorzMesh *Mesh, + const VertCoord *VCoord, Clock *ModelClock) { + + using Helper = TestHelper; + + std::vector Dims = Helper::getDims(Mesh, VCoord); + Config EmptyConfig; + + // Chain: SpatialMean then TimeMean of a 2D field with units "m" + std::string FieldName = "TestFieldChain"; + Helper::createField(FieldName, Dims, + [](I4 i, I4 j) -> R8 { return static_cast(1); }); + + auto MeanOp = + AnalysisOpFactory::createOp("SpatialMean", {FieldName}, EmptyConfig); + MeanOp->initialize(Env, Mesh, VCoord, EmptyConfig); + + std::string MeanName = FieldName + "_SpatialMean"; + std::string Period = "1day"; + auto ChainOp = AnalysisOpFactory::createOp( + "TimeMean", {MeanName}, makeOpConfig(opParam("Period", Period))); + ChainOp->initialize(Env, Mesh, VCoord, EmptyConfig); + + checkOutputMetadata("TimeMean of SpatialMean", + MeanName + "_TimeMean" + Period, TestUnits, TestStdName, + "area: depth: mean time: mean"); + + // A field with no units and no standard name: the outputs have none either + std::string BareName = "TestFieldBare"; + std::vector DimNames = Helper::getDimNames(); + auto BareField = Field::create(BareName, "Test field without units", "", "", + -1.0e30, 1.0e30, 2, DimNames); + Array2DR8 BareData(BareName + "_data", Dims[0], Dims[1]); + BareField->attachData(BareData); + + auto BareMaxOp = + AnalysisOpFactory::createOp("SpatialMax", {BareName}, EmptyConfig); + BareMaxOp->initialize(Env, Mesh, VCoord, EmptyConfig); + + checkOutputMetadata("SpatialMax of field without units", + BareName + "_SpatialMax", "", "", + "area: depth: maximum"); +} + //===----------------------------------------------------------------------===// // Main Test Functions //===----------------------------------------------------------------------===// @@ -3010,6 +3127,8 @@ int main(int argc, char *argv[]) { testTransectAccumulatorOp(DefEnv, Mesh, VCoord); + testInheritedMetadata(DefEnv, Mesh, VCoord, ModelClock); + if (NumFailed > 0) { Err = 1; LOG_ERROR("AnalysisOperatorTest failure"); From fc45da526625e105067f9abd517b77d7107b8802 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 13 Sep 2026 02:55:42 -0500 Subject: [PATCH 3/5] Keep braces in error message arguments from being read as placeholders Error::createErrMsg resumed its scan for {} placeholders at the position of the placeholder it had just replaced, so an argument that itself contained braces (a units string such as m^{-2}) was scanned again and its braces taken for further placeholders, logging "Not enough arguments for placeholders in error msg". Resume the scan after the inserted argument instead. ErrorTest now covers an argument with braces. Co-Authored-By: Claude Opus 5 (1M context) --- components/omega/src/infra/Error.cpp | 3 +++ components/omega/test/infra/ErrorTest.cpp | 18 ++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/components/omega/src/infra/Error.cpp b/components/omega/src/infra/Error.cpp index f364b1d44ed6..e77c4dee9edc 100644 --- a/components/omega/src/infra/Error.cpp +++ b/components/omega/src/infra/Error.cpp @@ -109,6 +109,9 @@ std::string Error::createErrMsg( int Length = EndPos - StartPos + 1; if (ArgIndex < MsgArgs.size()) { Result.replace(StartPos, Length, MsgArgs[ArgIndex]); + // Resume scanning after the inserted argument so that braces + // inside it are not taken for placeholders + EndPos = StartPos + MsgArgs[ArgIndex].size(); ArgIndex++; } else { LOG_ERROR("Not enough arguments for placeholders in error msg"); diff --git a/components/omega/test/infra/ErrorTest.cpp b/components/omega/test/infra/ErrorTest.cpp index 1cc4ea4e9dd3..83718bb28bd2 100644 --- a/components/omega/test/infra/ErrorTest.cpp +++ b/components/omega/test/infra/ErrorTest.cpp @@ -132,6 +132,24 @@ int main(int argc, char **argv) { LOG_ERROR(" Error constructor with message actual: {}", TotalError.Msg); } + // Test constructor with an argument that itself contains braces, which + // must be inserted verbatim and not be taken for a placeholder + { + std::string BraceLine = std::to_string(__LINE__ + 3); + std::string BraceExpect = " [error] [ErrorTest.cpp:" + BraceLine + + "] Bad units m^{-2} in field\n"; + Error BraceError(ErrorCode::Fail, __LINE__, __FILE__, + "Bad units {} in field", "m^{-2}"); + if (BraceError.isFail() and BraceError.Msg == BraceExpect) { + LOG_INFO("Error constructor with braces in argument: PASS"); + } else { + ++RetVal; + LOG_ERROR("Error constructor with braces in argument: FAIL"); + LOG_ERROR(" expected: {}", BraceExpect); + LOG_ERROR(" actual: {}", BraceError.Msg); + } + } + // Test reset of error code TotalError.reset(); if (TotalError.isSuccess()) { From 6c7e67a67ce9bc101eb3a1ab2566846229203d12 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 4 Oct 2026 07:42:18 -0500 Subject: [PATCH 4/5] Add scale factors to CFUnits A field whose values are multiplied by a constant, such as a transport reported in Sverdrups by multiplying by 1e-6, has its units scaled by the reciprocal. CFUnits could not express that, so an operator doing it could only drop the units or keep wrong ones. A units string may now start with a positive number, the scale factor of a scaled unit as udunits writes it (1e+06 m3 s-1). The factor multiplies through products, quotients and powers, takes part in equality to within rounding, and is formatted with the fewest digits that reproduce it, so 1/1e-6 is written 1e+06. CFUnits::scale() and scaled() multiply units by a positive factor. A number anywhere but before the unit symbols, or one that is not positive and finite, is still a parse error. Co-Authored-By: Claude Opus 5.5 (1M context) --- components/omega/doc/devGuide/CFUnits.md | 16 ++- components/omega/src/infra/CFUnits.cpp | 109 +++++++++++++++++++- components/omega/src/infra/CFUnits.h | 38 ++++++- components/omega/test/infra/CFUnitsTest.cpp | 52 +++++++++- 4 files changed, 199 insertions(+), 16 deletions(-) diff --git a/components/omega/doc/devGuide/CFUnits.md b/components/omega/doc/devGuide/CFUnits.md index f550619cf7a5..407ad51de2f3 100644 --- a/components/omega/doc/devGuide/CFUnits.md +++ b/components/omega/doc/devGuide/CFUnits.md @@ -23,6 +23,18 @@ a symbol whose exponent reaches zero is dropped, so `kg m-3` times `m3` is `kg`. In the result, symbols with positive exponents come first, then those with negative exponents, each in the order they first appeared. +A units string may also start with a positive number, the scale factor of a +scaled unit as udunits writes it: `1e+06 m3 s-1` is a million cubic meters +per second. Multiplying a field's values by a constant, for example to report +a transport in units of 10^6 m3 s-1, scales its units by the reciprocal: +```c++ +std::string Transport = CFUnits::scale("m3 s-1", 1.0 / 1.0e-6); // "1e+06 m3 s-1" +``` +Scale factors multiply through products, quotients and powers, and are +written with the fewest digits that reproduce them, so the rounding error in +`1.0 / 1.0e-6` does not show. A number is only accepted before the unit +symbols, and must be positive and finite. + A field with no units attribute has *unknown* units, represented by the empty string. Unknown units propagate through every operation, so a field derived from one without units has no units either, rather than wrong ones. @@ -38,8 +50,8 @@ CFUnits Area; CFUnits::parse("m2", Area); std::string Flux = (Units * Area).str(); // "m3 s-1" ``` -`CFUnits` also provides `operator/`, `pow(int)`, `isUnknown()`, -`isDimensionless()` and an order-insensitive `operator==`. +`CFUnits` also provides `operator/`, `pow(int)`, `scaled(double)`, +`isUnknown()`, `isDimensionless()` and an order-insensitive `operator==`. A malformed units string is a programming error in a field definition, and aborting when the derived field is created stops the run at initialization diff --git a/components/omega/src/infra/CFUnits.cpp b/components/omega/src/infra/CFUnits.cpp index 25d8246cac7d..ad8abd692ab8 100644 --- a/components/omega/src/infra/CFUnits.cpp +++ b/components/omega/src/infra/CFUnits.cpp @@ -5,22 +5,63 @@ // factors, and each factor into a unit symbol and an optional integer // exponent. Factors are kept in order of first appearance so that the // formatted result reads like the conventional spelling (kg m-3, not -// m-3 kg). +// m-3 kg). A leading number is the scale factor of a scaled unit. // //===----------------------------------------------------------------------===// #include "CFUnits.h" #include "Logging.h" +#include #include +#include +#include +#include #include #include namespace OMEGA { +namespace { + +// Relative tolerance for comparing scale factors, which pick up rounding +// error in products and powers (1/1e-6 is not exactly 1e6) +constexpr double ScaleTol = 1.0e-12; + +bool sameScale(double A, double B) { + return std::abs(A - B) <= ScaleTol * std::max(std::abs(A), std::abs(B)); +} + +// Formats a scale factor with the fewest significant digits that reproduce +// it to within rounding, so 1/1e-6 is written 1e+06, not 999999.9999999999, +// and without an exponent when that is no longer (1000, not 1e+03) +std::string formatScale(double Value) { + char Buffer[32]; + int Digits = 1; + for (; Digits < 17; ++Digits) { + std::snprintf(Buffer, sizeof(Buffer), "%.*g", Digits, Value); + if (sameScale(std::strtod(Buffer, nullptr), Value)) + break; + } + std::string Shortest = Buffer; + + // %g uses an exponent once the integer part has more digits than the + // precision; asking for all of them gives the fixed spelling instead + int Exponent = + static_cast(std::floor(std::log10(std::strtod(Buffer, nullptr)))); + if (Exponent >= Digits) { + std::snprintf(Buffer, sizeof(Buffer), "%.*g", Exponent + 1, Value); + if (std::string(Buffer).size() <= Shortest.size()) + return Buffer; + } + return Shortest; +} + +} // namespace + //------------------------------------------------------------------------------ // Default constructor creates unknown units -CFUnits::CFUnits() : Known(false) {} +CFUnits::CFUnits() : Known(false), Scale(1.0) {} //------------------------------------------------------------------------------ // Parses a units string into its factors @@ -56,6 +97,7 @@ Error CFUnits::parse(const std::string &UnitsStr, // [in] units string std::string Segment; std::istringstream Segments(Trimmed); bool FirstSegment = true; + bool FirstToken = true; while (std::getline(Segments, Segment, '/')) { // A slash with nothing before or after it is malformed @@ -64,6 +106,27 @@ Error CFUnits::parse(const std::string &UnitsStr, // [in] units string bool AnyToken = false; while (Tokens >> Token) { AnyToken = true; + + // A number other than 1 is a scale factor, which udunits allows + // only before the unit symbols + bool IsNumber = std::isdigit(static_cast(Token[0])) || + Token[0] == '.'; + if (IsNumber && Token != "1") { + double Value = 1.0; + if (!FirstToken || !parseScale(Token, Value)) { + Result = CFUnits(); + RETURN_ERROR(Err, ErrorCode::Fail, + "CFUnits: cannot parse units string '{}': " + "'{}' is not a positive scale factor before the " + "unit symbols", + UnitsStr, Token); + } + Result.Scale = Value; + FirstToken = false; + continue; + } + FirstToken = false; + if (!Result.parseFactor(Token, Sign)) { Result = CFUnits(); RETURN_ERROR(Err, ErrorCode::Fail, @@ -135,6 +198,19 @@ bool CFUnits::parseFactor(const std::string &Token, // [in] one factor } // end parseFactor +//------------------------------------------------------------------------------ +// Parses a scale factor: the whole token must be a positive, finite number +bool CFUnits::parseScale(const std::string &Token, // [in] numeric token + double &Value // [out] scale factor +) { + + char *End = nullptr; + Value = std::strtod(Token.c_str(), &End); + return End == Token.c_str() + Token.size() && std::isfinite(Value) && + Value > 0.0; + +} // end parseScale + //------------------------------------------------------------------------------ // Multiplies in one symbol raised to a power void CFUnits::multiplyFactor(const std::string &Symbol, // [in] unit symbol @@ -186,6 +262,12 @@ std::string CFUnits::power(const std::string &UnitsA, // [in] units return parseOrAbort(UnitsA).pow(Exponent).str(); } +std::string CFUnits::scale(const std::string &UnitsA, // [in] units + double Factor // [in] positive scale factor +) { + return parseOrAbort(UnitsA).scaled(Factor).str(); +} + //------------------------------------------------------------------------------ // Queries bool CFUnits::isUnknown() const { return !Known; } @@ -198,10 +280,12 @@ std::string CFUnits::str() const { if (!Known) return ""; + + bool Unscaled = sameScale(Scale, 1.0); if (Factors.empty()) - return "1"; + return Unscaled ? "1" : formatScale(Scale); - std::string Result; + std::string Result = Unscaled ? "" : formatScale(Scale); for (int Pass = 0; Pass < 2; ++Pass) { for (const auto &Factor : Factors) { bool Positive = Factor.second > 0; @@ -226,6 +310,7 @@ CFUnits CFUnits::operator*(const CFUnits &Other) const { return CFUnits(); CFUnits Result = *this; + Result.Scale *= Other.Scale; for (const auto &Factor : Other.Factors) Result.multiplyFactor(Factor.first, Factor.second); return Result; @@ -242,11 +327,25 @@ CFUnits CFUnits::pow(int Exponent) const { CFUnits Result; Result.Known = true; + Result.Scale = std::pow(Scale, Exponent); for (const auto &Factor : Factors) Result.multiplyFactor(Factor.first, Factor.second * Exponent); return Result; } +CFUnits CFUnits::scaled(double Factor) const { + + OMEGA_REQUIRE(std::isfinite(Factor) && Factor > 0.0, + "CFUnits: scale factor {} is not positive and finite", Factor); + + if (!Known) + return CFUnits(); + + CFUnits Result = *this; + Result.Scale *= Factor; + return Result; +} + //------------------------------------------------------------------------------ // Comparison, ignoring the order of factors bool CFUnits::operator==(const CFUnits &Other) const { @@ -256,7 +355,7 @@ bool CFUnits::operator==(const CFUnits &Other) const { std::map Mine(Factors.begin(), Factors.end()); std::map Theirs(Other.Factors.begin(), Other.Factors.end()); - return Mine == Theirs; + return Mine == Theirs && sameScale(Scale, Other.Scale); } bool CFUnits::operator!=(const CFUnits &Other) const { diff --git a/components/omega/src/infra/CFUnits.h b/components/omega/src/infra/CFUnits.h index cc5b498484a2..786b9257f9cb 100644 --- a/components/omega/src/infra/CFUnits.h +++ b/components/omega/src/infra/CFUnits.h @@ -14,6 +14,12 @@ /// factors, supports the product, quotient and integer power of two such /// expressions, and formats the result back in the plain form. /// +/// A units string may start with a positive number, the scale factor of a +/// scaled unit in udunits (1e+06 m3 s-1 is a million cubic meters per +/// second). Rescaling a field's values, for example to report a transport in +/// units of 10^6 m3 s-1, rescales its units the other way, so CFUnits carries +/// the scale factor through every operation. +/// /// A field with no units attribute has unknown units, represented by an empty /// string. Unknown units propagate: any expression involving them is unknown, /// so a derived field never claims units its inputs did not have. Parsing @@ -69,16 +75,26 @@ class CFUnits { int Exponent ///< [in] power ); + /// Returns a units string multiplied by a positive scale factor in the CF + /// plain form (scale("m3 s-1", 1e6) is "1e+06 m3 s-1"), or an empty string + /// if it is unknown. Aborts if the string is invalid or the factor is not + /// positive and finite. + static std::string scale(const std::string &UnitsA, ///< [in] units + double Factor ///< [in] positive scale factor + ); + //--------------------------------------------------------------------------- /// True if the units are unknown (the field has no units attribute) bool isUnknown() const; - /// True if the units are known and dimensionless + /// True if the units are known and dimensionless (a pure number, which + /// may be scaled) bool isDimensionless() const; - /// Formats the units in the CF plain form: symbols with positive exponents - /// first, then those with negative exponents, each in order of first - /// appearance; "1" if dimensionless; an empty string if unknown + /// Formats the units in the CF plain form: the scale factor if it is not + /// 1, then symbols with positive exponents, then those with negative + /// exponents, each in order of first appearance; "1" if dimensionless and + /// unscaled; an empty string if unknown std::string str() const; //--------------------------------------------------------------------------- @@ -91,8 +107,12 @@ class CFUnits { /// Units expression raised to an integer power; unknown if unknown CFUnits pow(int Exponent) const; + /// Units expression multiplied by a positive scale factor; unknown if + /// unknown. Aborts if the factor is not positive and finite. + CFUnits scaled(double Factor) const; + /// True if both are unknown, or both are known with the same factors in - /// any order + /// any order and the same scale factor to within rounding bool operator==(const CFUnits &Other) const; bool operator!=(const CFUnits &Other) const; @@ -103,6 +123,9 @@ class CFUnits { /// Unit symbols and their non-zero exponents, in order of first appearance std::vector> Factors; + /// Positive scale factor multiplying the unit symbols, 1 if unscaled + double Scale; + /// Parses a units string, aborting on failure, for the string interfaces static CFUnits parseOrAbort(const std::string &UnitsStr); @@ -116,6 +139,11 @@ class CFUnits { /// the exponent. Returns false if the token is malformed. bool parseFactor(const std::string &Token, int Sign); + /// Parses a token that starts with a digit or a decimal point as a scale + /// factor. Returns false unless the whole token is a positive, finite + /// number. + static bool parseScale(const std::string &Token, double &Value); + }; // end class CFUnits } // end namespace OMEGA diff --git a/components/omega/test/infra/CFUnitsTest.cpp b/components/omega/test/infra/CFUnitsTest.cpp index 3c9092679e0f..8b4d34f321d3 100644 --- a/components/omega/test/infra/CFUnitsTest.cpp +++ b/components/omega/test/infra/CFUnitsTest.cpp @@ -5,8 +5,9 @@ /// /// Checks that units strings in the CF plain form and the udunits caret and /// slash spellings parse and format as expected, that products, quotients -/// and powers give the plain form with merged exponents, that unknown -/// (empty) units propagate, and that malformed strings are rejected. +/// and powers give the plain form with merged exponents, that scale factors +/// of scaled units are carried through, that unknown (empty) units +/// propagate, and that malformed strings are rejected. // //===----------------------------------------------------------------------===// @@ -111,6 +112,16 @@ int main(int argc, char **argv) { checkParse("kg m-3/s", "kg m-3 s-1"); checkParse("m^2/s", "m2 s-1"); + // A leading positive number is the scale factor of a scaled unit + checkParse("1e6 m3 s-1", "1e+06 m3 s-1"); + checkParse("1e+06 m3/s", "1e+06 m3 s-1"); + checkParse("2 m", "2 m"); + checkParse("0.001 kg", "0.001 kg"); + checkParse(".5 m", "0.5 m"); + checkParse("1.0 m", "m"); + checkParse("1e6", "1e+06"); + checkParse("1e3/s", "1000 s-1"); + // Malformed strings are rejected and leave unknown units checkInvalid("N m^{-2}"); checkInvalid("m*s"); @@ -119,8 +130,13 @@ int main(int argc, char **argv) { checkInvalid("m^+2"); checkInvalid("m--1"); checkInvalid("m2s"); - checkInvalid("2 m"); - checkInvalid("1e6 m3 s-1"); + checkInvalid("m 2"); + checkInvalid("m/2"); + checkInvalid("1e6 1e6 m"); + checkInvalid("0 m"); + checkInvalid("-1 m"); + checkInvalid("1e400 m"); + checkInvalid("1e6e m"); checkInvalid("m/"); checkInvalid("/s"); checkInvalid("m//s"); @@ -152,6 +168,18 @@ int main(int argc, char **argv) { checkTrue("m s-1 != m s-2", A != C); checkTrue("m s-1 != unknown", A != Unknown); checkTrue("unknown == unknown", Unknown == CFUnits()); + + // Scale factors must agree, to within rounding + CFUnits Mega, Kilo; + CFUnits::parse("1e6 m s-1", Mega); + CFUnits::parse("1e3 m s-1", Kilo); + checkTrue("1e6 m s-1 != m s-1", Mega != A); + checkTrue("1e6 m s-1 != 1e3 m s-1", Mega != Kilo); + checkTrue("1e3 m s-1 * 1e3 == 1e6 m s-1", Kilo.scaled(1.0e3) == Mega); + checkTrue("m s-1 / 1e-6 == 1e6 m s-1", A.scaled(1.0 / 1.0e-6) == Mega); + CFUnits Number; + CFUnits::parse("1e6", Number); + checkTrue("1e6 is dimensionless", Number.isDimensionless()); } // Products, quotients and powers of strings @@ -178,6 +206,22 @@ int main(int argc, char **argv) { checkResult("m2 / unknown", CFUnits::divide("m2", ""), ""); checkResult("unknown^3", CFUnits::power("", 3), ""); checkResult("unknown^0", CFUnits::power("", 0), ""); + checkResult("unknown scaled", CFUnits::scale("", 1.0e6), ""); + + // Scale factors multiply through products, quotients and powers. Values + // multiplied by 1e-6 have units scaled by 1/1e-6, which is not exactly + // 1e6 in floating point but is written as 1e+06. + checkResult("scale m3 s-1 by 1/1e-6", CFUnits::scale("m3 s-1", 1.0 / 1.0e-6), + "1e+06 m3 s-1"); + checkResult("scale m by 0.5", CFUnits::scale("m", 0.5), "0.5 m"); + checkResult("scale 1e6 m by 1e-6", CFUnits::scale("1e6 m", 1.0e-6), "m"); + checkResult("scale 1 by 1e3", CFUnits::scale("1", 1.0e3), "1000"); + checkResult("1e3 m * 1e3 s-1", CFUnits::multiply("1e3 m", "1e3 s-1"), + "1e+06 m s-1"); + checkResult("1e6 m3 / 1e3 m", CFUnits::divide("1e6 m3", "1e3 m"), "1000 m2"); + checkResult("(1e3 m)^2", CFUnits::power("1e3 m", 2), "1e+06 m2"); + checkResult("(1e3 m)^-1", CFUnits::power("1e3 m", -1), "0.001 m-1"); + checkResult("(2 m)^0", CFUnits::power("2 m", 0), "1"); // The same algebra on parsed values { From 150e2437856eff7c5caedc3cc77de5ad7596692c Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 4 Oct 2026 07:50:31 -0500 Subject: [PATCH 5/5] Derive units for the remaining analysis operators The operators added for the MOC set their output units by hand, and several got them wrong: ScalarMultiply and HorzMean wrote no units, and BinaryMultiply wrote "A * B", which is not a units string, or the units of one input when the other had none. The MOC group then overwrote the streamfunction units with "Sv", which udunits reads as the sievert. All the operators now derive their CF metadata with the AnalysisOperator helpers. BinaryMultiply multiplies its input units with CFUnits, so a velocity times an area is m3 s-1, and is unknown if either input's units are. ScalarMultiply scales the units by the reciprocal of its scalar, so the MOC streamfunction is in 1e+06 m3 s-1, and drops the standard name for a negative scalar. HorzMean inherits the units and standard name and records "area: mean". The operators that sum, mask or rescale by a dimensionless factor inherit their input's units, standard name and cell methods, and PrefixSum aborts if its boundary condition field has different units from its input. CoordinateBinning keeps its dimensionless bin index. The MOC group no longer overwrites the units the operators derive. The operator CTest checks the units every operator writes and a chain like the MOC's from a velocity and an area to Sverdrups. The Analysis guides describe the units each operator derives and those of the MOC output. Co-Authored-By: Claude Opus 5.5 (1M context) --- components/omega/doc/devGuide/Analysis.md | 16 ++++ components/omega/doc/userGuide/Analysis.md | 4 +- .../omega/src/analysis/analysisGroups/MOC.cpp | 13 +-- .../src/analysis/operators/BinaryMultiplyOp.h | 56 ++++------- .../analysis/operators/BinnedAccumulatorOp.h | 27 +++--- .../src/analysis/operators/ExtractRegionOp.h | 21 ++--- .../omega/src/analysis/operators/HorzMeanOp.h | 11 ++- .../src/analysis/operators/PrefixSumOp.h | 43 ++++++--- .../analysis/operators/PseudoToGeometricOp.h | 28 +++--- .../src/analysis/operators/ScalarMultiplyOp.h | 42 ++++++--- .../operators/TransectAccumulatorOp.h | 21 ++--- .../test/analysis/AnalysisOperatorTest.cpp | 92 +++++++++++++++++++ 12 files changed, 248 insertions(+), 126 deletions(-) diff --git a/components/omega/doc/devGuide/Analysis.md b/components/omega/doc/devGuide/Analysis.md index 0b4aec0d3b15..d369d11fb0a6 100644 --- a/components/omega/doc/devGuide/Analysis.md +++ b/components/omega/doc/devGuide/Analysis.md @@ -291,6 +291,22 @@ auto Meta = inheritMetadata(InputNames[0], Meta.Units = CFUnits::multiply(Meta.Units, "m2"); Meta.StdName.clear(); ``` +The other operators derive their units the same way: + +- `BinaryMultiply` multiplies the units of its two inputs (a velocity in + `m s-1` times an area in `m2` is `m3 s-1`) and has no standard name or cell + methods, since the product is a new quantity. +- `ScalarMultiply` scales the units by the reciprocal of its scalar with + `CFUnits::scale()`, so a flux in `m3 s-1` multiplied by `1e-6` is in + `1e+06 m3 s-1` (Sverdrups). A negative scalar also clears the standard + name, since it reverses the quantity. +- `HorzMean` inherits the units and standard name and appends `area: mean`. +- `BinnedAccumulator`, `TransectAccumulator` and `PrefixSum` sum their input, + and `ExtractRegion` and `PseudoToGeometric` multiply it by a mask or a + dimensionless factor, so all of them inherit its units. `PrefixSum` aborts + if its boundary condition field has different units from its input. +- `CoordinateBinning` outputs a bin index, which is dimensionless (`1`). + A field without units yields outputs without units, never with invented ones, and a malformed units string aborts the run when the operator is constructed. diff --git a/components/omega/doc/userGuide/Analysis.md b/components/omega/doc/userGuide/Analysis.md index d08540aeee5b..945a3f555413 100644 --- a/components/omega/doc/userGuide/Analysis.md +++ b/components/omega/doc/userGuide/Analysis.md @@ -161,7 +161,9 @@ Computes the Meridional Overturning Circulation (MOC) streamfunction using two methods: latitude-binned regional MOC and transect-based MOC. The MOC represents zonally integrated meridional mass transport as a function of latitude and depth (regional MOC) or depth alone (transect MOC). Output is in -Sverdrups (Sv), where 1 Sv = 10⁶ m³/s. +Sverdrups (Sv), where 1 Sv = 10⁶ m³/s, so its `units` attribute is +`1e+06 m3 s-1`, the udunits spelling of a Sverdrup (udunits reads `Sv` as the +sievert). **Example:** diff --git a/components/omega/src/analysis/analysisGroups/MOC.cpp b/components/omega/src/analysis/analysisGroups/MOC.cpp index 2467cec6bf10..676728cdcb51 100644 --- a/components/omega/src/analysis/analysisGroups/MOC.cpp +++ b/components/omega/src/analysis/analysisGroups/MOC.cpp @@ -121,9 +121,9 @@ MOC::MOC(const std::string &GroupName, Config &AnalysisGroupOptions, buildTemporalChains(ChainStems, AnalysisGroupOptions, AnalysisManager, ChainConfigs); - // Stamp Sverdrup units and a readable long_name on every streamfunction - // output field. Both the CF-compliant "units" and the legacy capitalized - // "Units" are set. + // Stamp a readable long_name on every streamfunction output field. The + // operators derive the units, 1e+06 m3 s-1 (Sverdrups) after the final + // ScalarMultiply(1e-6); "Sv" would be read by udunits as the sievert. for (size_t i = 0; i < ChainStems.size(); ++i) { const auto &StemStr = ChainStems[i]; const auto &RegionName = RegionList[i]; @@ -137,8 +137,6 @@ MOC::MOC(const std::string &GroupName, Config &AnalysisGroupOptions, for (const auto &SFName : SFFieldNames) { if (Field::exists(SFName)) { auto SFField = Field::get(SFName); - SFField->updateMetadata("units", std::string("Sv")); - SFField->updateMetadata("Units", std::string("Sv")); SFField->updateMetadata("long_name", LongName); } } @@ -460,15 +458,14 @@ std::string MOC::buildDepthCoordChain(const std::string &RegionName, // intermediate fields shared with other chains). AnalysisManager->parseChainAndBuildOps(ChainStr, Cfg); - // Attach descriptive depth metadata for output (once per unique field). + // Attach a descriptive long_name for output (once per unique field). The + // HorzMean operator derives the units from GeomZInterface. if (Field::exists(DepthFieldName)) { auto DepthField = Field::get(DepthFieldName); DepthField->updateMetadata( "long_name", std::string("Area-weighted representative interface depth (") + RegionName + ")"); - DepthField->updateMetadata("units", std::string("m")); - DepthField->updateMetadata("Units", std::string("m")); } return DepthFieldName; diff --git a/components/omega/src/analysis/operators/BinaryMultiplyOp.h b/components/omega/src/analysis/operators/BinaryMultiplyOp.h index f27d337c8d09..0d923ca54fdf 100644 --- a/components/omega/src/analysis/operators/BinaryMultiplyOp.h +++ b/components/omega/src/analysis/operators/BinaryMultiplyOp.h @@ -139,34 +139,29 @@ template class BinaryMultiplyOp : public AnalysisOperator { OutputNames = {OutputFieldName}; InstanceName = OutputFieldName; - // Get metadata from both input fields - std::string Descr1, Units1, StdName1; - std::string Descr2, Units2, StdName2; - + // Get descriptions from both input fields + std::string Descr1, Descr2; Field1->getMetadata("Description", Descr1); - Field1->getMetadata("Units", Units1); - Field1->getMetadata("StdName", StdName1); - Field2->getMetadata("Description", Descr2); - Field2->getMetadata("Units", Units2); - Field2->getMetadata("StdName", StdName2); + std::string OutputDescr = "Product of " + Descr1 + " and " + Descr2; - // Combine metadata for output - std::string OutputDescr = "Product of " + Descr1 + " and " + Descr2; - std::string OutputUnits = combineUnits(Units1, Units2); - std::string OutputStdName = ""; // No standard name for generic product + // A product has the product of the input units (m s-1 times m2 is + // m3 s-1), which are unknown if either input's are. It is a new + // quantity, so it has no standard name or cell methods. + InheritedMetadata Meta; + Meta.Units = CFUnits::multiply(getFieldAttribute(InputNames[0], "units"), + getFieldAttribute(InputNames[1], "units")); // Create output Field - auto OutputField = - Field::create(OutputNames[0], - OutputDescr, // Description - OutputUnits, // Units (combined) - OutputStdName, // Standard name - -std::numeric_limits::max(), // Min valid value - std::numeric_limits::max(), // Max valid value - NDims, // Rank - DimNames // Dimension names - ); + auto OutputField = createOutputField( + OutputNames[0], + OutputDescr, // Description + Meta, // CF metadata + -std::numeric_limits::max(), // Min valid value + std::numeric_limits::max(), // Max valid value + NDims, // Rank + DimNames // Dimension names + ); // Allocate output data array matching input layout OutputData = OutputArrayT(OutputNames[0] + "_out", Data1.layout()); @@ -379,21 +374,6 @@ template class BinaryMultiplyOp : public AnalysisOperator { } // end compute private: - /// Combines units from two fields for the product - std::string combineUnits(const std::string &Units1, - const std::string &Units2) { - if (Units1.empty() && Units2.empty()) { - return ""; - } else if (Units1.empty()) { - return Units2; - } else if (Units2.empty()) { - return Units1; - } else { - // Combine with multiplication notation - return Units1 + " * " + Units2; - } - } - /// Output data array holding the product field values OutputArrayT OutputData; diff --git a/components/omega/src/analysis/operators/BinnedAccumulatorOp.h b/components/omega/src/analysis/operators/BinnedAccumulatorOp.h index b55582eeea80..0f413ffce36a 100644 --- a/components/omega/src/analysis/operators/BinnedAccumulatorOp.h +++ b/components/omega/src/analysis/operators/BinnedAccumulatorOp.h @@ -125,14 +125,14 @@ template class BinnedAccumulatorOp : public AnalysisOperator { OutputNames = {OutputFieldName}; InstanceName = OutputFieldName; - // Get input metadata - std::string ValueDescr, ValueUnits, ValueStdName; + // Get input metadata. A sum over the entities in each bin has the + // units of the values summed. + std::string ValueDescr; ScalarT ValueValidMin, ValueValidMax; ValueField->getMetadata("Description", ValueDescr); - ValueField->getMetadata("Units", ValueUnits); - ValueField->getMetadata("StdName", ValueStdName); ValueField->getMetadata("ValidMin", ValueValidMin); ValueField->getMetadata("ValidMax", ValueValidMax); + auto Meta = inheritMetadata(InputNames[0]); // Create output Field dimensions std::string NumBinsDimName = "NumBins" + InputNames[1]; @@ -147,16 +147,15 @@ template class BinnedAccumulatorOp : public AnalysisOperator { } // Create output Field - auto OutputField = - Field::create(OutputNames[0], - "Binned accumulation of " + ValueDescr, // Description - ValueUnits, // Units - ValueStdName, // Standard name - ValueValidMin, // Min valid - ValueValidMax, // Max valid - InputRank, // Rank - OutputDimNames // Dimension names - ); + auto OutputField = createOutputField( + OutputNames[0], + "Binned accumulation of " + ValueDescr, // Description + Meta, // CF metadata + ValueValidMin, // Min valid + ValueValidMax, // Max valid + InputRank, // Rank + OutputDimNames // Dimension names + ); // Allocate output and local accumulation arrays if (InputRank == 1) { diff --git a/components/omega/src/analysis/operators/ExtractRegionOp.h b/components/omega/src/analysis/operators/ExtractRegionOp.h index 191c540e4ba4..d4ef6fa4c615 100644 --- a/components/omega/src/analysis/operators/ExtractRegionOp.h +++ b/components/omega/src/analysis/operators/ExtractRegionOp.h @@ -155,26 +155,25 @@ template class ExtractRegionOp : public AnalysisOperator { std::vector DimNames; InputField->getDimNames(DimNames); - // Get input metadata - std::string InputDescr, InputUnits, InputStdName; + // Get input metadata. Masking a region leaves the units, standard name + // and cell methods unchanged. + std::string InputDescr; ScalarT InputValidMin, InputValidMax; InputField->getMetadata("Description", InputDescr); - InputField->getMetadata("Units", InputUnits); - InputField->getMetadata("StdName", InputStdName); InputField->getMetadata("ValidMin", InputValidMin); InputField->getMetadata("ValidMax", InputValidMax); + auto Meta = inheritMetadata(InputNames[0]); // Create output Field with same dimensions as input - auto OutputField = Field::create( + auto OutputField = createOutputField( OutputNames[0], InputDescr + " (region: " + RegionName + ")", // Description - InputUnits, // Units (unchanged) - InputStdName, // Standard name (unchanged) - InputValidMin, // Min valid (unchanged) - InputValidMax, // Max valid (unchanged) - NDims, // Rank - DimNames // Dimension names + Meta, // CF metadata + InputValidMin, // Min valid (unchanged) + InputValidMax, // Max valid (unchanged) + NDims, // Rank + DimNames // Dimension names ); // Allocate output data array matching input layout diff --git a/components/omega/src/analysis/operators/HorzMeanOp.h b/components/omega/src/analysis/operators/HorzMeanOp.h index c438b4f14834..a2a9754a1eb3 100644 --- a/components/omega/src/analysis/operators/HorzMeanOp.h +++ b/components/omega/src/analysis/operators/HorzMeanOp.h @@ -106,13 +106,16 @@ template class HorzMeanOp : public AnalysisOperator { // Allocate 1D output data array (always Real type) OutputData = Array1DReal(OutputNames[0], NVertSize); - // Register output Field with metadata over the vertical dimension + // Register output Field with metadata over the vertical dimension. A + // horizontal mean has the units and standard name of the field it + // averages; its cell_methods record the reduction after any the input + // already carries. I4 NDims = 1; std::vector DimNames = {VertDimName}; - auto OutputField = Field::create( + auto Meta = inheritMetadata(InputNames[0], "area: mean"); + auto OutputField = createOutputField( OutputNames[0], "Area-weighted horizontal mean of " + InputNames[0], - "", // Units - "", // Standard name + Meta, // CF metadata -std::numeric_limits::max(), // Min valid value std::numeric_limits::max(), // Max valid value NDims, // Rank diff --git a/components/omega/src/analysis/operators/PrefixSumOp.h b/components/omega/src/analysis/operators/PrefixSumOp.h index 4b45f7dcaedc..f5b5c6bdf284 100644 --- a/components/omega/src/analysis/operators/PrefixSumOp.h +++ b/components/omega/src/analysis/operators/PrefixSumOp.h @@ -143,14 +143,34 @@ template class PrefixSumOp : public AnalysisOperator { OutputNames = {OutputFieldName}; InstanceName = OutputFieldName; - // Get input metadata - std::string InputDescr, InputUnits, InputStdName; + // Get input metadata. A cumulative sum has the units of the values + // summed. + std::string InputDescr; ScalarT InputValidMin, InputValidMax; InputField->getMetadata("Description", InputDescr); - InputField->getMetadata("Units", InputUnits); - InputField->getMetadata("StdName", InputStdName); InputField->getMetadata("ValidMin", InputValidMin); InputField->getMetadata("ValidMax", InputValidMax); + auto Meta = inheritMetadata(InputNames[0]); + + // The BC field seeds the sum, so it must be in the same units as the + // input. If either has unknown units, so does the sum. + if (!BCFieldName.empty()) { + CFUnits InputUnits, BCUnits; + Error Err = CFUnits::parse(Meta.Units, InputUnits); + CHECK_ERROR_ABORT(Err, "PrefixSumOp: invalid units of field {}", + InputNames[0]); + Err = CFUnits::parse(getFieldAttribute(BCFieldName, "units"), BCUnits); + CHECK_ERROR_ABORT(Err, "PrefixSumOp: invalid units of BC field {}", + BCFieldName); + if (InputUnits.isUnknown() || BCUnits.isUnknown()) { + Meta.Units = ""; + } else if (InputUnits != BCUnits) { + ABORT_ERROR("PrefixSumOp: BC field {} has units '{}' but field {} " + "has units '{}'", + BCFieldName, BCUnits.str(), InputNames[0], + InputUnits.str()); + } + } // Allow callers to override the output field name (e.g., for a short // alias) via the "OutputName" config key. @@ -158,14 +178,13 @@ template class PrefixSumOp : public AnalysisOperator { // Create output Field with same dimensions as input auto OutputField = - Field::create(OutputNames[0], - "Cumulative sum of " + InputDescr, // Description - InputUnits, // Units - InputStdName, // Standard name - InputValidMin, // Min valid - InputValidMax, // Max valid - NDims, // Rank - DimNames // Dimension names + createOutputField(OutputNames[0], + "Cumulative sum of " + InputDescr, // Description + Meta, // CF metadata + InputValidMin, // Min valid + InputValidMax, // Max valid + NDims, // Rank + DimNames // Dim names ); // Allocate output data array matching input layout diff --git a/components/omega/src/analysis/operators/PseudoToGeometricOp.h b/components/omega/src/analysis/operators/PseudoToGeometricOp.h index 138b38cd66da..7812b01c73d9 100644 --- a/components/omega/src/analysis/operators/PseudoToGeometricOp.h +++ b/components/omega/src/analysis/operators/PseudoToGeometricOp.h @@ -97,23 +97,23 @@ template class PseudoToGeometricOp : public AnalysisOperator { OutputNames = {OutputFieldName}; InstanceName = OutputFieldName; - // Get input metadata - std::string PseudoDescr, PseudoUnits, PseudoStdName; + // Get input metadata. The conversion factor RhoSw * SpecVol is + // dimensionless, so the geometric quantity has the units of the pseudo + // quantity. + std::string PseudoDescr; PseudoField->getMetadata("Description", PseudoDescr); - PseudoField->getMetadata("Units", PseudoUnits); - PseudoField->getMetadata("StdName", PseudoStdName); + auto Meta = inheritMetadata(InputNames[0]); // Create output Field with updated metadata - auto OutputField = - Field::create(OutputNames[0], - "Geometric conversion of" + PseudoDescr, // Description - PseudoUnits, // Units - PseudoStdName, // Standard name - -std::numeric_limits::max(), - std::numeric_limits::max(), - NDims, // Rank - DimNames // Dimension names - ); + auto OutputField = createOutputField( + OutputNames[0], + "Geometric conversion of" + PseudoDescr, // Description + Meta, // CF metadata + -std::numeric_limits::max(), // Min valid + std::numeric_limits::max(), // Max valid + NDims, // Rank + DimNames // Dimension names + ); // Allocate output data array matching input layout OutputData = OutputArrayT(OutputNames[0] + "_out", PseudoData.layout()); diff --git a/components/omega/src/analysis/operators/ScalarMultiplyOp.h b/components/omega/src/analysis/operators/ScalarMultiplyOp.h index c0d2183a59f1..8b13665bd440 100644 --- a/components/omega/src/analysis/operators/ScalarMultiplyOp.h +++ b/components/omega/src/analysis/operators/ScalarMultiplyOp.h @@ -32,6 +32,8 @@ #include "AnalysisOperator.h" +#include + namespace OMEGA { /// ScalarMultiplyOp multiplies all elements of an input field by a scalar @@ -65,10 +67,14 @@ template class ScalarMultiplyOp : public AnalysisOperator { // Read required scalar parameter from configuration. // Read as string first to preserve the original representation for field // naming (std::to_string loses precision, e.g. "1.0e-6" -> "0.000001"). + // The units are scaled with the scalar in double precision, so that in + // a single-precision build 1e-6 still gives units of 1e+06 m3 s-1. std::string ScalarStr; + double ScalarR8; Error Err = Options.get("Scalar", ScalarStr); if (Err.isSuccess()) { - Scalar = static_cast(std::stod(ScalarStr)); + ScalarR8 = std::stod(ScalarStr); + Scalar = static_cast(ScalarR8); } else { // Fallback: try reading directly as Real Err = Options.get("Scalar", Scalar); @@ -78,6 +84,7 @@ template class ScalarMultiplyOp : public AnalysisOperator { "in configuration"); } ScalarStr = std::to_string(Scalar); + ScalarR8 = static_cast(Scalar); } // Read optional InPlace parameter (default: false) @@ -104,24 +111,33 @@ template class ScalarMultiplyOp : public AnalysisOperator { InstanceName = OutputFieldName; // Get input metadata - std::string InputDescr, InputUnits, InputStdName; - ScalarT InputValidMin, InputValidMax, InputFillValue; + std::string InputDescr; + ScalarT InputValidMin, InputValidMax; InputField->getMetadata("Description", InputDescr); - InputField->getMetadata("StdName", InputStdName); InputField->getMetadata("ValidMin", InputValidMin); InputField->getMetadata("ValidMax", InputValidMax); + // Multiplying the values by the scalar scales the units by its + // reciprocal, so a transport in m3 s-1 multiplied by 1e-6 is in + // 1e+06 m3 s-1 (Sverdrups). A negative scalar also reverses the sign of + // the quantity, which then no longer matches the input standard name. A + // zero scalar gives zeros, which are valid in any units. + auto Meta = inheritMetadata(InputNames[0]); + if (ScalarR8 != 0) + Meta.Units = CFUnits::scale(Meta.Units, 1.0 / std::abs(ScalarR8)); + if (ScalarR8 < 0) + Meta.StdName = ""; + // Create output Field with same dimensions as input auto OutputField = - Field::create(OutputNames[0], - InputDescr + " multiplied by " + - std::to_string(Scalar), // Description - "", // Units - InputStdName, // Standard name - InputValidMin, // Min valid - InputValidMax, // Max valid - NDims, // Rank - DimNames // Dimension names + createOutputField(OutputNames[0], + InputDescr + " multiplied by " + + std::to_string(Scalar), // Description + Meta, // CF metadata + InputValidMin, // Min valid + InputValidMax, // Max valid + NDims, // Rank + DimNames // Dimension names ); // Store array size for parallel iteration diff --git a/components/omega/src/analysis/operators/TransectAccumulatorOp.h b/components/omega/src/analysis/operators/TransectAccumulatorOp.h index 79cc4fdcd79b..8ec288f0a0ab 100644 --- a/components/omega/src/analysis/operators/TransectAccumulatorOp.h +++ b/components/omega/src/analysis/operators/TransectAccumulatorOp.h @@ -144,28 +144,27 @@ class TransectAccumulatorOp : public AnalysisOperator { OutputNames = {OutputFieldName}; InstanceName = OutputFieldName; - // Get input metadata - std::string TransportDescr, TransportUnits, TransportStdName; + // Get input metadata. A signed sum of the transport through the + // transect edges has the units of the transport. + std::string TransportDescr; ScalarT TransportValidMin, TransportValidMax; TransportField->getMetadata("Description", TransportDescr); - TransportField->getMetadata("Units", TransportUnits); - TransportField->getMetadata("StdName", TransportStdName); TransportField->getMetadata("ValidMin", TransportValidMin); TransportField->getMetadata("ValidMax", TransportValidMax); + auto Meta = inheritMetadata(InputNames[0]); // Create output Field dimensions (1D: NVertLayers only) std::vector OutputDimNames = {TransportDimNames[1]}; // Create output Field - auto OutputField = Field::create( + auto OutputField = createOutputField( OutputNames[0], "Transect " + TransectName + " accumulated transport", // Description - TransportUnits, // Units - TransportStdName, // Standard name - TransportValidMin, // Min valid - TransportValidMax, // Max valid - 1, // Rank (1D) - OutputDimNames // Dimension names + Meta, // CF metadata + TransportValidMin, // Min valid + TransportValidMax, // Max valid + 1, // Rank (1D) + OutputDimNames // Dim names ); // Allocate output and local accumulation arrays (both 1D) diff --git a/components/omega/test/analysis/AnalysisOperatorTest.cpp b/components/omega/test/analysis/AnalysisOperatorTest.cpp index f9593219fb84..4a2903678767 100644 --- a/components/omega/test/analysis/AnalysisOperatorTest.cpp +++ b/components/omega/test/analysis/AnalysisOperatorTest.cpp @@ -1231,6 +1231,68 @@ void testInheritedMetadata(const MachEnv *Env, const HorzMesh *Mesh, "area: depth: maximum"); } +//------------------------------------------------------------------------------ +// Tests the units operators derive along a chain like the MOC's: a velocity +// times an area is a volume flux, and multiplying the flux by 1e-6 to get +// Sverdrups scales its units by 1e6. Also tests that a product with a field +// without units has no units, and that a negative scalar drops the standard +// name of the quantity it reverses. +void testDerivedUnits(const MachEnv *Env, const HorzMesh *Mesh, + const VertCoord *VCoord) { + + std::vector DimNames = {"NCells"}; + auto createUnitsField = [&](const std::string &Name, + const std::string &Units, + const std::string &StdName) { + auto NewField = Field::create(Name, "Derived units test field", Units, + StdName, -1.0e30, 1.0e30, 1, DimNames); + Array1DReal Data(Name + "_data", Mesh->NCellsSize); + NewField->attachData(Data); + }; + createUnitsField("TestUnitsVelocity", "m s-1", "upward_sea_water_velocity"); + createUnitsField("TestUnitsArea", "m2", "cell_area"); + createUnitsField("TestUnitsBare", "", ""); + + Config EmptyConfig; + auto scalarConfig = [](const std::string &Scalar) { + return makeOpConfig(opParam("Scalar", Scalar)); + }; + + // A velocity times an area is a volume flux with no standard name + std::string FluxName = "TestUnitsVelocity_BinaryMultiply(TestUnitsArea)"; + auto FluxOp = AnalysisOpFactory::createOp( + "BinaryMultiply", {"TestUnitsVelocity", "TestUnitsArea"}, EmptyConfig); + checkOutputMetadata("BinaryMultiply of velocity and area", FluxName, + "m3 s-1", "", ""); + + // Values multiplied by 1e-6 are in units of 1e6 m3 s-1 + auto SvOp = AnalysisOpFactory::createOp("ScalarMultiply", {FluxName}, + scalarConfig("1.0e-6")); + checkOutputMetadata("ScalarMultiply of flux to Sverdrups", + FluxName + "_ScalarMultiply(1.0e-6)", "1e+06 m3 s-1", "", + ""); + + // A scalar of 1 keeps the units and standard name; -1 keeps the units but + // reverses the quantity, so it drops the standard name + auto OneOp = AnalysisOpFactory::createOp( + "ScalarMultiply", {"TestUnitsVelocity"}, scalarConfig("1")); + checkOutputMetadata("ScalarMultiply by 1", + "TestUnitsVelocity_ScalarMultiply(1)", "m s-1", + "upward_sea_water_velocity", ""); + auto NegOp = AnalysisOpFactory::createOp( + "ScalarMultiply", {"TestUnitsVelocity"}, scalarConfig("-1")); + checkOutputMetadata("ScalarMultiply by -1", + "TestUnitsVelocity_ScalarMultiply(-1)", "m s-1", "", ""); + + // A product with a field without units has no units, rather than the + // units of the other field + auto BareOp = AnalysisOpFactory::createOp( + "BinaryMultiply", {"TestUnitsVelocity", "TestUnitsBare"}, EmptyConfig); + checkOutputMetadata("BinaryMultiply with a field without units", + "TestUnitsVelocity_BinaryMultiply(TestUnitsBare)", "", + "", ""); +} + //===----------------------------------------------------------------------===// // Main Test Functions //===----------------------------------------------------------------------===// @@ -1452,6 +1514,9 @@ void testScalarMultiplyOpType(const std::string &TypeName, const MachEnv *Env, } reportTest("ScalarMultiplyOp: " + TypeName, Passed); + checkOutputMetadata("ScalarMultiplyOp: " + TypeName, + FieldName + "_ScalarMultiply(" + ScalarStr + ")", + "0.4 m", TestStdName, ""); if (!Passed) { LOG_ERROR(" Expected: {}, Scalar: {}", static_cast(ExpectedValue), @@ -1583,6 +1648,9 @@ void testBinaryMultiplyOpSameRank(const std::string &TypeName, } reportTest("BinaryMultiplyOp (same-rank): " + TypeName, Passed); + checkOutputMetadata("BinaryMultiplyOp (same-rank): " + TypeName, + Field1Name + "_BinaryMultiply(" + Field2Name + ")", "m2", + "", ""); if (!Passed) { LOG_ERROR(" Expected: {}, Value1: {}, Value2: {}", @@ -1663,6 +1731,9 @@ void testBinaryMultiplyOpVerticalExpansion(const std::string &TypeName, } reportTest("BinaryMultiplyOp (vertical expansion): " + TypeName, Passed); + checkOutputMetadata("BinaryMultiplyOp (vertical expansion): " + TypeName, + Field2DName + "_BinaryMultiply(" + Field1DName + ")", + "m2", "", ""); if (!Passed) { LOG_ERROR(" Expected: {}, Value2D: {}, Value1D: {}", @@ -1857,6 +1928,8 @@ void testPrefixSumOpType(const std::string &TypeName, const MachEnv *Env, break; } reportTest("PrefixSumOp: " + TypeName, Passed); + checkOutputMetadata("PrefixSumOp: " + TypeName, FieldName + "_PrefixSum", + "units", "", ""); } } @@ -1989,6 +2062,9 @@ void testPrefixSumOpWithBCType(const std::string &TypeName, const MachEnv *Env, } reportTest("PrefixSumOpWithBC: " + TypeName, Passed); + checkOutputMetadata("PrefixSumOpWithBC: " + TypeName, + InputFieldName + "_PrefixSum(BC=" + BCFieldName + ")", + "m3/s", "", ""); } //------------------------------------------------------------------------------ @@ -2082,6 +2158,8 @@ void testCoordinateBinningOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("CoordinateBinningOp", Passed); + checkOutputMetadata("CoordinateBinningOp", FieldName + "_BinIndex", "1", "", + ""); } //------------------------------------------------------------------------------ @@ -2182,6 +2260,10 @@ void testBinnedAccumulatorOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("BinnedAccumulatorOp", Passed); + checkOutputMetadata("BinnedAccumulatorOp", + ValueFieldName + "_BinnedAccumulator(" + + BinIndexFieldName + ")", + "m3/s", "", ""); } //------------------------------------------------------------------------------ @@ -2343,6 +2425,9 @@ void testPseudoToGeometricOpType(const std::string &TypeName, } reportTest("PseudoToGeometricOp: " + TypeName, Passed); + checkOutputMetadata("PseudoToGeometricOp: " + TypeName, + FieldName + "_PseudoToGeometric", TestUnits, TestStdName, + ""); if (!Passed) { LOG_ERROR("PseudoToGeometricOp {} test failed after verifying {} points", @@ -2545,6 +2630,8 @@ void testExtractRegionOpType(const std::string &TypeName, const MachEnv *Env, } reportTest("ExtractRegionOp: " + TypeName, Passed); + checkOutputMetadata("ExtractRegionOp: " + TypeName, OutputNames[0], + TestUnits, TestStdName, ""); if (!Passed) { LOG_ERROR("ExtractRegionOp {} test failed after verifying {} points", @@ -2664,6 +2751,8 @@ void testHorzMeanOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("HorzMeanOp", Passed); + checkOutputMetadata("HorzMeanOp", OutputNames[0], TestUnits, TestStdName, + "area: mean"); Field::destroy(OutputNames[0]); @@ -2949,6 +3038,7 @@ void testTransectAccumulatorOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("TransectAccumulatorOp", Passed); + checkOutputMetadata("TransectAccumulatorOp", OutputNames[0], "m3/s", "", ""); } //===----------------------------------------------------------------------===// @@ -3129,6 +3219,8 @@ int main(int argc, char *argv[]) { testInheritedMetadata(DefEnv, Mesh, VCoord, ModelClock); + testDerivedUnits(DefEnv, Mesh, VCoord); + if (NumFailed > 0) { Err = 1; LOG_ERROR("AnalysisOperatorTest failure");