Add 3DVar and error covariance toolbox (dirac test) main programs. (#106

)

* Add 3DVar and error covariance toolbox (dirac test) main programs.
* Add interpolator increment apply and adjoint methods and add increment to
state method. 
* Update ci container (#117)
* Add reference data for 3dvar and dirac application tests. Turn off
creation of json schema again for 3dvar and errorcovariancetoolbox
applications.

Co-authored-by: Toby Searle <[email protected]>

.github/workflows/ci.yml

@@ -5,9 +5,10 @@ on:
     branches: [develop]
   pull_request:
     branches: [develop]
+  workflow_dispatch:
 env:
     REGISTRY: ghcr.io
-    IMAGE_NAME: twsearle/orca-jedi/ci-almalinux9:v1.3.0
+    IMAGE_NAME: twsearle/orca-jedi/ci-almalinux9:v1.4.0
 jobs:
   build:
     runs-on: ubuntu-latest
@@ -32,6 +33,20 @@ jobs:
           repository: JCSDA-internal/oops
           token: ${{ secrets.GHCR_PAT }}
 
+      - name: checkout vader
+        uses: actions/checkout@v3
+        with:
+          path: ci/vader
+          repository: JCSDA-internal/vader
+          token: ${{ secrets.GHCR_PAT }}
+
+      - name: checkout saber
+        uses: actions/checkout@v3
+        with:
+          path: ci/saber
+          repository: JCSDA-internal/saber
+          token: ${{ secrets.GHCR_PAT }}
+
       - name: checkout ioda
         uses: actions/checkout@v3
         with:

CMakeLists.txt

@@ -36,6 +36,9 @@ ecbuild_debug( "   atlas_orca_FEATURES: [${atlas_orca_FEATURES}]" )
 find_package( oops REQUIRED )
 ecbuild_debug( "   oops_FEATURES : [${oops_FEATURES}]" )
 
+find_package( saber REQUIRED)
+ecbuild_debug( "   saber_FEATURES : [${saber_FEATURES}]" )
+
 find_package( ufo REQUIRED)
 ecbuild_debug( "   ufo_FEATURES : [${ufo_FEATURES}]" )
 

VERSION

@@ -1 +1 @@
-1.3.0
+1.4.0

ci/CMakeLists.txt

@@ -18,6 +18,8 @@ if(NOT DEFINED jedicmake_DIR)
 endif()
 
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/oops" EXCLUDE_FROM_ALL)
+add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/vader" EXCLUDE_FROM_ALL)
+add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/saber" EXCLUDE_FROM_ALL)
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/ioda" EXCLUDE_FROM_ALL)
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/ufo" EXCLUDE_FROM_ALL)
 add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/orca-jedi")

ci/hpccm_recipe_almalinux9.py

@@ -50,6 +50,7 @@ def gitlab_url(repo, vn):
 odc_vn = USERARG.get('odc_vn', '1.5.2')
 openmpi_vn = USERARG.get('openmpi_vn', '4.1.5')
 pycodestyle_vn = USERARG.get('pycodestyle_vn', '2.10')
+qhull_vn = USERARG.get('qhull_vn', '8.0.2')  # no qhull-devel rpm
 udunits_vn = USERARG.get('udunits_vn', '2.2.28')
 yaxt_vn = USERARG.get('yaxt_vn', '528-0.10.0')  # URL has a number and a version
 
@@ -203,6 +204,13 @@ def gitlab_url(repo, vn):
     cmake_opts=['-DCMAKE_BUILD_TYPE=Release', '-DJSON_BuildTests=OFF'],
 )
 
+Stage0 += generic_cmake(
+    prefix='/usr/local',
+    url=github_url('qhull/qhull', f'v{qhull_vn}'),
+    directory=f'qhull-{qhull_vn}',
+    cmake_opts=['-DCMAKE_BUILD_TYPE=Release'],
+)
+
 Stage0 += generic_cmake(
     prefix='/usr/local',
     url=github_url('NOAA-EMC/NCEPLIBS-bufr', f'v{nceplibs_bufr_vn}'),

src/mains/CMakeLists.txt

@@ -11,6 +11,12 @@ if( ${nemo-feedback_FOUND} )
                           LIBS    orcamodel nemo_feedback
                           DEFINITIONS NEMO_FEEDBACK_EXISTS
                         )
