Add use of atlas float fields

README.md

@@ -9,7 +9,7 @@ JEDI model interface interface for the NEMO ocean model configurations on ORCA g
 
 ## Description
 
-_orca-jedi_ includes executables to calculate model values at observation locations and to perform quality control. _orca-jedi_ is a JEDI "psuedomodel", meaning that rather than interfacing directly with NEMO or NEMOVAR, the model state is derived from input files. Further applications may be developed based on _orca-jedi_ in the future. (These would be implementations of JEDI OOPS apps, such as for observation generation applications, and various DA applications).
+_orca-jedi_ includes executables to calculate model values at observation locations and to perform quality control. _orca-jedi_ is a JEDI "pseudo-model", meaning that rather than interfacing directly with NEMO or NEMOVAR, the model state is derived from input files. Further applications may be developed based on _orca-jedi_ in the future. (These would be implementations of JEDI OOPS apps, such as for observation generation applications, and various DA applications).
 
 ## Getting Started
 

examples/hofx3d_nc_orca1_checker_ice_mpi.yaml

@@ -34,7 +34,7 @@ observations:
       time interpolation: linear
       atlas-interpolator:
         type: unstructured-bilinear-lonlat
-        non_linear: missing-if-all-missing
+        non_linear: missing-if-all-missing-real32
         max_fraction_elems_to_try: 0.0
     obs operator:
       name: Identity

examples/hofx3d_odb_prof_example.yaml

@@ -42,7 +42,7 @@ observations:
     get values:
       atlas-interpolator:
         type: unstructured-bilinear-lonlat
-        non_linear: missing-if-all-missing
+        non_linear: missing-if-all-missing-real32
     obs operator:
       name: VertInterp
       observation alias file: testinput/test_name_map.yaml

examples/hofx_nc_ice_mpi.yaml

@@ -57,7 +57,7 @@ observations:
       time interpolation: linear
       atlas-interpolator:
         type: unstructured-bilinear-lonlat
-        non_linear: missing-if-all-missing
+        non_linear: missing-if-all-missing-real32
         max_fraction_elems_to_try: 0.0
     obs operator:
       name: Composite

src/orca-jedi/geometry/Geometry.cc

@@ -14,6 +14,7 @@
 #include "eckit/mpi/Comm.h"
 #include "eckit/config/Configuration.h"
 #include "eckit/exception/Exceptions.h"
+#include "eckit/system/ResourceUsage.h"
 
 #include "oops/base/Variables.h"
 #include "oops/util/DateTime.h"
@@ -44,21 +45,28 @@ Geometry::Geometry(const eckit::Configuration & config,
                    const eckit::mpi::Comm & comm) :
                       comm_(comm), vars_(orcaVariableFactory(config)),
                       n_levels_(config.getInt("number levels")),
-                      grid_(config.getString("grid name"))
+                      grid_(config.getString("grid name")),
+                      eckit_timer_(new eckit::Timer("Geometry(ORCA): ", oops::Log::trace()))
 {
+    eckit_timer_->start();
+    log_status();
     params_.validateAndDeserialize(config);
     int64_t halo = params_.sourceMeshHalo.value();
     auto meshgen_config = grid_.meshgenerator()
                           | atlas::option::halo(halo);
 
     atlas::MeshGenerator meshgen(meshgen_config);
+    log_status();
     auto partitioner_config = grid_.partitioner();
     partitioner_config.set("type",
         params_.partitioner.value());
     partitioner_ = atlas::grid::Partitioner(partitioner_config);
+    log_status();
     mesh_ = meshgen.generate(grid_, partitioner_);
+    log_status();
     funcSpace_ = atlas::functionspace::NodeColumns(
         mesh_, atlas::option::halo(halo));
+    log_status();
 }
 
 // -----------------------------------------------------------------------------
@@ -185,7 +193,7 @@ const bool Geometry::variable_in_variable_type(std::string variable_name,
   auto nemoFields = params_.nemoFields.value();
   for (const auto & nemoField : nemoFields) {
     if (nemoField.name.value() == variable_name) {
-      std::string type = nemoField.variableType.value().value_or("background");
+      std::string type = nemoField.variableType.value();
       return type == variable_type;
     }
   }
@@ -196,9 +204,31 @@ const bool Geometry::variable_in_variable_type(std::string variable_name,
   throw eckit::BadValue(err_stream.str(), Here());
 }
 
+/// \brief Data type of the atlas field holding the named variable data.
+/// \param[in]     variable_name  Name of variable.
+/// \return        orcamodel::FieldDType enum.
+FieldDType Geometry::fieldPrecision(std::string variable_name) const {
+  auto nemoFields = params_.nemoFields.value();
+  for (const auto & nemoField : nemoFields) {
+    if (nemoField.name.value() == variable_name) {
+       return nemoField.fieldPrecision.value();
+    }
+  }
+
+  std::stringstream err_stream;
+  err_stream << "orcamodel::Geometry::fieldPrecision variable name ";
+  err_stream << "\"" << variable_name << "\" not recognised. " << std::endl;
+  throw eckit::BadValue(err_stream.str(), Here());
+}
+
 void Geometry::print(std::ostream & os) const {
   os << "Not Implemented";
 }
 
+void Geometry::log_status() const {
+  oops::Log::trace() << "orcamodel::log_status " << eckit_timer_->elapsed() << " "
+      << static_cast<double>(eckit::system::ResourceUsage().maxResidentSetSize()) / 1.0e+9
+      << " Gb" << std::endl;
+}
 
 }  // namespace orcamodel

