diff --git a/components/omega/doc/devGuide/Analysis.md b/components/omega/doc/devGuide/Analysis.md index 6199fb3b7bfb..d369d11fb0a6 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,69 @@ 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(); +``` +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. + ## Operator factory and type dispatch Fields carry runtime metadata describing their scalar type, rank, and memory @@ -493,9 +558,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/devGuide/CFUnits.md b/components/omega/doc/devGuide/CFUnits.md new file mode 100644 index 000000000000..407ad51de2f3 --- /dev/null +++ b/components/omega/doc/devGuide/CFUnits.md @@ -0,0 +1,58 @@ +(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 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. + +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)`, `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 +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/doc/userGuide/Analysis.md b/components/omega/doc/userGuide/Analysis.md index e760f15b3e5f..945a3f555413 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:** @@ -157,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/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/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/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/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/src/infra/CFUnits.cpp b/components/omega/src/infra/CFUnits.cpp new file mode 100644 index 000000000000..ad8abd692ab8 --- /dev/null +++ b/components/omega/src/infra/CFUnits.cpp @@ -0,0 +1,367 @@ +//===-- 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). 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), Scale(1.0) {} + +//------------------------------------------------------------------------------ +// 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; + bool FirstToken = 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; + + // 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, + "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 + +//------------------------------------------------------------------------------ +// 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 + 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(); +} + +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; } + +bool CFUnits::isDimensionless() const { return Known && Factors.empty(); } + +//------------------------------------------------------------------------------ +// Formats the units in the CF plain form +std::string CFUnits::str() const { + + if (!Known) + return ""; + + bool Unscaled = sameScale(Scale, 1.0); + if (Factors.empty()) + return Unscaled ? "1" : formatScale(Scale); + + std::string Result = Unscaled ? "" : formatScale(Scale); + 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; + Result.Scale *= Other.Scale; + 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; + 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 { + + 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 && sameScale(Scale, Other.Scale); +} + +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..786b9257f9cb --- /dev/null +++ b/components/omega/src/infra/CFUnits.h @@ -0,0 +1,152 @@ +#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 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 +/// 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 + ); + + /// 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 (a pure number, which + /// may be scaled) + bool isDimensionless() const; + + /// 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; + + //--------------------------------------------------------------------------- + /// 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; + + /// 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 and the same scale factor to within rounding + 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; + + /// 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); + + /// 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); + + /// 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 + +//===----------------------------------------------------------------------===// +#endif // OMEGA_CFUNITS_H 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/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/analysis/AnalysisOperatorTest.cpp b/components/omega/test/analysis/AnalysisOperatorTest.cpp index f324147d61d0..4a2903678767 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,117 @@ 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"); +} + +//------------------------------------------------------------------------------ +// 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 //===----------------------------------------------------------------------===// @@ -1335,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), @@ -1466,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: {}", @@ -1546,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: {}", @@ -1740,6 +1928,8 @@ void testPrefixSumOpType(const std::string &TypeName, const MachEnv *Env, break; } reportTest("PrefixSumOp: " + TypeName, Passed); + checkOutputMetadata("PrefixSumOp: " + TypeName, FieldName + "_PrefixSum", + "units", "", ""); } } @@ -1872,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", "", ""); } //------------------------------------------------------------------------------ @@ -1965,6 +2158,8 @@ void testCoordinateBinningOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("CoordinateBinningOp", Passed); + checkOutputMetadata("CoordinateBinningOp", FieldName + "_BinIndex", "1", "", + ""); } //------------------------------------------------------------------------------ @@ -2065,6 +2260,10 @@ void testBinnedAccumulatorOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("BinnedAccumulatorOp", Passed); + checkOutputMetadata("BinnedAccumulatorOp", + ValueFieldName + "_BinnedAccumulator(" + + BinIndexFieldName + ")", + "m3/s", "", ""); } //------------------------------------------------------------------------------ @@ -2226,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", @@ -2428,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", @@ -2547,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]); @@ -2832,6 +3038,7 @@ void testTransectAccumulatorOp(const MachEnv *Env, const HorzMesh *Mesh, } reportTest("TransectAccumulatorOp", Passed); + checkOutputMetadata("TransectAccumulatorOp", OutputNames[0], "m3/s", "", ""); } //===----------------------------------------------------------------------===// @@ -3010,6 +3217,10 @@ int main(int argc, char *argv[]) { testTransectAccumulatorOp(DefEnv, Mesh, VCoord); + testInheritedMetadata(DefEnv, Mesh, VCoord, ModelClock); + + testDerivedUnits(DefEnv, Mesh, VCoord); + if (NumFailed > 0) { Err = 1; LOG_ERROR("AnalysisOperatorTest failure"); diff --git a/components/omega/test/infra/CFUnitsTest.cpp b/components/omega/test/infra/CFUnitsTest.cpp new file mode 100644 index 000000000000..8b4d34f321d3 --- /dev/null +++ b/components/omega/test/infra/CFUnitsTest.cpp @@ -0,0 +1,248 @@ +//===-- 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 scale factors +/// of scaled units are carried through, 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"); + + // 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"); + checkInvalid("m.s"); + checkInvalid("m^"); + checkInvalid("m^+2"); + checkInvalid("m--1"); + checkInvalid("m2s"); + 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"); + 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()); + + // 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 + 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), ""); + 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 + { + 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; +} 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()) {