+
+  ecbuild_add_executable( TARGET  orcamodel_3DVar.x
+                          SOURCES orcamodel3DVar.cc
+                          LIBS    orcamodel nemo_feedback
+                          DEFINITIONS NEMO_FEEDBACK_EXISTS
+                        )
 else()
   message("      compiling without feedback file support")
   ecbuild_add_executable( TARGET  orcamodel_hofx.x
@@ -22,6 +28,21 @@ else()
                           SOURCES orcamodelHofX3D.cc
                           LIBS    orcamodel
                         )
+
+  ecbuild_add_executable( TARGET  orcamodel_3DVar.x
+                          SOURCES orcamodel3DVar.cc
+                          LIBS    orcamodel
+                        )
 endif()
+
+ecbuild_add_executable( TARGET  orcamodel_ErrorCovarianceToolbox.x
+                        SOURCES orcamodelErrorCovarianceToolbox.cc
+                        LIBS    orcamodel
+                      )
+
 oops_output_json_schema( "orcamodel_hofx.x" )
 oops_output_json_schema( "orcamodel_hofx3D.x" )
+# The applications below throw an error generating json schema.
+# This appears to be a generic issue see also oops/qg/mains/CMakeLists.txt
+#oops_output_json_schema( "orcamodel_3DVar.x" )
+#oops_output_json_schema( "orcamodel_ErrorCovarianceToolbox.x" )

src/mains/orcamodel3DVar.cc

@@ -0,0 +1,36 @@
+/*
+ * (C) British Crown Copyright 2024 MetOffice
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#include "atlas/library/Library.h"
+#include "oops/generic/instantiateModelFactory.h"
+
+#include "oops/runs/Run.h"
+#include "oops/runs/Variational.h"
+#include "saber/oops/instantiateCovarFactory.h"
+
+#include "ufo/ObsTraits.h"
+#include "ufo/instantiateObsFilterFactory.h"
+#if defined(NEMO_FEEDBACK_EXISTS)
+  #include "nemo-feedback/instantiateObsFilterFactory.h"
+#endif
+
+#include "orca-jedi/utilities/OrcaModelTraits.h"
+
+int main(int argc,  char ** argv) {
+  oops::Run run(argc, argv);
+  oops::instantiateModelFactory<orcamodel::OrcaModelTraits>();
+  atlas::Library::instance().initialise();
+  saber::instantiateCovarFactory<orcamodel::OrcaModelTraits>();
+  ufo::instantiateObsFilterFactory();
+#if defined(NEMO_FEEDBACK_EXISTS)
+  nemo_feedback::instantiateObsFilterFactory<ufo::ObsTraits>();
+#endif
+  oops::Variational<orcamodel::OrcaModelTraits , ufo::ObsTraits> var;
+  int i = run.execute(var);
+  atlas::Library::instance().finalise();
+  return i;
+}

src/mains/orcamodelErrorCovarianceToolbox.cc

@@ -0,0 +1,25 @@
+/*
+ * (C) British Crown Copyright 2024 MetOffice
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ */
+
+#include "atlas/library/Library.h"
+
+#include "saber/oops/ErrorCovarianceToolbox.h"
+#include "saber/oops/instantiateCovarFactory.h"
+
+#include "oops/runs/Run.h"
+#include "orca-jedi/utilities/OrcaModelTraits.h"
+
+int main(int argc,  char ** argv) {
+  oops::Run run(argc, argv);
+  atlas::Library::instance().initialise();
+  saber::instantiateCovarFactory<orcamodel::OrcaModelTraits>();
+
+  saber::ErrorCovarianceToolbox<orcamodel::OrcaModelTraits> var;
+  int i = run.execute(var);
+  atlas::Library::instance().finalise();
+  return i;
+}

src/orca-jedi/CMakeLists.txt

@@ -35,15 +35,20 @@ utilities/Types.h
 utilities/Types.cc
 utilities/IOUtils.h
 utilities/IOUtils.cc
+utilities/ModelData.h
+utilities/ModelData.cc
 variablechanges/VariableChangeParameters.h
 variablechanges/VariableChange.h
