test/ufo/GeoVaLs.h

Go to the documentation of this file.

/*
 * (C) Copyright 2019 UK Met Office
 *
 * 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.
 */
 
#ifndef TEST_UFO_GEOVALS_H_
#define TEST_UFO_GEOVALS_H_
 
#include <cmath>
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
 
#define ECKIT_TESTING_SELF_REGISTER_CASES 0
 
#include "eckit/config/LocalConfiguration.h"
#include "eckit/testing/Test.h"
#include "ioda/ObsSpace.h"
#include "oops/mpi/mpi.h"
#include "oops/runs/Test.h"
#include "oops/util/FloatCompare.h"
#include "oops/util/Logger.h"
#include "test/TestEnvironment.h"
#include "ufo/GeoVaLs.h"
#include "ufo/Locations.h"
#include "ufo/ObsOperator.h"
 
namespace ufo {
namespace test {
 
// -----------------------------------------------------------------------------
 
void testGeoVaLs() {
  const eckit::LocalConfiguration conf(::test::TestEnvironment::config());
  util::DateTime bgn(conf.getString("window begin"));
  util::DateTime end(conf.getString("window end"));
 
  std::vector<eckit::LocalConfiguration> confs;
  conf.get("geovals test", confs);
  for (size_t jconf = 0; jconf < confs.size(); ++jconf) {
/// Setup ObsSpace
    const eckit::LocalConfiguration obsconf(confs[jconf], "obs space");
    ioda::ObsSpace ospace(obsconf, oops::mpi::world(), bgn, end, oops::mpi::myself());
 
/// Setup GeoVaLs
    const eckit::LocalConfiguration gconf(confs[jconf], "geovals");
    const oops::Variables ingeovars(gconf, "state variables");
    const GeoVaLs gval(gconf, ospace, ingeovars);
 
    const double tol = gconf.getDouble("tolerance");
 
/// Check that GeoVaLs default constructor works
     oops::Log::trace() <<
      "GeoVaLs default constructor - does not allocate fields" << std::endl;
    GeoVaLs gv_temp(ospace.comm());
 
/// Check that GeoVaLs constructor to create a GeoVaLs with one location works
    if (gconf.has("one location check")) {
      oops::Log::trace() << "Check that GeoVaLs constructor for one location works" << std::endl;
      const eckit::LocalConfiguration gconfone(gconf, "one location check");
      const std::string var = gconfone.getString("variable");
      const std::vector<int> ind = gconfone.getIntVector("indices");
      const std::vector<float> values = gconfone.getFloatVector("values");
      const float oneloctol = gconfone.getFloat("tolerance");
 
      // Loop over each location and test just the lowest level
      oops::TestVerbosity verbosity = oops::TestVerbosity::LOG_SUCCESS_AND_FAILURE;
      for (std::size_t i = 0; i < ind.size(); ++i) {
        GeoVaLs gv_one(gval, ind[i]);
        std::vector<float> gv_val(1);
        gv_one.get(gv_val, var, 1);
        EXPECT(oops::is_close_absolute(gv_val[0], values[i], oneloctol, verbosity));
      }
    } else {
      oops::Log::trace() << "Test just the constructor for a one location GeoVaLs" << std::endl;
      int index = 0;
      GeoVaLs gv_one(gval, index);
    }
 
    GeoVaLs gv(gval);
    if (gconf.has("reorderzdir check")) {
      const eckit::LocalConfiguration gconfchk(gconf, "reorderzdir check");
      const std::string flipto = gconfchk.getString("direction");
      std::string flipback = (flipto == "bottom2top") ? "top2bottom" : "bottom2top";
      std::size_t nobs = ospace.nlocs();
      gv.reorderzdir("air_pressure_levels", flipto);
      std::vector<float> gvar(nobs);
      std::vector<float> gvarref(nobs);
      float sum;
      for (size_t i = 0; i < ingeovars.size(); ++i) {
        size_t nlevs = gv.nlevs(ingeovars[i]);
        sum = 0;
        for (size_t k = 0; k < nlevs; ++k) {
          size_t kk = nlevs - k;
          gv.get(gvar, ingeovars[i], k+1);
          gval.get(gvarref, ingeovars[i], kk);
          for (size_t  j = 0; j < nobs; ++j) {
            gvar[j] = gvar[j] - gvarref[j];
            sum += sum + gvar[j];
          }
        }
      }
      gv.reorderzdir("air_pressure_levels", flipback);
      const double tol = gconfchk.getDouble("tolerance");
      EXPECT(abs(sum) < tol);
    }
 
/// Check that GeoVaLs merge followed by a split gives back the original geovals
    oops::Log::trace() <<
      "GeoVaLs merge followed by a split gives back the original geovals" << std::endl;
 
    double dp_gval = gval.dot_product_with(gval);
    oops::Log::debug()<< "initial gval dot product with itself " << dp_gval << std::endl;
 
    gv.zero();
    gv.merge(gval, gval);
 
    double dp_gv_merged = gv.dot_product_with(gv);
    oops::Log::debug()<< "gv dot product with itself after merge " << dp_gv_merged << std::endl;
    EXPECT(abs(dp_gv_merged - 2.0 * dp_gval)/dp_gv_merged < tol);
 
    GeoVaLs gv1(ospace.comm()), gv2(ospace.comm());
    gv.split(gv1, gv2);
 
    double dp_gv1_split = gv1.dot_product_with(gv1);
    double dp_gv2_split = gv2.dot_product_with(gv2);
 
    oops::Log::debug()<< "gv1 gv2 dot products with itself after split "
                      << dp_gv1_split << " " << dp_gv2_split << std::endl;
 
    EXPECT(gv1.rms() == gv2.rms());
    EXPECT(gv1.rms() == gval.rms());
    EXPECT(abs(dp_gv1_split - dp_gv2_split)/dp_gv1_split < tol);
    EXPECT(abs(dp_gv1_split - dp_gval)/dp_gv1_split < tol);
 
    oops::Log::trace() <<
      "GeoVaLs merge followed by a split test succeeded" << std::endl;
 
///  Check that  GeoVaLs & operator *= (const std::vector<float>);
    oops::Log::trace() <<
      "Check that GeoVaLs & operator *= (const std::vector<float>);" << std::endl;
 
    std::size_t nlocs = ospace.nlocs();
    {
      std::vector<float> tw(nlocs, 1.0f);
      gv1 *= tw;
      EXPECT(gv1.rms() == gval.rms());
    }
    {
      std::vector<float> tw(nlocs, 2.0f);
      gv1 *= tw;
      double rms1, rms2;
      rms1 = gv1.rms();
      rms2 = 2.0 * gval.rms();
      oops::Log::debug()<< "rms1, rms2 = " <<  rms1  << "  " << rms2 << std::endl;
      EXPECT(rms1 == rms2);
    }
    oops::Log::trace() <<
      "GeoVaLs & operator *= (const std::vector<float>); test succeeded" << std::endl;
  }
}
// -----------------------------------------------------------------------------
 
class GeoVaLs : public oops::Test {
 public:
  GeoVaLs() {}
  virtual ~GeoVaLs() {}
 private:
  std::string testid() const override {return "ufo::test::GeoVaLs";}
 
  void register_tests() const override {
    std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
 
    ts.emplace_back(CASE("ufo/GeoVaLs/testGeoVaLs")
      { testGeoVaLs(); });
  }
 
  void clear() const override {}
};
 
// ----------;
 
}  // namespace test
}  // namespace ufo
 
#endif  // TEST_UFO_GEOVALS_H_