src/orca-jedi/geometry/Geometry.h

@@ -8,6 +8,7 @@
 #include <iostream>
 #include <string>
 #include <vector>
+#include <memory>
 
 #include "atlas/field/Field.h"
 #include "atlas/field/FieldSet.h"
@@ -20,19 +21,20 @@
 #include "atlas/runtime/Log.h"
 
 #include "eckit/mpi/Comm.h"
+#include "eckit/log/Timer.h"
 
 #include "oops/base/Variables.h"
 #include "oops/util/ObjectCounter.h"
 #include "oops/util/parameters/Parameters.h"
 #include "oops/util/Printable.h"
 
 #include "orca-jedi/geometry/GeometryParameters.h"
+#include "orca-jedi/utilities/Types.h"
 
 namespace atlas {
   class Field;
   class FieldSet;
   class Mesh;
-
 }
 
 namespace orcamodel {
@@ -65,6 +67,9 @@ class Geometry : public util::Printable {
   bool levelsAreTopDown() const {return true;}
   std::string distributionType() const {
       return params_.partitioner.value();}
+  FieldDType fieldPrecision(std::string variable_name) const;
+  std::shared_ptr<eckit::Timer> timer() const {return eckit_timer_;}
+  void log_status() const;
 
  private:
   void print(std::ostream &) const;
@@ -77,6 +82,7 @@ class Geometry : public util::Printable {
   atlas::Mesh mesh_;
   atlas::functionspace::NodeColumns funcSpace_;
   atlas::FieldSet nofields_;
+  std::shared_ptr<eckit::Timer> eckit_timer_;
 };
 // -----------------------------------------------------------------------------
 

src/orca-jedi/geometry/GeometryParameterTraitsFieldDType.cc

@@ -0,0 +1,13 @@
+/*
+ * (C) British Crown Copyright 2024 Met Office
+ */
+
+#include "orca-jedi/geometry/GeometryParameterTraitsFieldDType.h"
+
+namespace orcamodel {
+
+constexpr char FieldDTypeParameterTraitsHelper::enumTypeName[];
+constexpr util::NamedEnumerator<FieldDType>
+  FieldDTypeParameterTraitsHelper::namedValues[];
+
+}  // namespace orcamodel

src/orca-jedi/geometry/GeometryParameterTraitsFieldDType.h

@@ -0,0 +1,37 @@
+/*
+ * (C) British Crown Copyright 2024 Met Office
+ */
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "eckit/exception/Exceptions.h"
+#include "oops/util/parameters/ParameterTraits.h"
+
+#include "orca-jedi/utilities/Types.h"
+
+namespace orcamodel {
+
+/// Helps with the conversion of FieldDType values to/from strings.
+struct FieldDTypeParameterTraitsHelper {
+  typedef FieldDType EnumType;
+  static constexpr char enumTypeName[] = "FieldDType";
+  static constexpr util::NamedEnumerator<EnumType> namedValues[] = {
+    { EnumType::Float, "float" },
+    { EnumType::Double, "double" }
+  };
+};
+
+}  // namespace orcamodel
+
+namespace oops {
+
+/// Specialization of ParameterTraits for FieldDType.
+template <>
+struct ParameterTraits<orcamodel::FieldDType> :
+    public EnumParameterTraits<orcamodel::FieldDTypeParameterTraitsHelper>
+{};
+
+}  // namespace oops

src/orca-jedi/geometry/GeometryParameters.h

@@ -13,7 +13,8 @@
 #include "oops/base/ParameterTraitsVariables.h"
 #include "oops/util/parameters/Parameter.h"
 #include "oops/util/parameters/RequiredParameter.h"
-#include "oops/util/parameters/OptionalParameter.h"
+#include "orca-jedi/geometry/GeometryParameterTraitsFieldDType.h"
+#include "orca-jedi/utilities/Types.h"
 
 namespace orcamodel {
 
@@ -24,7 +25,14 @@ class NemoFieldParameters : public oops::Parameters {
   oops::RequiredParameter<std::string> name {"name", this};
   oops::RequiredParameter<std::string> nemoName {"nemo field name", this};
   oops::RequiredParameter<std::string> modelSpace {"model space", this};
-  oops::OptionalParameter<std::string> variableType {"variable type", this};
+  oops::Parameter<std::string> variableType {"variable type",
+    "type of variable (default is 'background' other option is 'background variance')",
+    "background",
+    this};
+  oops::Parameter<FieldDType> fieldPrecision{"field precision",
+    "Precision to store atlas fields (float (default) or double).",
+    FieldDType::Float,
+    this};
 };
 
 class OrcaGeometryParameters : public oops::Parameters {

src/orca-jedi/interpolator/Interpolator.cc

@@ -2,6 +2,7 @@
  * (C) British Crown Copyright 2024 Met Office
  */
 
+#include <cstddef>
 #include <fstream>
 #include <memory>
 #include <ostream>
@@ -111,34 +112,68 @@ namespace orcamodel {
       nvals += nlocs_ * varSizes[jvar];
       oops::Log::debug() << " varSizes[" << jvar << "] = " << varSizes[jvar];
     }
-    result.resize(nvals);
     oops::Log::debug() << " nvals = " << nvals << std::endl;
+    result.resize(nvals);
 
-    std::size_t out_idx = 0;
+    auto res_iter = result.begin();
     for (size_t jvar=0; jvar < nvars; ++jvar) {
-      auto gv_varname = vars[jvar];
-      atlas::Field tgt_field = atlasObsFuncSpace_.createField<double>(
-          atlas::option::name(gv_varname) |
-          atlas::option::levels(varSizes[jvar]));
-      interpolator_.execute(state.stateFields()[gv_varname], tgt_field);
-      auto field_view = atlas::array::make_view<double, 2>(tgt_field);
-      atlas::field::MissingValue mv(state.stateFields()[gv_varname]);
-      bool has_mv = static_cast<bool>(mv);
-      for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) {
-        for (std::size_t iloc=0; iloc < nlocs_; iloc++) {
-          if (has_mv && mv(field_view(iloc, klev))) {
-            result[out_idx] = util::missingValue<double>();
-          } else {
-            result[out_idx] = field_view(iloc, klev);
-          }
-          ++out_idx;
-        }
+      switch (state.geometry()->fieldPrecision(vars[jvar])) {
+        case FieldDType::Double:
+          executeInterpolation<double>(vars[jvar], varSizes[jvar], state, res_iter);
+          break;
+        case FieldDType::Float:
+          executeInterpolation<float>(vars[jvar], varSizes[jvar], state, res_iter);
+          break;
+        default:
+          throw eckit::BadParameter("orcamodel::Interpolator::apply '"
+              + vars[jvar] + "' field type not recognised This line should never run!");
       }
     }
+    assert(result.size() == nvals);
     oops::Log::trace() << "orcamodel::Interpolator::apply done "
                        << std::endl;
   }
 
+  /// \brief Execute the atlas interpolation on a single field.
+  /// \param gv_varname The GeoVaLs variable name of the source field.
+  /// \param var_size The number of observation locations in the interpolant.
+  /// \param state The state object containing the model state.
+  /// \param iter Reference to the interator into the output vector.
+  template<class T> void Interpolator::executeInterpolation(
+      const std::string& gv_varname,
+      size_t var_size,
+      const State& state,
+      std::vector<double>::iterator& iter) const {
+    atlas::Field tgt_field = atlasObsFuncSpace_.createField<T>(
+        atlas::option::name(gv_varname) |
+        atlas::option::levels(var_size));
+    interpolator_.execute(state.stateFields()[gv_varname], tgt_field);
+    auto field_view = atlas::array::make_view<T, 2>(tgt_field);
+    atlas::field::MissingValue mv(state.stateFields()[gv_varname]);
+    bool has_mv = static_cast<bool>(mv);
+    for (std::size_t klev=0; klev < var_size; ++klev) {
+      for (std::size_t iloc=0; iloc < nlocs_; iloc++) {
+        if (has_mv && mv(field_view(iloc, klev))) {
+          *iter = util::missingValue<double>();
+        } else {
+          *iter = static_cast<double>(field_view(iloc, klev));
+        }
+        ++iter;
+      }
+    }
+  }
+
+  template void Interpolator::executeInterpolation<double>(
+      const std::string& gv_varname,
+      size_t var_size,
+      const State& state,
+      std::vector<double>::iterator& iter) const;
+  template void Interpolator::executeInterpolation<float>(
+      const std::string& gv_varname,
+      size_t var_size,
+      const State& state,
+      std::vector<double>::iterator& iter) const;
+
   void Interpolator::print(std::ostream & os) const {
     os << "orcamodel::Interpolator: " << std::endl;
     os << "  Obs function space " << atlasObsFuncSpace_ << std::endl;

src/orca-jedi/interpolator/Interpolator.h

@@ -70,6 +70,11 @@ class Interpolator : public util::Printable,
   }
 
  private:
+  template<class T> void executeInterpolation(
+      const std::string& gv_varname,
+      size_t var_size,
+      const State& state,
+      std::vector<double>::iterator& result) const;
   void print(std::ostream &) const override;
   int64_t nlocs_;
   atlas::functionspace::PointCloud atlasObsFuncSpace_;