SABER
OoBump.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2017 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 SABER_OOPS_OOBUMP_H_
9 #define SABER_OOPS_OOBUMP_H_
10 
11 #include <algorithm>
12 #include <memory>
13 #include <sstream>
14 #include <string>
15 #include <utility>
16 #include <vector>
17 
18 #include "atlas/field.h"
19 #include "atlas/functionspace.h"
20 
21 #include "eckit/config/Configuration.h"
22 #include "eckit/mpi/Comm.h"
23 
24 #include "oops/assimilation/GMRESR.h"
25 #include "oops/base/IdentityMatrix.h"
26 #include "oops/base/Variables.h"
27 #include "oops/interface/Increment.h"
28 #include "oops/util/DateTime.h"
29 #include "oops/util/Logger.h"
30 
31 #include "saber/bump/type_bump.h"
32 
33 namespace eckit {
34  class Configuration;
35 }
36 
37 namespace oops {
38  class Variables;
39  class UnstructuredGrid;
40 }
41 
42 namespace saber {
43 
44 // -----------------------------------------------------------------------------
45 /// OoBump C++ interface
46 
47 template<typename MODEL> class OoBump {
48  typedef oops::Geometry<MODEL> Geometry_;
49  typedef oops::Increment<MODEL> Increment_;
50 
51  public:
52  OoBump(const Geometry_ &, const oops::Variables &, const util::DateTime &,
53  const eckit::LocalConfiguration);
54  explicit OoBump(OoBump &);
55  ~OoBump();
56 
57  // C++ interfaces
58  size_t getSize() {return keyOoBump_.size();}
59  int getKey(int igrid) const {return keyOoBump_[igrid];}
60  void clearKey() {keyOoBump_.clear();}
61 
62  // Fortran interfaces
63  void addMember(const atlas::FieldSet & atlasFieldSet, const int &, const int &) const;
64  void runDrivers() const;
65  void multiplyVbal(const Increment_ &, Increment_ &) const;
66  void multiplyVbalInv(const Increment_ &, Increment_ &) const;
67  void multiplyVbalAd(const Increment_ &, Increment_ &) const;
68  void multiplyVbalInvAd(const Increment_ &, Increment_ &) const;
69  void multiplyStdDev(const Increment_ &, Increment_ &) const;
70  void multiplyStdDevInv(const Increment_ &, Increment_ &) const;
71  void multiplyNicas(Increment_ &) const;
72  void multiplyNicas(const Increment_ &, Increment_ &) const;
73  void inverseMultiplyNicas(const Increment_ &, Increment_ &) const;
74  void randomize(Increment_ &) const;
75  void getParameter(const std::string &, Increment_ &) const;
76  void setParameter(const std::string &, const Increment_ &) const;
77 
78  // Aliases for inversion with GMRESR
79  void multiply(const Increment_ & dxi, Increment_ & dxo) const {multiplyNicas(dxi, dxo);}
80 
81  private:
82  std::vector<int> keyOoBump_;
83 };
84 
85 // -----------------------------------------------------------------------------
86 template<typename MODEL>
88  const oops::Variables & vars,
89  const util::DateTime & time,
90  const eckit::LocalConfiguration conf) : keyOoBump_() {
91  // Grids
92  std::vector<eckit::LocalConfiguration> grids;
93 
94  // Get global prefix
95  std::string prefix;
96  conf.get("prefix", prefix);
97 
98  // Get the grids configuration from input configuration and complete it
99  if (conf.has("grids")) {
100  // Get grids from input configuration
101  conf.get("grids", grids);
102  ASSERT(grids.size() > 0);
103  } else {
104  // Create one empty configuration
105  eckit::LocalConfiguration emptyConf;
106  grids.push_back(emptyConf);
107  }
108 
109  // Loop over grids
110  for (unsigned int jgrid = 0; jgrid < grids.size(); ++jgrid) {
111  // Add prefix
112  if (!grids[jgrid].has("prefix")) {
113  std::ostringstream ss;
114  ss << std::setw(2) << std::setfill('0') << jgrid;
115  grids[jgrid].set("prefix", prefix + "_" + ss.str());
116  }
117 
118  // Add input variables to the grid configuration
119  std::vector<std::string> vars_str;
120  if (grids[jgrid].has("variables")) {
121  grids[jgrid].get("variables", vars_str);
122  } else {
123  for (unsigned int jvar = 0; jvar < vars.size(); ++jvar) {
124  vars_str.push_back(vars[jvar]);
125  }
126  grids[jgrid].set("variables", vars_str);
127  }
128  grids[jgrid].set("nv", vars_str.size());
129 
130  // Get the required number of levels add it to the grid configuration
131  Increment_ dx(resol, vars, time);
132  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
133  dx.setAtlas(atlasFieldSet.get());
134  int nl = 0;
135  for (unsigned int jvar = 0; jvar < vars_str.size(); ++jvar) {
136  atlas::Field atlasField = atlasFieldSet->field(vars_str[jvar]);
137  nl = std::max(nl, atlasField.levels());
138  }
139  grids[jgrid].set("nl", nl);
140 
141  // Add level index for 2D fields (first or last, first by default)
142  if (!grids[jgrid].has("lev2d")) {
143  grids[jgrid].set("lev2d", "first");
144  }
145  }
146 
147  // Check grids number
148  ASSERT(grids.size() > 0);
149 
150  // Print configuration
151  oops::Log::info() << "Configuration: " << conf << std::endl;
152 
153  for (unsigned int jgrid = 0; jgrid < grids.size(); ++jgrid) {
154  // Print configuration for this grid
155  oops::Log::info() << "Grid " << jgrid << ": " << grids[jgrid] << std::endl;
156 
157  // Create OoBump instance
158  int keyOoBump = 0;
159  bump_create_f90(keyOoBump, &resol.getComm(),
160  resol.atlasFunctionSpace()->get(),
161  resol.atlasFieldSet()->get(),
162  conf, grids[jgrid]);
163  keyOoBump_.push_back(keyOoBump);
164  }
165 }
166 // -----------------------------------------------------------------------------
167 template<typename MODEL>
168 OoBump<MODEL>::OoBump(OoBump & other) : keyOoBump_() {
169  for (unsigned int jgrid = 0; jgrid < other.getSize(); ++jgrid) {
170  keyOoBump_.push_back(other.getKey(jgrid));
171  }
172  other.clearKey();
173 }
174 // -----------------------------------------------------------------------------
175 template<typename MODEL>
177  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
178  if (keyOoBump_[jgrid] > 0) bump_dealloc_f90(keyOoBump_[jgrid]);
179  }
180 }
181 // -----------------------------------------------------------------------------
182 template<typename MODEL>
183 void OoBump<MODEL>::addMember(const atlas::FieldSet & atlasFieldSet, const int & ie,
184  const int & iens) const {
185  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
186  bump_add_member_f90(keyOoBump_[jgrid], atlasFieldSet.get(), ie+1, iens);
187  }
188 }
189 // -----------------------------------------------------------------------------
190 template<typename MODEL>
192  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
193  bump_run_drivers_f90(keyOoBump_[jgrid]);
194  }
195 }
196 // -----------------------------------------------------------------------------
197 template<typename MODEL>
198 void OoBump<MODEL>::multiplyVbal(const Increment_ & dxi, Increment_ & dxo) const {
199  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
200  dxi.toAtlas(atlasFieldSet.get());
201  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
202  bump_apply_vbal_f90(keyOoBump_[jgrid], atlasFieldSet->get());
203  }
204  dxo.fromAtlas(atlasFieldSet.get());
205 }
206 // -----------------------------------------------------------------------------
207 template<typename MODEL>
208 void OoBump<MODEL>::multiplyVbalInv(const Increment_ & dxi, Increment_ & dxo) const {
209  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
210  dxi.toAtlas(atlasFieldSet.get());
211  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
212  bump_apply_vbal_inv_f90(keyOoBump_[jgrid], atlasFieldSet->get());
213  }
214  dxo.fromAtlas(atlasFieldSet.get());
215 }
216 // -----------------------------------------------------------------------------
217 template<typename MODEL>
218 void OoBump<MODEL>::multiplyVbalAd(const Increment_ & dxi, Increment_ & dxo) const {
219  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
220  dxi.toAtlas(atlasFieldSet.get());
221  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
222  bump_apply_vbal_ad_f90(keyOoBump_[jgrid], atlasFieldSet->get());
223  }
224  dxo.fromAtlas(atlasFieldSet.get());
225 }
226 // -----------------------------------------------------------------------------
227 template<typename MODEL>
229  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
230  dxi.toAtlas(atlasFieldSet.get());
231  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
232  bump_apply_vbal_inv_ad_f90(keyOoBump_[jgrid], atlasFieldSet->get());
233  }
234  dxo.fromAtlas(atlasFieldSet.get());
235 }
236 // -----------------------------------------------------------------------------
237 template<typename MODEL>
238 void OoBump<MODEL>::multiplyStdDev(const Increment_ & dxi, Increment_ & dxo) const {
239  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
240  dxi.toAtlas(atlasFieldSet.get());
241  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
242  bump_apply_stddev_f90(keyOoBump_[jgrid], atlasFieldSet->get());
243  }
244  dxo.fromAtlas(atlasFieldSet.get());
245 }
246 // -----------------------------------------------------------------------------
247 template<typename MODEL>
249  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
250  dxi.toAtlas(atlasFieldSet.get());
251  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
252  bump_apply_stddev_inv_f90(keyOoBump_[jgrid], atlasFieldSet->get());
253  }
254  dxo.fromAtlas(atlasFieldSet.get());
255 }
256 // -----------------------------------------------------------------------------
257 template<typename MODEL>
259  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
260  dx.toAtlas(atlasFieldSet.get());
261  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
262  bump_apply_nicas_f90(keyOoBump_[jgrid], atlasFieldSet->get());
263  }
264  dx.fromAtlas(atlasFieldSet.get());
265 }
266 // -----------------------------------------------------------------------------
267 template<typename MODEL>
268 void OoBump<MODEL>::multiplyNicas(const Increment_ & dxi, Increment_ & dxo) const {
269  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
270  dxi.toAtlas(atlasFieldSet.get());
271  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
272  bump_apply_nicas_f90(keyOoBump_[jgrid], atlasFieldSet->get());
273  }
274  dxo.fromAtlas(atlasFieldSet.get());
275 }
276 // -----------------------------------------------------------------------------
277 template<typename MODEL>
279  oops::IdentityMatrix<Increment_> Id;
280  dxo.zero();
281  GMRESR(dxo, dxi, *this, Id, 10, 1.0e-3);
282 }
283 // -----------------------------------------------------------------------------
284 template<typename MODEL>
286  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
287  dx.setAtlas(atlasFieldSet.get());
288  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
289  bump_randomize_f90(keyOoBump_[jgrid], atlasFieldSet->get());
290  }
291  dx.fromAtlas(atlasFieldSet.get());
292 }
293 // -----------------------------------------------------------------------------
294 template<typename MODEL>
295 void OoBump<MODEL>::getParameter(const std::string & param, Increment_ & dx) const {
296  const int nstr = param.size();
297  const char *cstr = param.c_str();
298  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
299  dx.setAtlas(atlasFieldSet.get());
300  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
301  bump_get_parameter_f90(keyOoBump_[jgrid], nstr, cstr, atlasFieldSet->get());
302  }
303  dx.fromAtlas(atlasFieldSet.get());
304 }
305 // -----------------------------------------------------------------------------
306 template<typename MODEL>
307 void OoBump<MODEL>::setParameter(const std::string & param, const Increment_ & dx) const {
308  const int nstr = param.size();
309  const char *cstr = param.c_str();
310  std::unique_ptr<atlas::FieldSet> atlasFieldSet(new atlas::FieldSet());
311  dx.toAtlas(atlasFieldSet.get());
312  for (unsigned int jgrid = 0; jgrid < keyOoBump_.size(); ++jgrid) {
313  bump_set_parameter_f90(keyOoBump_[jgrid], nstr, cstr, atlasFieldSet->get());
314  }
315 }
316 // -----------------------------------------------------------------------------
317 
318 } // namespace saber
319 
320 #endif // SABER_OOPS_OOBUMP_H_
oops
Definition: ErrorCovarianceBUMP.h:33
saber::bump_get_parameter_f90
void bump_get_parameter_f90(const int &, const int &, const char *, const atlas::field::FieldSetImpl *)
saber::OoBump::multiplyStdDevInv
void multiplyStdDevInv(const Increment_ &, Increment_ &) const
Definition: OoBump.h:248
type_bump.h
saber::OoBump::inverseMultiplyNicas
void inverseMultiplyNicas(const Increment_ &, Increment_ &) const
Definition: OoBump.h:278
saber::OoBump::randomize
void randomize(Increment_ &) const
Definition: OoBump.h:285
saber::OoBump::clearKey
void clearKey()
Definition: OoBump.h:60
saber::OoBump::runDrivers
void runDrivers() const
Definition: OoBump.h:191
saber::bump_add_member_f90
void bump_add_member_f90(const int &, const atlas::field::FieldSetImpl *, const int &, const int &)
saber::bump_dealloc_f90
void bump_dealloc_f90(const int &)
saber::bump_apply_nicas_f90
void bump_apply_nicas_f90(const int &, const atlas::field::FieldSetImpl *)
saber::OoBump::getSize
size_t getSize()
Definition: OoBump.h:58
saber::OoBump::multiplyVbal
void multiplyVbal(const Increment_ &, Increment_ &) const
Definition: OoBump.h:198
saber
Definition: type_bump.h:22
saber::OoBump
OoBump C++ interface.
Definition: OoBump.h:47
saber::bump_create_f90
void bump_create_f90(int &, const eckit::mpi::Comm *, const atlas::functionspace::FunctionSpaceImpl *, const atlas::field::FieldSetImpl *, const eckit::Configuration &, const eckit::Configuration &)
saber::bump_randomize_f90
void bump_randomize_f90(const int &, const atlas::field::FieldSetImpl *)
eckit
Definition: type_bump.h:18
saber::OoBump::multiplyVbalInv
void multiplyVbalInv(const Increment_ &, Increment_ &) const
Definition: OoBump.h:208
saber::OoBump::multiplyVbalAd
void multiplyVbalAd(const Increment_ &, Increment_ &) const
Definition: OoBump.h:218
saber::OoBump::getKey
int getKey(int igrid) const
Definition: OoBump.h:59
saber::bump_apply_vbal_ad_f90
void bump_apply_vbal_ad_f90(const int &, const atlas::field::FieldSetImpl *)
saber::OoBump::multiplyStdDev
void multiplyStdDev(const Increment_ &, Increment_ &) const
Definition: OoBump.h:238
saber::bump_apply_vbal_inv_f90
void bump_apply_vbal_inv_f90(const int &, const atlas::field::FieldSetImpl *)
saber::OoBump::setParameter
void setParameter(const std::string &, const Increment_ &) const
Definition: OoBump.h:307
saber::OoBump::Increment_
oops::Increment< MODEL > Increment_
Definition: OoBump.h:49
saber::OoBump::multiplyVbalInvAd
void multiplyVbalInvAd(const Increment_ &, Increment_ &) const
Definition: OoBump.h:228
saber::OoBump::multiply
void multiply(const Increment_ &dxi, Increment_ &dxo) const
Definition: OoBump.h:79
saber::bump_apply_stddev_f90
void bump_apply_stddev_f90(const int &, const atlas::field::FieldSetImpl *)
saber::OoBump::addMember
void addMember(const atlas::FieldSet &atlasFieldSet, const int &, const int &) const
Definition: OoBump.h:183
saber::bump_set_parameter_f90
void bump_set_parameter_f90(const int &, const int &, const char *, const atlas::field::FieldSetImpl *)
saber::OoBump::Geometry_
oops::Geometry< MODEL > Geometry_
Definition: OoBump.h:48
saber::OoBump::keyOoBump_
std::vector< int > keyOoBump_
Definition: OoBump.h:82
saber::OoBump::~OoBump
~OoBump()
Definition: OoBump.h:176
saber::bump_apply_vbal_f90
void bump_apply_vbal_f90(const int &, const atlas::field::FieldSetImpl *)
saber::bump_apply_vbal_inv_ad_f90
void bump_apply_vbal_inv_ad_f90(const int &, const atlas::field::FieldSetImpl *)
saber::bump_run_drivers_f90
void bump_run_drivers_f90(const int &)
saber::OoBump::OoBump
OoBump(const Geometry_ &, const oops::Variables &, const util::DateTime &, const eckit::LocalConfiguration)
Definition: OoBump.h:87
saber::bump_apply_stddev_inv_f90
void bump_apply_stddev_inv_f90(const int &, const atlas::field::FieldSetImpl *)
saber::OoBump::multiplyNicas
void multiplyNicas(Increment_ &) const
Definition: OoBump.h:258
saber::OoBump::getParameter
void getParameter(const std::string &, Increment_ &) const
Definition: OoBump.h:295