UFO
test/ufo/GeoVaLs.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2019 UK Met Office
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
8 #ifndef TEST_UFO_GEOVALS_H_
9 #define TEST_UFO_GEOVALS_H_
10 
11 #include <cmath>
12 #include <iomanip>
13 #include <iostream>
14 #include <memory>
15 #include <string>
16 #include <vector>
17 
18 #define ECKIT_TESTING_SELF_REGISTER_CASES 0
19 
20 #include "eckit/config/LocalConfiguration.h"
21 #include "eckit/testing/Test.h"
22 #include "ioda/ObsSpace.h"
23 #include "oops/mpi/mpi.h"
24 #include "oops/runs/Test.h"
25 #include "oops/util/FloatCompare.h"
26 #include "oops/util/Logger.h"
27 #include "test/TestEnvironment.h"
28 #include "ufo/GeoVaLs.h"
29 #include "ufo/Locations.h"
30 #include "ufo/ObsOperator.h"
31 
32 namespace ufo {
33 namespace test {
34 
35 // -----------------------------------------------------------------------------
36 
37 void testGeoVaLs() {
38  const eckit::LocalConfiguration conf(::test::TestEnvironment::config());
39  util::DateTime bgn(conf.getString("window begin"));
40  util::DateTime end(conf.getString("window end"));
41 
42  std::vector<eckit::LocalConfiguration> confs;
43  conf.get("geovals test", confs);
44  for (size_t jconf = 0; jconf < confs.size(); ++jconf) {
45 /// Setup ObsSpace
46  const eckit::LocalConfiguration obsconf(confs[jconf], "obs space");
47  ioda::ObsSpace ospace(obsconf, oops::mpi::world(), bgn, end, oops::mpi::myself());
48 
49 /// Setup GeoVaLs
50  const eckit::LocalConfiguration gconf(confs[jconf], "geovals");
51  const oops::Variables ingeovars(gconf, "state variables");
52  const GeoVaLs gval(gconf, ospace, ingeovars);
53 
54  const double tol = gconf.getDouble("tolerance");
55 
56 /// Check that GeoVaLs default constructor works
57  oops::Log::trace() <<
58  "GeoVaLs default constructor - does not allocate fields" << std::endl;
59  GeoVaLs gv_temp(ospace.comm());
60 
61 /// Check that GeoVaLs constructor to create a GeoVaLs with one location works
62  if (gconf.has("one location check")) {
63  oops::Log::trace() << "Check that GeoVaLs constructor for one location works" << std::endl;
64  const eckit::LocalConfiguration gconfone(gconf, "one location check");
65  const std::string var = gconfone.getString("variable");
66  const std::vector<int> ind = gconfone.getIntVector("indices");
67  const std::vector<float> values = gconfone.getFloatVector("values");
68  const float oneloctol = gconfone.getFloat("tolerance");
69 
70  // Loop over each location and test just the lowest level
71  oops::TestVerbosity verbosity = oops::TestVerbosity::LOG_SUCCESS_AND_FAILURE;
72  for (std::size_t i = 0; i < ind.size(); ++i) {
73  GeoVaLs gv_one(gval, ind[i]);
74  std::vector<float> gv_val(1);
75  gv_one.get(gv_val, var, 1);
76  EXPECT(oops::is_close_absolute(gv_val[0], values[i], oneloctol, verbosity));
77  }
78  } else {
79  oops::Log::trace() << "Test just the constructor for a one location GeoVaLs" << std::endl;
80  int index = 0;
81  GeoVaLs gv_one(gval, index);
82  }
83 
84  GeoVaLs gv(gval);
85  if (gconf.has("reorderzdir check")) {
86  const eckit::LocalConfiguration gconfchk(gconf, "reorderzdir check");
87  const std::string flipto = gconfchk.getString("direction");
88  std::string flipback = (flipto == "bottom2top") ? "top2bottom" : "bottom2top";
89  std::size_t nobs = ospace.nlocs();
90  gv.reorderzdir("air_pressure_levels", flipto);
91  std::vector<float> gvar(nobs);
92  std::vector<float> gvarref(nobs);
93  float sum;
94  for (size_t i = 0; i < ingeovars.size(); ++i) {
95  size_t nlevs = gv.nlevs(ingeovars[i]);
96  sum = 0;
97  for (size_t k = 0; k < nlevs; ++k) {
98  size_t kk = nlevs - k;
99  gv.get(gvar, ingeovars[i], k+1);
100  gval.get(gvarref, ingeovars[i], kk);
101  for (size_t j = 0; j < nobs; ++j) {
102  gvar[j] = gvar[j] - gvarref[j];
103  sum += sum + gvar[j];
104  }
105  }
106  }
107  gv.reorderzdir("air_pressure_levels", flipback);
108  const double tol = gconfchk.getDouble("tolerance");
109  EXPECT(abs(sum) < tol);
110  }
111 
112 /// Check that GeoVaLs merge followed by a split gives back the original geovals
113  oops::Log::trace() <<
114  "GeoVaLs merge followed by a split gives back the original geovals" << std::endl;
115 
116  double dp_gval = gval.dot_product_with(gval);
117  oops::Log::debug()<< "initial gval dot product with itself " << dp_gval << std::endl;
118 
119  gv.zero();
120  gv.merge(gval, gval);
121 
122  double dp_gv_merged = gv.dot_product_with(gv);
123  oops::Log::debug()<< "gv dot product with itself after merge " << dp_gv_merged << std::endl;
124  EXPECT(abs(dp_gv_merged - 2.0 * dp_gval)/dp_gv_merged < tol);
125 
126  GeoVaLs gv1(ospace.comm()), gv2(ospace.comm());
127  gv.split(gv1, gv2);
128 
129  double dp_gv1_split = gv1.dot_product_with(gv1);
130  double dp_gv2_split = gv2.dot_product_with(gv2);
131 
132  oops::Log::debug()<< "gv1 gv2 dot products with itself after split "
133  << dp_gv1_split << " " << dp_gv2_split << std::endl;
134 
135  EXPECT(gv1.rms() == gv2.rms());
136  EXPECT(gv1.rms() == gval.rms());
137  EXPECT(abs(dp_gv1_split - dp_gv2_split)/dp_gv1_split < tol);
138  EXPECT(abs(dp_gv1_split - dp_gval)/dp_gv1_split < tol);
139 
140  oops::Log::trace() <<
141  "GeoVaLs merge followed by a split test succeeded" << std::endl;
142 
143 /// Check that GeoVaLs & operator *= (const std::vector<float>);
144  oops::Log::trace() <<
145  "Check that GeoVaLs & operator *= (const std::vector<float>);" << std::endl;
146 
147  std::size_t nlocs = ospace.nlocs();
148  {
149  std::vector<float> tw(nlocs, 1.0f);
150  gv1 *= tw;
151  EXPECT(gv1.rms() == gval.rms());
152  }
153  {
154  std::vector<float> tw(nlocs, 2.0f);
155  gv1 *= tw;
156  double rms1, rms2;
157  rms1 = gv1.rms();
158  rms2 = 2.0 * gval.rms();
159  oops::Log::debug()<< "rms1, rms2 = " << rms1 << " " << rms2 << std::endl;
160  EXPECT(rms1 == rms2);
161  }
162  oops::Log::trace() <<
163  "GeoVaLs & operator *= (const std::vector<float>); test succeeded" << std::endl;
164  }
165 }
166 // -----------------------------------------------------------------------------
167 
168 class GeoVaLs : public oops::Test {
169  public:
170  GeoVaLs() {}
171  virtual ~GeoVaLs() {}
172  private:
173  std::string testid() const override {return "ufo::test::GeoVaLs";}
174 
175  void register_tests() const override {
176  std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
177 
178  ts.emplace_back(CASE("ufo/GeoVaLs/testGeoVaLs")
179  { testGeoVaLs(); });
180  }
181 
182  void clear() const override {}
183 };
184 
185 // ----------;
186 
187 } // namespace test
188 } // namespace ufo
189 
190 #endif // TEST_UFO_GEOVALS_H_
ufo::test::GeoVaLs
Definition: test/ufo/GeoVaLs.h:168
ufo::test::CASE
CASE("ufo/MetOfficeBuddyPairFinder/" "Duplicates, constraints on buddy counts, legacy pair collector")
Definition: test/ufo/MetOfficeBuddyPairFinder.h:171
ufo_radiancerttov_utils_mod::debug
logical, public debug
Definition: ufo_radiancerttov_utils_mod.F90:100
ObsOperator.h
ufo::test::GeoVaLs::~GeoVaLs
virtual ~GeoVaLs()
Definition: test/ufo/GeoVaLs.h:171
ufo::test::GeoVaLs::testid
std::string testid() const override
Definition: test/ufo/GeoVaLs.h:173
ufo
Definition: RunCRTM.h:27
ufo::abs
util::Duration abs(const util::Duration &duration)
Definition: TrackCheckUtils.h:31
ufo::test::GeoVaLs::register_tests
void register_tests() const override
Definition: test/ufo/GeoVaLs.h:175
ufo::test::GeoVaLs::GeoVaLs
GeoVaLs()
Definition: test/ufo/GeoVaLs.h:170
ufo::test::GeoVaLs::clear
void clear() const override
Definition: test/ufo/GeoVaLs.h:182
ufo::test::testGeoVaLs
void testGeoVaLs()
Definition: test/ufo/GeoVaLs.h:37
conf
Definition: conf.py:1