FV3 Bundle
oops/src/test/interface/LinearObsOperator.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2017-2018 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_LINEAROBSOPERATOR_H_
9 #define TEST_INTERFACE_LINEAROBSOPERATOR_H_
10 
11 #include <string>
12 #include <vector>
13 
14 #define BOOST_TEST_NO_MAIN
15 #define BOOST_TEST_ALTERNATIVE_INIT_API
16 #define BOOST_TEST_DYN_LINK
17 #include <boost/test/unit_test.hpp>
18 
19 #include <boost/scoped_ptr.hpp>
20 
21 #include "oops/interface/LinearObsOperator.h"
22 #include "oops/interface/ObsAuxControl.h"
23 #include "oops/interface/ObsAuxIncrement.h"
24 #include "oops/interface/ObsOperator.h"
25 #include "oops/runs/Test.h"
26 #include "oops/util/dot_product.h"
27 #include "oops/util/Logger.h"
28 #include "test/interface/ObsTestsFixture.h"
29 #include "test/TestEnvironment.h"
30 
31 namespace test {
32 
33 // -----------------------------------------------------------------------------
34 
35 template <typename MODEL> void testConstructor() {
36  typedef ObsTestsFixture<MODEL> Test_;
37  typedef oops::LinearObsOperator<MODEL> LinearObsOperator_;
38 
39  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
40  boost::scoped_ptr<LinearObsOperator_> ov(new LinearObsOperator_(Test_::obspace()[jj]));
41  BOOST_CHECK(ov.get());
42 
43  ov.reset();
44  BOOST_CHECK(!ov.get());
45  }
46 }
47 
48 // -----------------------------------------------------------------------------
49 
50 template <typename MODEL> void testLinearity() {
51  typedef ObsTestsFixture<MODEL> Test_;
52  typedef oops::GeoVaLs<MODEL> GeoVaLs_;
53  typedef oops::Locations<MODEL> Locations_;
54  typedef oops::ObsAuxControl<MODEL> ObsAuxCtrl_;
55  typedef oops::ObsAuxIncrement<MODEL> ObsAuxIncr_;
56  typedef oops::ObsOperator<MODEL> ObsOperator_;
57  typedef oops::LinearObsOperator<MODEL> LinearObsOperator_;
58  typedef oops::ObsVector<MODEL> ObsVector_;
59 
60  const double zero = 0.0;
61  const double coef = 3.14;
62  const double tol = TestEnvironment::config().getDouble("LinearObsOpTest.toleranceAD");
63  const eckit::LocalConfiguration obsconf(TestEnvironment::config(), "Observations");
64  std::vector<eckit::LocalConfiguration> conf;
65  obsconf.get("ObsTypes", conf);
66 
67  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
68  ObsOperator_ hop(Test_::obspace()[jj]);
69  LinearObsOperator_ hoptl(Test_::obspace()[jj]);
70 
71  const eckit::LocalConfiguration gconf(conf[jj], "GeoVaLs");
72  Locations_ locs(Test_::obspace()[jj].locations(Test_::tbgn(), Test_::tend()));
73  const GeoVaLs_ gval(gconf, hop.variables());
74 
75  eckit::LocalConfiguration biasConf;
76  conf[jj].get("ObsBias", biasConf);
77  const ObsAuxCtrl_ ybias(biasConf);
78  hoptl.setTrajectory(gval, ybias);
79 
80  const ObsAuxIncr_ ybinc(biasConf);
81  ObsVector_ dy1(Test_::obspace()[jj]);
82  GeoVaLs_ gv(gconf, hoptl.variables());
83 
84  gv.zero();
85  hoptl.simulateObsTL(gv, dy1, ybinc);
86 
87  BOOST_CHECK_EQUAL(dy1.rms(), zero);
88 
89  gv.random();
90  hoptl.simulateObsTL(gv, dy1, ybinc);
91  dy1 *= coef;
92  BOOST_CHECK(dy1.rms() > zero);
93 
94  gv *= coef;
95  ObsVector_ dy2(Test_::obspace()[jj]);
96  hoptl.simulateObsTL(gv, dy2, ybinc);
97 
98  dy1 -= dy2;
99 
100  BOOST_CHECK_SMALL(dy1.rms(), tol);
101  }
102 }
103 
104 // -----------------------------------------------------------------------------
105 
106 template <typename MODEL> void testAdjoint() {
107  typedef ObsTestsFixture<MODEL> Test_;
108  typedef oops::GeoVaLs<MODEL> GeoVaLs_;
109  typedef oops::Locations<MODEL> Locations_;
110  typedef oops::ObsOperator<MODEL> ObsOperator_;
111  typedef oops::LinearObsOperator<MODEL> LinearObsOperator_;
112  typedef oops::ObsAuxControl<MODEL> ObsAuxCtrl_;
113  typedef oops::ObsAuxIncrement<MODEL> ObsAuxIncr_;
114  typedef oops::ObsVector<MODEL> ObsVector_;
115 
116  const double zero = 0.0;
117  const double tol = TestEnvironment::config().getDouble("LinearObsOpTest.toleranceAD");
118  const eckit::LocalConfiguration obsconf(TestEnvironment::config(), "Observations");
119  std::vector<eckit::LocalConfiguration> conf;
120  obsconf.get("ObsTypes", conf);
121 
122  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
123  ObsOperator_ hop(Test_::obspace()[jj]);
124  LinearObsOperator_ hoptl(Test_::obspace()[jj]);
125  eckit::LocalConfiguration gconf(conf[jj], "GeoVaLs");
126  Locations_ locs(Test_::obspace()[jj].locations(Test_::tbgn(), Test_::tend()));
127  const GeoVaLs_ gval(gconf, hop.variables());
128 
129  eckit::LocalConfiguration biasConf;
130  conf[jj].get("ObsBias", biasConf);
131  const ObsAuxCtrl_ ybias(biasConf);
132 
133  hoptl.setTrajectory(gval, ybias);
134 
135  ObsAuxIncr_ ybinc(biasConf);
136 
137  ObsVector_ dy1(Test_::obspace()[jj]);
138  ObsVector_ dy2(Test_::obspace()[jj]);
139  GeoVaLs_ gv1(gconf, hoptl.variables());
140  GeoVaLs_ gv2(gconf, hoptl.variables());
141 
142  gv1.random();
143  BOOST_REQUIRE(dot_product(gv1, gv1) > zero);
144  hoptl.simulateObsTL(gv1, dy1, ybinc);
145  BOOST_CHECK(dot_product(dy1, dy1) > zero);
146 
147  dy2.random();
148  BOOST_REQUIRE(dot_product(dy2, dy2) > zero);
149  gv2.zero();
150  hoptl.simulateObsAD(gv2, dy2, ybinc);
151  BOOST_CHECK(dot_product(gv2, gv2) > zero);
152 
153  const double zz1 = dot_product(gv1, gv2);
154  const double zz2 = dot_product(dy1, dy2);
155 
156  oops::Log::debug() << "Adjoint test result: (<x,HTy>-<Hx,y>)/<Hx,y> = "
157  << (zz1-zz2)/zz2 << std::endl;
158 
159  BOOST_CHECK(zz1 != zero);
160  BOOST_CHECK(zz2 != zero);
161  BOOST_CHECK_CLOSE(zz1, zz2, tol);
162  }
163 }
164 
165 // -----------------------------------------------------------------------------
166 
167 template <typename MODEL> void testTangentLinear() {
168  // Test ||(hop(x+alpha*dx)-hop(x)) - hoptl(alpha*dx)|| < tol
169  typedef ObsTestsFixture<MODEL> Test_;
170  typedef oops::GeoVaLs<MODEL> GeoVaLs_;
171  typedef oops::Locations<MODEL> Locations_;
172  typedef oops::ObsAuxControl<MODEL> ObsAuxCtrl_;
173  typedef oops::ObsAuxIncrement<MODEL> ObsAuxIncr_;
174  typedef oops::LinearObsOperator<MODEL> LinearObsOperator_;
175  typedef oops::ObsOperator<MODEL> ObsOperator_;
176  typedef oops::ObsVector<MODEL> ObsVector_;
177 
178  const eckit::LocalConfiguration obsconf(TestEnvironment::config(), "Observations");
179  std::vector<eckit::LocalConfiguration> conf;
180  obsconf.get("ObsTypes", conf);
181 
182  const double tol = TestEnvironment::config().getDouble("LinearObsOpTest.toleranceTL");
183  const int iter = TestEnvironment::config().getDouble("LinearObsOpTest.testiterTL");
184 
185  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
186  LinearObsOperator_ hoptl(Test_::obspace()[jj]);
187  ObsOperator_ hop(Test_::obspace()[jj]);
188 
189  const eckit::LocalConfiguration gconf(conf[jj], "GeoVaLs");
190  Locations_ locs(Test_::obspace()[jj].locations(Test_::tbgn(), Test_::tend()));
191 
192  eckit::LocalConfiguration biasConf;
193  conf[jj].get("ObsBias", biasConf);
194  const ObsAuxCtrl_ ybias(biasConf);
195 
196  const ObsAuxIncr_ ybinc(biasConf);
197 
198  ObsVector_ y1(Test_::obspace()[jj]); // y1 = hop(x)
199  ObsVector_ y2(Test_::obspace()[jj]); // y2 = hop(x+alpha*dx)
200  ObsVector_ y3(Test_::obspace()[jj]); // y3 = hoptl(alpha*dx)
201 
202  GeoVaLs_ gv(gconf, hop.variables()); // Background
203 
204  hoptl.setTrajectory(gv, ybias);
205 
206  hop.simulateObs(gv, y1, ybias);
207 
208  GeoVaLs_ dgv(gconf, hoptl.variables());
209  dgv.random();
210 
211  GeoVaLs_ gv0(gconf, hop.variables());
212  gv0 = gv;
213  ObsVector_ y3_init(Test_::obspace()[jj]);
214  y3_init = y3;
215  double alpha = 0.1;
216  for (int jter = 0; jter < iter; ++jter) {
217  gv = gv0;
218  dgv *= alpha;
219  gv += dgv;
220 
221  hop.simulateObs(gv, y2, ybias);
222  y2 -= y1;
223  hoptl.simulateObsTL(dgv, y3, ybinc);
224  y2 -= y3;
225  double test_norm = y2.rms();
226  y3 = y3_init;
227  oops::Log::debug() << "Iter:" << jter << " ||(h(x+alpha*dx)-h(x))/h'(alpha*dx)||="
228  << test_norm << std::endl;
229  }
230  BOOST_CHECK(y2.rms() < tol);
231  }
232 }
233 
234 // -----------------------------------------------------------------------------
235 
236 template <typename MODEL> class LinearObsOperator : public oops::Test {
237  public:
238  LinearObsOperator() {}
239  virtual ~LinearObsOperator() {}
240  private:
241  std::string testid() const {return "test::LinearObsOperator<" + MODEL::name() + ">";}
242 
243  void register_tests() const {
244  boost::unit_test::test_suite * ts = BOOST_TEST_SUITE("interface/LinearObsOperator");
245 
246  ts->add(BOOST_TEST_CASE(&testConstructor<MODEL>));
247  ts->add(BOOST_TEST_CASE(&testLinearity<MODEL>));
248  ts->add(BOOST_TEST_CASE(&testTangentLinear<MODEL>));
249  ts->add(BOOST_TEST_CASE(&testAdjoint<MODEL>));
250  boost::unit_test::framework::master_test_suite().add(ts);
251  }
252 };
253 
254 // -----------------------------------------------------------------------------
255 
256 } // namespace test
257 
258 #endif // TEST_INTERFACE_LINEAROBSOPERATOR_H_
TestEnvironment.h
conf
Definition: conf.py:1
test::TestEnvironment::config
static const eckit::Configuration & config()
Definition: TestEnvironment.h:39
interpolator_mod::name
character(len=32) name
Definition: interpolator.F90:329
test_cases_nlm_mod::alpha
real, public alpha
Definition: test_cases_nlm.F90:98
accoeff_define::zero
real(double), parameter zero
Definition: ACCoeff_Define.f90:68
mpp_mod::debug
logical debug
Definition: mpp.F90:1297
ObsTestsFixture.h
test::testConstructor
void testConstructor()
Definition: test/interface/Geometry.h:52
cloud_interpolator_mod::tol
real, parameter tol
Definition: cloud_interpolator.F90:35
Test.h