+variablechanges/LinearVariableChange.h
+variablechanges/LinearVariableChange.cc
+variablechanges/LinearVariableChangeParameters.h
 )
 
 
 
 ecbuild_add_library( TARGET   orcamodel
                      SOURCES  ${oops_orcamodel_src_files}
-                     PUBLIC_LIBS  oops ufo atlas atlas-orca eckit netcdf
+                     PUBLIC_LIBS  oops saber ufo atlas atlas-orca eckit netcdf
                      INSTALL_HEADERS LISTED
                      LINKER_LANGUAGE ${OOPS_LINKER_LANGUAGE}
                     )

src/orca-jedi/increment/Increment.cc

@@ -207,7 +207,7 @@ Increment & Increment::operator*=(const double & zz) {
   return *this;
 }
 
-/// \brief Create increment from the different of two state objects.
+/// \brief Create increment from the difference of two state objects.
 /// \param x1 State object.
 /// \param x2 State object subtracted.
 void Increment::diff(const State & x1, const State & x2) {
@@ -222,6 +222,12 @@ void Increment::diff(const State & x1, const State & x2) {
     atlas::Field field2 = x2.getField(i);
     atlas::Field fieldi = incrementFields_[i];
 
+    atlas::field::MissingValue mv1(field1);
+    bool has_mv1 = static_cast<bool>(mv1);
+    atlas::field::MissingValue mv2(field2);
+    bool has_mv2 = static_cast<bool>(mv2);
+    bool has_mv = has_mv1 || has_mv2;
+
     std::string fieldName1 = field1.name();
     std::string fieldName2 = field2.name();
     std::string fieldNamei = fieldi.name();
@@ -234,8 +240,15 @@ void Increment::diff(const State & x1, const State & x2) {
     auto field_viewi = atlas::array::make_view<double, 2>(fieldi);
     for (atlas::idx_t j = 0; j < field_viewi.shape(0); ++j) {
       for (atlas::idx_t k = 0; k < field_viewi.shape(1); ++k) {
-        if (!ghost(j)) { field_viewi(j, k) = field_view1(j, k) - field_view2(j, k);
-        } else { field_viewi(j, k) = 0; }
+        if (!ghost(j)) {
+          if (!has_mv1 || (has_mv1 && !mv1(field_view1(j, k)))) {
+            if (!has_mv2 || (has_mv2 && !mv2(field_view2(j, k)))) {
+              field_viewi(j, k) = field_view1(j, k) - field_view2(j, k);
+            }
+          }
+        } else {
+          field_viewi(j, k) = 0;
+        }
       }
     }
   }
@@ -532,7 +545,7 @@ void Increment::read(const eckit::Configuration & conf) {
   throw eckit::NotImplemented(err_message, Here());
 }
 
-/// \brief write out increments fields to a file using params specified filename.
+/// \brief Write out increments fields to a file using params specified filename.
 void Increment::write(const OrcaIncrementParameters & params) const {
   oops::Log::debug() << "orcamodel::increment::write" << std::endl;
 
@@ -545,10 +558,12 @@ void Increment::write(const OrcaIncrementParameters & params) const {
   oops::Log::debug() << "Increment::write to filename "
                      << nemo_field_path << std::endl;
 
+incrementFields_.haloExchange();
+
   writeFieldsToFile(nemo_field_path, *geom_, time_, incrementFields_);
 }
 
-/// \brief write out increments fields to a file using config specified filename.
+/// \brief Write out increments fields to a file using config specified filename.
 void Increment::write(const eckit::Configuration & config) const {
   write(oops::validateAndDeserialize<OrcaIncrementParameters>(config));
 }

src/orca-jedi/interpolator/Interpolator.cc

@@ -27,6 +27,7 @@
 
 #include "orca-jedi/geometry/Geometry.h"
 #include "orca-jedi/state/State.h"
+#include "orca-jedi/increment/Increment.h"
 
 #include "orca-jedi/interpolator/Interpolator.h"
 
@@ -177,6 +178,128 @@ template void Interpolator::executeInterpolation<float>(
     const std::vector<bool> & mask,
     std::vector<double>::iterator& iter) const;
 
+/// \brief Interpolate from model space to observation space
+/// \param vars Oops variables
+/// \param inc Increment object (input)
+/// \param mask Mask vector in observation space
+/// \param result Result (output) vector in observation space
+void Interpolator::apply(const oops::Variables& vars, const Increment& inc,
+           const std::vector<bool> & mask,
+           std::vector<double>& result) const {
+  // input is inc output is result
+  const size_t nvars = vars.size();
+
+  for (size_t j=0; j < nvars; ++j) {
+    if (!inc.variables().has(vars[j])) {
+      std::stringstream err_stream;
+      err_stream << "orcamodel::Interpolator::apply varname \" "
+                 << "\" " << vars[j]
+                 << " not found in the model increment." << std::endl;
+      err_stream << "    add the variable to the increment variables and "
+                 << "add a mapping from the geometry to that variable."
+                 << std::endl;
+      throw eckit::BadParameter(err_stream.str(), Here());
+    }
+  }
+
+  const std::vector<size_t> varSizes =
+    inc.geometry()->variableSizes(vars);
+  size_t nvals = 0;
+  for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar];
+  result.resize(nvals);
+
+  std::size_t out_idx = 0;
+  for (size_t jvar=0; jvar < nvars; ++jvar) {
+    auto gv_varname = vars[jvar].name();
+    atlas::Field tgt_field = atlasObsFuncSpace_.createField<double>(
+        atlas::option::name(gv_varname) |
+        atlas::option::levels(varSizes[jvar]));
+    interpolator_.execute(inc.incrementFields()[gv_varname], tgt_field);
+    auto field_view = atlas::array::make_view<double, 2>(tgt_field);
+    atlas::field::MissingValue mv(inc.incrementFields()[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;
+      }
+    }
+  }
+}
+
+/// \brief Interpolate from observation space to model space
+/// \param vars Oops variables
+/// \param inc Increment object (output)
+/// \param mask Mask (observation space) vector
+/// \param resultin Values (observation space) vector (input)
+void Interpolator::applyAD(const oops::Variables& vars, Increment& inc,
+             const std::vector<bool> & mask,
+             const std::vector<double> & resultin) const
+{
+  // input is resultin output is inc
+
+  oops::Log::trace() << "orcamodel::Interpolator::applyAD start "
+                     << std::endl;
+
+  const size_t nvars = vars.size();
+
+  for (size_t j=0; j < nvars; ++j) {
+    if (!inc.variables().has(vars[j])) {
+      std::stringstream err_stream;
+      err_stream << "orcamodel::Interpolator::apply varname \" "
+                 << "\" " << vars[j]
+                 << " not found in the model increment." << std::endl;
+      err_stream << "    add the variable to the increment variables and "
+                 << "add a mapping from the geometry to that variable."
+                 << std::endl;
+      throw eckit::BadParameter(err_stream.str(), Here());
+    }
+  }
+
+  const std::vector<size_t> varSizes =
+    inc.geometry()->variableSizes(vars);
+  size_t nvals = 0;
+
+  for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar];
+
+  std::size_t out_idx = 0;
+  for (size_t jvar=0; jvar < nvars; ++jvar) {
+    auto gv_varname = vars[jvar].name();
+    auto tgt_field = atlasObsFuncSpace_.createField<double>(
+      atlas::option::name(gv_varname) |
+      atlas::option::levels(varSizes[jvar]));
+
+// Copying observation array vector to an atlas observation field (tgt_field)
+    auto field_view = atlas::array::make_view<double, 2>(tgt_field);
+    atlas::field::MissingValue mv(inc.incrementFields()[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))) {
+             field_view(iloc, klev) = util::missingValue<double>();
+          } else {
+             field_view(iloc, klev) = resultin[out_idx];
+          }
+        ++out_idx;
+      }
+    }
+
+// halo exchange update ghost points
+    std::shared_ptr<const Geometry> geom = inc.geometry();
+    geom->functionSpace().haloExchange(inc.incrementFields()[gv_varname]);
+
+    interpolator_.execute_adjoint(inc.incrementFields()[gv_varname], tgt_field);
+  }    // jvar
+
+  oops::Log::trace() << "orcamodel::Interpolator::applyAD done "
+                     << std::endl;
+}
+
 void Interpolator::print(std::ostream & os) const {
   os << "orcamodel::Interpolator: " << std::endl;
   os << "  Obs function space " << atlasObsFuncSpace_ << std::endl;