OOPS
test/interface/LinearGetValues.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2020 UCAR
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_INTERFACE_LINEARGETVALUES_H_
9 #define TEST_INTERFACE_LINEARGETVALUES_H_
10 
11 #include <algorithm>
12 #include <cmath>
13 #include <iomanip>
14 #include <iostream>
15 #include <limits>
16 #include <memory>
17 #include <string>
18 #include <vector>
19 
20 #define ECKIT_TESTING_SELF_REGISTER_CASES 0
21 
22 #include <boost/noncopyable.hpp>
23 #include <boost/ptr_container/ptr_vector.hpp>
24 
25 #include "eckit/config/LocalConfiguration.h"
26 #include "eckit/testing/Test.h"
27 #include "oops/base/Variables.h"
28 #include "oops/interface/Geometry.h"
29 #include "oops/interface/GeoVaLs.h"
30 #include "oops/interface/GetValues.h"
31 #include "oops/interface/Increment.h"
32 #include "oops/interface/LinearGetValues.h"
33 #include "oops/interface/Locations.h"
34 #include "oops/interface/State.h"
35 #include "oops/mpi/mpi.h"
36 #include "oops/runs/Test.h"
37 #include "oops/util/DateTime.h"
38 #include "oops/util/dot_product.h"
39 #include "oops/util/Duration.h"
40 #include "test/TestEnvironment.h"
41 
42 namespace test {
43 
44 // =================================================================================================
45 
46 template <typename MODEL, typename OBS> class LinearGetValuesFixture : private boost::noncopyable {
47  typedef eckit::LocalConfiguration LocalConfig_;
48  typedef oops::GeoVaLs<OBS> GeoVaLs_;
49  typedef oops::Geometry<MODEL> Geometry_;
50  typedef oops::GetValues<MODEL, OBS> GetValues_;
51  typedef oops::LinearGetValues<MODEL, OBS> LinearGetValues_;
52  typedef oops::Locations<OBS> Locations_;
53  typedef oops::State<MODEL> State_;
54  typedef oops::Variables Variables_;
55  typedef util::DateTime DateTime_;
56 
57  public:
58  static const DateTime_ & time() {return *getInstance().time_;}
59  static const DateTime_ & timebeg() {return *getInstance().timebeg_;}
60  static const DateTime_ & timeend() {return *getInstance().timeend_;}
61  static const GeoVaLs_ & geovals() {return *getInstance().geovals_;}
62  static const Geometry_ & resol() {return *getInstance().resol_;}
63  static const GetValues_ & getvalues() {return *getInstance().getvalues_;}
64  static const LinearGetValues_ & lineargetvalues() {return *getInstance().lineargetvalues_;}
65  static const LocalConfig_ & testconf() {return *getInstance().testconf_;}
66  static const Locations_ & locs() {return *getInstance().locs_;}
67  static const State_ & state() {return *getInstance().state_;}
68  static const Variables_ & statevars() {return *getInstance().statevars_;}
69  static const Variables_ & geovalvars() {return *getInstance().geovalvars_;}
70 
71  private:
72  static LinearGetValuesFixture<MODEL, OBS>& getInstance() {
73  static LinearGetValuesFixture<MODEL, OBS> theLinearGetValuesFixture;
74  return theLinearGetValuesFixture;
75  }
76 
77  LinearGetValuesFixture<MODEL, OBS>() {
78  testconf_.reset(new LocalConfig_(TestEnvironment::config(), "linear getvalues test"));
79 
80  // Geometry
81  const LocalConfig_ resolConfig(TestEnvironment::config(), "geometry");
82  resol_.reset(new Geometry_(resolConfig, oops::mpi::world()));
83 
84  // Variables
85  geovalvars_.reset(new Variables_(TestEnvironment::config(), "state variables"));
86 
87  // Locations
88  const LocalConfig_ locsConfig(TestEnvironment::config(), "locations");
89  locs_.reset(new Locations_(locsConfig, oops::mpi::world()));
90 
91  // Window times
92  timebeg_.reset(new DateTime_(locsConfig.getString("window begin")));
93  timeend_.reset(new DateTime_(locsConfig.getString("window end")));
94 
95  // GeoVaLs
96  geovals_.reset(new GeoVaLs_(*locs_, *geovalvars_));
97 
98  // Nonlinear GetValues
99  getvalues_.reset(new GetValues_(*resol_, *locs_));
100 
101  // State
102  const LocalConfig_ stateConfig(TestEnvironment::config(), "background");
103  state_.reset(new State_(*resol_, stateConfig));
104  statevars_.reset(new Variables_(state_->variables()));
105 
106  // Valid time
107  time_.reset(new DateTime_(state_->validTime()));
108 
109  // LinearGetValues
110  lineargetvalues_.reset(new LinearGetValues_(*resol_, *locs_));
111 
112  // Set trajectory
113  GeoVaLs_ gvtraj(*locs_, *geovalvars_);
114  lineargetvalues_->setTrajectory(*state_, *timebeg_, *timeend_, gvtraj);
115  }
116 
117  ~LinearGetValuesFixture<MODEL, OBS>() {}
118 
119  std::unique_ptr<const DateTime_> time_;
120  std::unique_ptr<const DateTime_> timebeg_;
121  std::unique_ptr<const DateTime_> timeend_;
122  std::unique_ptr<const GeoVaLs_> geovals_;
123  std::unique_ptr<const Geometry_> resol_;
124  std::unique_ptr<const GetValues_> getvalues_;
125  std::unique_ptr<LinearGetValues_> lineargetvalues_;
126  std::unique_ptr<const LocalConfig_> testconf_;
127  std::unique_ptr<const Locations_> locs_;
128  std::unique_ptr<const State_> state_;
129  std::unique_ptr<const Variables_> statevars_;
130  std::unique_ptr<const Variables_> geovalvars_;
131 };
132 
133 // =================================================================================================
134 
135 template <typename MODEL, typename OBS> void testLinearGetValuesConstructor() {
136  typedef LinearGetValuesFixture<MODEL, OBS> Test_;
137  typedef oops::LinearGetValues<MODEL, OBS> LinearGetValues_;
138 
139  std::unique_ptr<const LinearGetValues_> lineargetvalues(new LinearGetValues_(Test_::resol(),
140  Test_::locs()));
141  EXPECT(lineargetvalues.get());
142 
143  lineargetvalues.reset();
144  EXPECT(!lineargetvalues.get());
145 }
146 
147 // -------------------------------------------------------------------------------------------------
148 
149 template <typename MODEL, typename OBS> void testLinearGetValuesZeroPert() {
150  typedef LinearGetValuesFixture<MODEL, OBS> Test_;
151  typedef oops::GeoVaLs<OBS> GeoVaLs_;
152  typedef oops::Increment<MODEL> Increment_;
153 
154  Increment_ dx(Test_::resol(), Test_::statevars(), Test_::time());
155  dx.zero();
156 
157  GeoVaLs_ gv(Test_::locs(), Test_::geovalvars());
158 
159  EXPECT(dx.norm() == 0.0);
160 
161  Test_::lineargetvalues().fillGeoVaLsTL(dx, Test_::timebeg(), Test_::timeend(), gv);
162 
163  EXPECT(dx.norm() == 0.0);
164  EXPECT(gv.rms() == 0.0);
165 
166  Test_::lineargetvalues().fillGeoVaLsAD(dx, Test_::timebeg(), Test_::timeend(), gv);
167 
168  EXPECT(dx.norm() == 0.0);
169  EXPECT(gv.rms() == 0.0);
170 }
171 
172 // -------------------------------------------------------------------------------------------------
173 
174 template <typename MODEL, typename OBS> void testLinearGetValuesLinearity() {
175  typedef LinearGetValuesFixture<MODEL, OBS> Test_;
176  typedef oops::GeoVaLs<OBS> GeoVaLs_;
177  typedef oops::Increment<MODEL> Increment_;
178 
179  const double zz = 3.1415;
180 
181  Increment_ dx1(Test_::resol(), Test_::statevars(), Test_::time());
182  dx1.random();
183  Increment_ dx2(dx1);
184 
185  EXPECT(dx1.norm() > 0.0);
186  EXPECT(dx2.norm() > 0.0);
187 
188  GeoVaLs_ gv1(Test_::locs(), Test_::geovalvars());
189  GeoVaLs_ gv2(Test_::locs(), Test_::geovalvars());
190 
191  // Compute geovals
192  Test_::lineargetvalues().fillGeoVaLsTL(dx1, Test_::timebeg(), Test_::timeend(), gv1);
193 
194  gv1 *= zz;
195  dx2 *= zz;
196 
197  // Compute geovals
198  Test_::lineargetvalues().fillGeoVaLsTL(dx2, Test_::timebeg(), Test_::timeend(), gv2);
199 
200  const double tol = Test_::testconf().getDouble("tolerance linearity", 1.0e-11);
201  EXPECT(oops::is_close(gv1.rms(), gv2.rms(), tol));
202 }
203 
204 // -------------------------------------------------------------------------------------------------
205 
206 template <typename MODEL, typename OBS> void testLinearGetValuesLinearApproximation() {
207  typedef LinearGetValuesFixture<MODEL, OBS> Test_;
208  typedef oops::GeoVaLs<OBS> GeoVaLs_;
209  typedef oops::Increment<MODEL> Increment_;
210  typedef oops::State<MODEL> State_;
211 
212  const unsigned int ntest = Test_::testconf().getInt("iterations TL", 10);
213  double zz = Test_::testconf().getDouble("first multiplier TL", 1.0e-2);
214 
215  // Compute nonlinear geovals
216  State_ xx0(Test_::state());
217  GeoVaLs_ gv0(Test_::locs(), Test_::geovalvars());
218  Test_::getvalues().fillGeoVaLs(xx0, Test_::timebeg(), Test_::timeend(), gv0);
219 
220  // Run tangent linear
221  Increment_ dx0(Test_::resol(), Test_::statevars(), Test_::time());
222  dx0.random();
223 
224  std::vector<double> errors;
225  for (unsigned int jtest = 0; jtest < ntest; ++jtest) {
226  Increment_ dxx(dx0);
227  State_ xx(xx0);
228  GeoVaLs_ gv(Test_::locs(), Test_::geovalvars());
229  GeoVaLs_ dgv(Test_::locs(), Test_::geovalvars());
230  dxx *= zz;
231  xx += dxx;
232 
233  // Nonlinear
234  Test_::getvalues().fillGeoVaLs(xx, Test_::timebeg(), Test_::timeend(), gv);
235 
236  // Tangent linear
237  Test_::lineargetvalues().fillGeoVaLsTL(dxx, Test_::timebeg(), Test_::timeend(), dgv);
238 
239  GeoVaLs_ gvdiff(gv);
240  gvdiff -= gv0;
241 
242  // Print the norms
243  const double nlnorm = gvdiff.rms();
244  const double tlnorm = dgv.rms();
245 
246  // Subtract the tlm pert
247  gvdiff -= dgv;
248 
249  double testdot = dot_product(gvdiff, gvdiff);
250  errors.push_back(testdot);
251 
252  oops::Log::test() << " ||g(x+dx) - g(x)|| = " << std::setprecision(16) << nlnorm << std::endl;
253  oops::Log::test() << " ||Gdx|| = " << std::setprecision(16) << tlnorm << std::endl;
254  oops::Log::test() << "||g(x+dx)-g(x)-Gdx|| = " << std::setprecision(16) << testdot << std::endl;
255 
256  zz /= 10.0;
257  }
258 
259  // Analyze results
260  const double approx = *std::min_element(errors.begin(), errors.end());
261  oops::Log::test() << "Test LinearGetValuesTL min error = " << approx << std::endl;
262  const double tol = Test_::testconf().getDouble("tolerance TL", 1.0e-11);
263  EXPECT(approx < tol);
264 }
265 
266 // -------------------------------------------------------------------------------------------------
267 
268 template <typename MODEL, typename OBS> void testLinearGetValuesAdjoint() {
269  typedef LinearGetValuesFixture<MODEL, OBS> Test_;
270  typedef oops::GeoVaLs<OBS> GeoVaLs_;
271  typedef oops::Increment<MODEL> Increment_;
272 
273  Increment_ dx_in(Test_::resol(), Test_::statevars(), Test_::time());
274  Increment_ dx_ou(Test_::resol(), Test_::statevars(), Test_::time());
275 
276  GeoVaLs_ gv_ou(Test_::locs(), Test_::geovalvars());
277 
278  // Tangent linear
279  dx_in.random();
280  EXPECT(dx_in.norm() > 0.0);
281  Test_::lineargetvalues().fillGeoVaLsTL(dx_in, Test_::timebeg(), Test_::timeend(), gv_ou);
282  EXPECT(gv_ou.rms() > 0.0);
283 
284  // Adjoint
285  GeoVaLs_ gv_in(gv_ou);
286  gv_in.random(); // No order dependency but need a copy
287  EXPECT(gv_in.rms() > 0.0);
288  dx_ou.zero();
289  Test_::lineargetvalues().fillGeoVaLsAD(dx_ou, Test_::timebeg(), Test_::timeend(), gv_in);
290  EXPECT(dx_ou.norm() > 0.0);
291 
292  // Dot products
293  const double dot1 = dot_product(dx_in, dx_ou);
294  const double dot2 = dot_product(gv_in, gv_ou);
295  const double tol = Test_::testconf().getDouble("tolerance AD", 1.0e-11);
296  EXPECT(oops::is_close(dot1, dot2, tol));
297 
298  oops::Log::test() << "Dot Product <dx, M^Tgv> = " << std::setprecision(16) << dot1 << std::endl;
299  oops::Log::test() << "Dot Product <gv, M dx> = " << std::setprecision(16) << dot2 << std::endl;
300  oops::Log::test() << "Relative diff: " << std::setprecision(16) << (dot1-dot2)/dot1 << std::endl;
301 }
302 
303 // =================================================================================================
304 
305 template <typename MODEL, typename OBS>
306 class LinearGetValues : public oops::Test {
307  public:
308  LinearGetValues() {}
309  virtual ~LinearGetValues() {}
310 
311  private:
312  std::string testid() const override {return "test::LinearGetValues<" + MODEL::name() + ", "
313  + OBS::name() + ">";}
314 
315  void register_tests() const override {
316  std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
317 
318  ts.emplace_back(CASE("interface/GeometryIterator/testLinearGetValuesConstructor")
319  { testLinearGetValuesConstructor<MODEL, OBS>(); });
320  ts.emplace_back(CASE("interface/GeometryIterator/testLinearGetValuesZeroPert")
321  { testLinearGetValuesZeroPert<MODEL, OBS>(); });
322  ts.emplace_back(CASE("interface/GeometryIterator/testLinearGetValuesLinearity")
323  { testLinearGetValuesLinearity<MODEL, OBS>(); });
324  ts.emplace_back(CASE("interface/GeometryIterator/testLinearGetValuesLinearApproximation")
325  { testLinearGetValuesLinearApproximation<MODEL, OBS>(); });
326  ts.emplace_back(CASE("interface/GeometryIterator/testLinearGetValuesAdjoint")
327  { testLinearGetValuesAdjoint<MODEL, OBS>(); });
328  }
329 
330  void clear() const override {}
331 };
332 
333 // =================================================================================================
334 
335 } // namespace test
336 
337 #endif // TEST_INTERFACE_LINEARGETVALUES_H_
test::LinearGetValuesFixture::timeend
static const DateTime_ & timeend()
Definition: test/interface/LinearGetValues.h:60
test::testLinearGetValuesConstructor
void testLinearGetValuesConstructor()
Definition: test/interface/LinearGetValues.h:135
test::LinearGetValuesFixture::statevars_
std::unique_ptr< const Variables_ > statevars_
Definition: test/interface/LinearGetValues.h:129
test::LinearGetValues::register_tests
void register_tests() const override
Definition: test/interface/LinearGetValues.h:315
test::LinearGetValuesFixture::time
static const DateTime_ & time()
Definition: test/interface/LinearGetValues.h:58
test::LinearGetValuesFixture::getInstance
static LinearGetValuesFixture< MODEL, OBS > & getInstance()
Definition: test/interface/LinearGetValues.h:72
test::LinearGetValuesFixture::resol
static const Geometry_ & resol()
Definition: test/interface/LinearGetValues.h:62
GetValues.h
Locations.h
test::LinearGetValuesFixture::lineargetvalues
static const LinearGetValues_ & lineargetvalues()
Definition: test/interface/LinearGetValues.h:64
test::LinearGetValuesFixture::testconf_
std::unique_ptr< const LocalConfig_ > testconf_
Definition: test/interface/LinearGetValues.h:126
mpi.h
test::LinearGetValuesFixture::getvalues
static const GetValues_ & getvalues()
Definition: test/interface/LinearGetValues.h:63
test::LinearGetValuesFixture::state_
std::unique_ptr< const State_ > state_
Definition: test/interface/LinearGetValues.h:128
test::LinearGetValuesFixture::State_
oops::State< MODEL > State_
Definition: test/interface/LinearGetValues.h:53
test::LinearGetValuesFixture::statevars
static const Variables_ & statevars()
Definition: test/interface/LinearGetValues.h:68
test::LinearGetValuesFixture
Definition: test/interface/LinearGetValues.h:46
test::LinearGetValues::clear
void clear() const override
Definition: test/interface/LinearGetValues.h:330
test::LinearGetValuesFixture::timebeg_
std::unique_ptr< const DateTime_ > timebeg_
Definition: test/interface/LinearGetValues.h:120
test::CASE
CASE("test_linearmodelparameterswrapper_valid_name")
Definition: LinearModelFactory.cc:22
test::LinearGetValuesFixture::GetValues_
oops::GetValues< MODEL, OBS > GetValues_
Definition: test/interface/LinearGetValues.h:50
test::LinearGetValues
Definition: test/interface/LinearGetValues.h:306
test::LinearGetValuesFixture::time_
std::unique_ptr< const DateTime_ > time_
Definition: test/interface/LinearGetValues.h:119
oops::GetValues
Gets values from model State to observation locations (fills GeoVaLs)
Definition: oops/interface/GetValues.h:32
test::LinearGetValuesFixture::resol_
std::unique_ptr< const Geometry_ > resol_
Definition: test/interface/LinearGetValues.h:123
test::LinearGetValues::~LinearGetValues
virtual ~LinearGetValues()
Definition: test/interface/LinearGetValues.h:309
test::LinearGetValuesFixture::LinearGetValues_
oops::LinearGetValues< MODEL, OBS > LinearGetValues_
Definition: test/interface/LinearGetValues.h:51
test
Definition: LinearModelFactory.cc:20
test::LinearGetValuesFixture::Locations_
oops::Locations< OBS > Locations_
Definition: test/interface/LinearGetValues.h:52
test::LinearGetValuesFixture::LocalConfig_
eckit::LocalConfiguration LocalConfig_
Definition: test/interface/LinearGetValues.h:47
test::LinearGetValuesFixture::Geometry_
oops::Geometry< MODEL > Geometry_
Definition: test/interface/LinearGetValues.h:49
test::LinearGetValuesFixture::locs_
std::unique_ptr< const Locations_ > locs_
Definition: test/interface/LinearGetValues.h:127
test::LinearGetValuesFixture::geovals_
std::unique_ptr< const GeoVaLs_ > geovals_
Definition: test/interface/LinearGetValues.h:122
test::LinearGetValuesFixture::state
static const State_ & state()
Definition: test/interface/LinearGetValues.h:67
test::LinearGetValuesFixture::geovalvars
static const Variables_ & geovalvars()
Definition: test/interface/LinearGetValues.h:69
oops::LinearGetValues
sets trajectory and computes TL and AD for GetValues
Definition: oops/interface/LinearGetValues.h:33
test::LinearGetValuesFixture::geovals
static const GeoVaLs_ & geovals()
Definition: test/interface/LinearGetValues.h:61
LinearGetValues.h
test::testLinearGetValuesLinearity
void testLinearGetValuesLinearity()
Definition: test/interface/LinearGetValues.h:174
test::LinearGetValuesFixture::timebeg
static const DateTime_ & timebeg()
Definition: test/interface/LinearGetValues.h:59
Test.h
test::TestEnvironment::config
static const eckit::Configuration & config()
Definition: TestEnvironment.h:40
test::LinearGetValuesFixture::testconf
static const LocalConfig_ & testconf()
Definition: test/interface/LinearGetValues.h:65
test::LinearGetValuesFixture::timeend_
std::unique_ptr< const DateTime_ > timeend_
Definition: test/interface/LinearGetValues.h:121
test::testLinearGetValuesLinearApproximation
void testLinearGetValuesLinearApproximation()
Definition: test/interface/LinearGetValues.h:206
test::testLinearGetValuesZeroPert
void testLinearGetValuesZeroPert()
Definition: test/interface/LinearGetValues.h:149
test::LinearGetValuesFixture::DateTime_
util::DateTime DateTime_
Definition: test/interface/LinearGetValues.h:55
test::LinearGetValues::LinearGetValues
LinearGetValues()
Definition: test/interface/LinearGetValues.h:308
test::LinearGetValuesFixture::Variables_
oops::Variables Variables_
Definition: test/interface/LinearGetValues.h:54
test::LinearGetValuesFixture::geovalvars_
std::unique_ptr< const Variables_ > geovalvars_
Definition: test/interface/LinearGetValues.h:130
test::LinearGetValues::testid
std::string testid() const override
Definition: test/interface/LinearGetValues.h:312
oops::Geometry
Geometry class used in oops; subclass of interface class above.
Definition: oops/interface/Geometry.h:189
test::LinearGetValuesFixture::locs
static const Locations_ & locs()
Definition: test/interface/LinearGetValues.h:66
test::LinearGetValuesFixture::lineargetvalues_
std::unique_ptr< LinearGetValues_ > lineargetvalues_
Definition: test/interface/LinearGetValues.h:125
test::LinearGetValuesFixture::getvalues_
std::unique_ptr< const GetValues_ > getvalues_
Definition: test/interface/LinearGetValues.h:124
TestEnvironment.h
test::testLinearGetValuesAdjoint
void testLinearGetValuesAdjoint()
Definition: test/interface/LinearGetValues.h:268
oops::State
Encapsulates the model state.
Definition: CostJbState.h:28
oops::Test
Definition: Test.h:39
State.h
oops::mpi::world
const eckit::mpi::Comm & world()
Default communicator with all MPI tasks (ie MPI_COMM_WORLD)
Definition: oops/mpi/mpi.cc:22
oops::Variables
Definition: oops/base/Variables.h:23
oops::Increment
Increment Class: Difference between two states.
Definition: CostJbState.h:27
oops::Locations
Locations of observations for observation operator.
Definition: oops/interface/Locations.h:36
GeoVaLs.h
oops::GeoVaLs
Definition: oops/interface/GeoVaLs.h:32
test::LinearGetValuesFixture::GeoVaLs_
oops::GeoVaLs< OBS > GeoVaLs_
Definition: test/interface/LinearGetValues.h:48
Variables.h
Geometry.h
Increment.h