OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
FD_op.hpp
1 /*
2  * FD_op.hpp
3  *
4  * Created on: May 1, 2020
5  * Author: Abhinav Singh
6  */
7 #ifndef FD_OP_HPP_
8 #define FD_OP_HPP_
9 #include "util/common.hpp"
10 #include "FD_expressions.hpp"
11 
12 namespace FD
13 {
14  constexpr int CENTRAL = 0;
15  constexpr int CENTRAL_ONE_SIDE_FORWARD = 1;
16  constexpr int CENTRAL_ONE_SIDE_BACKWARD = 2;
17  constexpr int CENTRAL_STAG = 3;
18  constexpr int CENTRAL_STAG_ONE_SIDE_FORWARD = 4;
19  constexpr int CENTRAL_STAG_ONE_SIDE_BACKWARD = 5;
20 
21 
22  // ord_d means order of the derivative and not order of convergence
23  template<unsigned int dir,unsigned int ord_d,unsigned int ord, unsigned int impl>
24  struct Derivative_impl
25  {
26  template<typename rtype,typename expr_type>
27  static rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key)
28  {
29  std::cerr << __FILE__ << ":" << __LINE__ << " Error the give derivative scheme with convergence order " << ord << " has not been implemented yet" << std::endl;
30  }
31  };
32 
33 
34  template<unsigned int dir>
35  struct Derivative_impl<dir,1,2,CENTRAL>
36  {
37  template<typename rtype,typename expr_type>
38  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key,comb<expr_type::gtype::dims> & c_where)
39  {
40  // x0, dx are defined in proper dir є(x, y, z)
41  auto dx = o1.getGrid().spacing(dir);
42  long int x0 = key.getKeyRef().get(dir);
43 
44  key.getKeyRef().set_d(dir, x0 + 1);
45  rtype ret = o1.value(key,c_where)*1.0/2.0;
46 
47  key.getKeyRef().set_d(dir, x0 - 1);
48  ret += o1.value(key,c_where)*(-1.0/2.0);
49 
50  key.getKeyRef().set_d(dir, x0);
51  return ret/dx;
52  }
53 
54  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
55  inline static void calculate_nz(o1_type o1,
56  const gmap_type & g_map,
57  grid_dist_key_dx<Sys_eqs::dims> & kmap,
58  const grid_sm<Sys_eqs::dims,void> & gs,
59  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
60  unordered_map_type & cols,
61  typename Sys_eqs::stype coeff,
62  unsigned int comp,
63  comb<Sys_eqs::dims> & c_where)
64  {
65  long int old_val = kmap.getKeyRef().get(dir);
66  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
67  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir]/2.0,comp,c_where);
68  kmap.getKeyRef().set_d(dir,old_val);
69 
70  old_val = kmap.getKeyRef().get(dir);
71  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
72  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir]/2.0,comp,c_where);
73  kmap.getKeyRef().set_d(dir,old_val);
74  }
75  };
76 
77  template<unsigned int dir>
78  struct Derivative_impl<dir,1,2,CENTRAL_STAG>
79  {
80  template<typename rtype,typename expr_type>
81  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key,comb<expr_type::gtype::dims> & c_where)
82  {
83  rtype ret;
84  // x0, dx are defined in proper dir є(x, y, z)
85  auto dx = o1.getGrid().spacing(dir);
86  long int x0 = key.getKeyRef().get(dir);
87 
88  if (c_where[dir] == -1)
89  {
90  char old_c = c_where[dir];
91  c_where[dir] = 0;
92  key.getKeyRef().set_d(dir, x0 - 1);
93  ret = o1.value(key,c_where)*(-1.0);
94 
95  key.getKeyRef().set_d(dir, x0);
96  ret += o1.value(key,c_where)*(1.0);
97 
98  c_where[dir] = old_c;
99  }
100  else
101  {
102  char old_c = c_where[dir];
103  c_where[dir] = -1;
104  key.getKeyRef().set_d(dir, x0);
105  ret = o1.value(key,c_where)*(-1.0);
106 
107  key.getKeyRef().set_d(dir, x0 + 1);
108  ret += o1.value(key,c_where)*(1.0);
109 
110  c_where[dir] = old_c;
111  }
112 
113  return ret/dx;
114  }
115 
116  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
117  inline static void calculate_nz(o1_type o1,
118  const gmap_type & g_map,
119  grid_dist_key_dx<Sys_eqs::dims> & kmap,
120  const grid_sm<Sys_eqs::dims,void> & gs,
121  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
122  unordered_map_type & cols,
123  typename Sys_eqs::stype coeff,
124  unsigned int comp,
125  comb<Sys_eqs::dims> & c_where)
126  {
127  // x0, dx are defined in proper dir є(x, y, z)
128  long int old_val = kmap.getKeyRef().get(dir);
129 
130  if (c_where[dir] == -1)
131  {
132  char old_c = c_where.c[dir];
133  c_where.c[dir] = 0;
134  kmap.getKeyRef().set_d(dir, old_val - 1);
135  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir],comp,c_where);
136 
137  kmap.getKeyRef().set_d(dir, old_val);
138  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir],comp,c_where);
139 
140  c_where.c[dir] = old_c;
141  }
142  else
143  {
144  char old_c = c_where.c[dir];
145  c_where.c[dir] = -1;
146  kmap.getKeyRef().set_d(dir, old_val);
147  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-coeff/spacing[dir],comp,c_where);
148 
149  kmap.getKeyRef().set_d(dir, old_val + 1);
150  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/spacing[dir],comp,c_where);
151 
152  c_where.c[dir] = old_c;
153  }
154 
155  kmap.getKeyRef().set_d(dir,old_val);
156  }
157  };
158 
159  template<unsigned int dir>
160  struct Derivative_impl<dir,2,2,CENTRAL>
161  {
162  template<typename rtype,typename expr_type>
163  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
164  {
165  // x0, dx are defined in proper dir є(x, y, z)
166  auto dx = o1.getGrid().spacing(dir);
167  long int x0 = key.getKeyRef().get(dir);
168 
169  rtype ret = o1.value(key,c_where)*(-2.0);
170 
171  key.getKeyRef().set_d(dir, x0 + 1);
172  ret += o1.value(key,c_where)*(1.0);
173 
174  key.getKeyRef().set_d(dir, x0 - 1);
175  ret += o1.value(key,c_where)*(1.0);
176 
177  key.getKeyRef().set_d(dir, x0);
178  return ret/(dx*dx);
179  }
180 
181  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
182  inline static void calculate_nz(o1_type o1,
183  const gmap_type & g_map,
184  grid_dist_key_dx<Sys_eqs::dims> & kmap,
185  const grid_sm<Sys_eqs::dims,void> & gs,
186  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
187  unordered_map_type & cols,
188  typename Sys_eqs::stype coeff,
189  unsigned int comp,
190  comb<Sys_eqs::dims> & c_where)
191  {
192  long int old_val = kmap.getKeyRef().get(dir);
193  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
194 
195  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
196  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/(spacing[dir]*spacing[dir]),comp,c_where);
197  kmap.getKeyRef().set_d(dir,old_val);
198 
199  old_val = kmap.getKeyRef().get(dir);
200  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
201  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff/(spacing[dir]*spacing[dir]),comp,c_where);
202  kmap.getKeyRef().set_d(dir,old_val);
203  }
204  };
205 
206  enum one_side_direction
207  {
208  OS_CENTRAL,
209  OS_FORWARD,
210  OS_BACKWARD,
211  };
212 
213  template<typename gtype, typename key_type>
214  one_side_direction use_one_side(gtype & g, unsigned int dir, key_type & key)
215  {
216  if (g.getDecomposition().periodicity()[dir] == true)
217  {
218  return one_side_direction::OS_CENTRAL;
219  }
220 
221  auto keyg = g.getGKey(key);
222  if (keyg.get(dir) == 0)
223  {
224  return one_side_direction::OS_FORWARD;
225  }
226  else if (keyg.get(dir) == g.getGridInfoVoid().size(dir) - 1)
227  {
228  return one_side_direction::OS_BACKWARD;
229  }
230 
231  return one_side_direction::OS_CENTRAL;
232  }
233 
234  template<unsigned int dir>
235  struct Derivative_impl<dir,1,2,CENTRAL_ONE_SIDE_FORWARD>
236  {
237  template<typename rtype,typename expr_type>
238  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
239  {
240  // x0, dx are defined in proper dir є(x, y, z)
241  auto dx = o1.getGrid().spacing(dir);
242  long int x0 = key.getKeyRef().get(dir);
243 
244  rtype ret = o1.value(key,c_where)*(-3.0/2.0);
245 
246  key.getKeyRef().set_d(dir, x0 + 1);
247  ret += o1.value(key,c_where)*2.0;
248 
249  key.getKeyRef().set_d(dir, x0 + 2);
250  ret += o1.value(key,c_where)*(-0.5);
251 
252  key.getKeyRef().set_d(dir, x0);
253  return ret/dx;
254  }
255 
256  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
257  inline static void calculate_nz(o1_type o1,
258  const gmap_type & g_map,
259  grid_dist_key_dx<Sys_eqs::dims> & kmap,
260  const grid_sm<Sys_eqs::dims,void> & gs,
261  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
262  unordered_map_type & cols,
263  typename Sys_eqs::stype coeff,
264  unsigned int comp,
265  comb<Sys_eqs::dims> & c_where)
266  { long int old_val = kmap.getKeyRef().get(dir);
267  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
268  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-3.0*coeff/spacing[dir]/2.0,comp,c_where);
269  kmap.getKeyRef().set_d(dir,old_val);
270  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
271  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[dir],comp,c_where);
272  kmap.getKeyRef().set_d(dir,old_val);
273  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 2);
274  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[dir],comp,c_where);
275  kmap.getKeyRef().set_d(dir,old_val);
276  }
277  };
278 
279  template<unsigned int dir>
280  struct Derivative_impl<dir,1,2,CENTRAL_STAG_ONE_SIDE_FORWARD>
281  {
282  template<typename rtype,typename expr_type>
283  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
284  {
285 
286  std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
287 
288  return 0.0;
289  }
290 
291  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
292  inline static void calculate_nz(o1_type o1,
293  const gmap_type & g_map,
294  grid_dist_key_dx<Sys_eqs::dims> & kmap,
295  const grid_sm<Sys_eqs::dims,void> & gs,
296  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
297  unordered_map_type & cols,
298  typename Sys_eqs::stype coeff,
299  unsigned int comp,
300  comb<Sys_eqs::dims> & c_where)
301  {
302  std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
303  }
304  };
305 
306  template<unsigned int dir>
307  struct Derivative_impl<dir,2,2,CENTRAL_ONE_SIDE_FORWARD>
308  {
309  template<typename rtype,typename expr_type>
310  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
311  {
312  // x0, dx are defined in proper dir є(x, y, z)
313  auto dx = o1.getGrid().spacing(dir);
314  long int x0 = key.getKeyRef().get(dir);
315 
316  rtype ret = o1.value(key,c_where)*(2.0);
317 
318  key.getKeyRef().set_d(dir, x0 + 1);
319  ret += o1.value(key,c_where)*(-5.0);
320 
321  key.getKeyRef().set_d(dir, x0 + 2);
322  ret += o1.value(key,c_where)*(4.0);
323 
324  key.getKeyRef().set_d(dir, x0 + 3);
325  ret += o1.value(key,c_where)*(-1.0);
326 
327  key.getKeyRef().set_d(dir, x0);
328  return ret/(dx*dx);
329  }
330 
331  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
332  inline static void calculate_nz(o1_type o1,
333  const gmap_type & g_map,
334  grid_dist_key_dx<Sys_eqs::dims> & kmap,
335  const grid_sm<Sys_eqs::dims,void> & gs,
336  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
337  unordered_map_type & cols,
338  typename Sys_eqs::stype coeff,
339  unsigned int comp,
340  comb<Sys_eqs::dims> & c_where)
341  {
342  long int old_val = kmap.getKeyRef().get(dir);
343  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
344  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
345  kmap.getKeyRef().set_d(dir,old_val);
346  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 1);
347  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-5.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
348  kmap.getKeyRef().set_d(dir,old_val);
349  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 2);
350  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,4.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
351  kmap.getKeyRef().set_d(dir,old_val);
352  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) + 3);
353  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-1.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
354  kmap.getKeyRef().set_d(dir,old_val);
355  }
356  };
357 
358 
359 
360  template<unsigned int dir>
361  struct Derivative_impl<dir,1,2,CENTRAL_ONE_SIDE_BACKWARD>
362  {
363  template<typename rtype,typename expr_type>
364  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
365  {
366  // x0, dx are defined in proper dir є(x, y, z)
367  auto dx = o1.getGrid().spacing(dir);
368  long int x0 = key.getKeyRef().get(dir);
369 
370  rtype ret = o1.value(key,c_where)*(3.0/2.0);
371 
372  key.getKeyRef().set_d(dir, x0 - 1);
373  ret += o1.value(key,c_where)*(-2.0);
374 
375  key.getKeyRef().set_d(dir, x0 - 2);
376  ret += o1.value(key,c_where)*0.5;
377 
378  key.getKeyRef().set_d(dir, x0);
379  return ret/dx;
380  }
381 
382  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
383  inline static void calculate_nz(o1_type o1,
384  const gmap_type & g_map,
385  grid_dist_key_dx<Sys_eqs::dims> & kmap,
386  const grid_sm<Sys_eqs::dims,void> & gs,
387  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
388  unordered_map_type & cols,
389  typename Sys_eqs::stype coeff,
390  unsigned int comp,
391  comb<Sys_eqs::dims> & c_where)
392  { long int old_val = kmap.getKeyRef().get(dir);
393  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
394  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,3.0*coeff/spacing[dir]/2.0,comp,c_where);
395  kmap.getKeyRef().set_d(dir,old_val);
396  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
397  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[dir],comp,c_where);
398  kmap.getKeyRef().set_d(dir,old_val);
399  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 2);
400  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[dir],comp,c_where);
401  kmap.getKeyRef().set_d(dir,old_val);
402  }
403  };
404 
405  template<unsigned int dir>
406  struct Derivative_impl<dir,1,2,CENTRAL_STAG_ONE_SIDE_BACKWARD>
407  {
408  template<typename rtype,typename expr_type>
409  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
410  {
411 
412  std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
413 
414  return 0.0;
415  }
416 
417  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
418  inline static void calculate_nz(o1_type o1,
419  const gmap_type & g_map,
420  grid_dist_key_dx<Sys_eqs::dims> & kmap,
421  const grid_sm<Sys_eqs::dims,void> & gs,
422  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
423  unordered_map_type & cols,
424  typename Sys_eqs::stype coeff,
425  unsigned int comp,
426  comb<Sys_eqs::dims> & c_where)
427  {
428  std::cout << __FILE__ << ":" << __LINE__ << " error we do not have implemented yet one sided staggered grid" << std::endl;
429  }
430  };
431 
432  template<unsigned int dir>
433  struct Derivative_impl<dir,2,2,CENTRAL_ONE_SIDE_BACKWARD>
434  {
435  template<typename rtype,typename expr_type>
436  static inline rtype calculate(expr_type & o1, grid_dist_key_dx<expr_type::gtype::dims> & key, comb<expr_type::gtype::dims> & c_where)
437  {
438  // x0, dx are defined in proper dir є(x, y, z)
439  auto dx = o1.getGrid().spacing(dir);
440  long int x0 = key.getKeyRef().get(dir);
441 
442  rtype ret = o1.value(key,c_where)*(2.0);
443 
444  key.getKeyRef().set_d(dir, x0 - 1);
445  ret += o1.value(key,c_where)*(-5.0);
446 
447  key.getKeyRef().set_d(dir, x0 - 2);
448  ret += o1.value(key,c_where)*(4.0);
449 
450  key.getKeyRef().set_d(dir, x0 - 3);
451  ret += o1.value(key,c_where)*(-1.0);
452 
453  key.getKeyRef().set_d(dir, x0);
454  return ret/(dx*dx);
455  }
456 
457  template<typename Sys_eqs, typename o1_type, typename gmap_type, typename unordered_map_type>
458  inline static void calculate_nz(o1_type o1,
459  const gmap_type & g_map,
460  grid_dist_key_dx<Sys_eqs::dims> & kmap,
461  const grid_sm<Sys_eqs::dims,void> & gs,
462  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
463  unordered_map_type & cols,
464  typename Sys_eqs::stype coeff,
465  unsigned int comp,
466  comb<Sys_eqs::dims> & c_where)
467  { long int old_val = kmap.getKeyRef().get(dir);
468  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir));
469  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,2.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
470  kmap.getKeyRef().set_d(dir,old_val);
471  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 1);
472  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-5.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
473  kmap.getKeyRef().set_d(dir,old_val);
474  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 2);
475  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,4.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
476  kmap.getKeyRef().set_d(dir,old_val);
477  kmap.getKeyRef().set_d(dir, kmap.getKeyRef().get(dir) - 3);
478  o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,-1.0*coeff/(spacing[dir]*spacing[dir]),comp,c_where);
479  kmap.getKeyRef().set_d(dir,old_val);
480  }
481  };
482 
483  template<unsigned dir_, unsigned int ord_d_, unsigned int ord_, unsigned int impl_>
484  struct GRID_DERIVATIVE
485  {
486  typedef boost::mpl::int_<dir_> dir;
487  typedef boost::mpl::int_<ord_> ord;
488  typedef boost::mpl::int_<ord_d_> ord_d;
489  typedef boost::mpl::int_<impl_> impl;
490  };
491 
492  template <typename exp1, unsigned int dir, unsigned int ord_d,unsigned int ord, unsigned int impl>
493  class grid_dist_expression_op<exp1,void,GRID_DERIVATIVE<dir,ord_d,ord,impl> >
494  {
496  const exp1 o1;
497 
498  public:
499 
500  typedef typename exp1::gtype gtype;
501 
503  inline grid_dist_expression_op(const exp1 & o1)
504  :o1(o1)
505  {}
506 
512  inline void init() const
513  {
514  o1.init();
515  }
516 
524  inline auto value(grid_dist_key_dx<gtype::dims> & key, comb<gtype::dims> & c_where) const -> typename std::remove_reference<decltype(o1.value(key,c_where))>::type
525  {
526  typedef typename std::remove_reference<decltype(o1.value(key,c_where))>::type r_type;
527 
528  r_type val;
529 
530  one_side_direction os = use_one_side(getGrid(),dir,key);
531 
532  if (os == one_side_direction::OS_CENTRAL || getGrid().is_staggered() == true)
533  {
534  val = Derivative_impl<dir,ord_d,ord,impl>::template calculate<r_type>(o1,key,c_where);
535  }
536  else if (os == one_side_direction::OS_FORWARD)
537  {
538  val = Derivative_impl<dir,ord_d,ord,impl+1>::template calculate<r_type>(o1,key,c_where);
539  }
540  else
541  {
542  val = Derivative_impl<dir,ord_d,ord,impl+2>::template calculate<r_type>(o1,key,c_where);
543  }
544 
545  return val;
546  }
547 
548 
549  template<typename Sys_eqs, typename gmap_type, typename unordered_map_type>
550  inline void value_nz(const gmap_type & g_map,
551  grid_dist_key_dx<Sys_eqs::dims> & kmap,
552  const grid_sm<Sys_eqs::dims,void> & gs,
553  typename Sys_eqs::stype (& spacing )[Sys_eqs::dims],
554  unordered_map_type & cols,
555  typename Sys_eqs::stype coeff,
556  unsigned int comp,
557  comb<Sys_eqs::dims> & c_where) const
558  {
559  //o1.template value_nz<Sys_eqs>(g_map,kmap,gs,spacing,cols,coeff,comp);
560  long int old_val = kmap.getKeyRef().get(dir);
561  one_side_direction os = use_one_side(getGrid(),dir,kmap);
562  if (os == one_side_direction::OS_CENTRAL || getGrid().is_staggered() == true)
563  {
564  Derivative_impl<dir,ord_d,ord,impl>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
565  }
566  else if (os == one_side_direction::OS_FORWARD)
567  {
568  Derivative_impl<dir,ord_d,ord,impl+1>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
569  }
570  else
571  {
572  Derivative_impl<dir,ord_d,ord,impl+2>::template calculate_nz<Sys_eqs>(o1,g_map,kmap,gs,spacing,cols,coeff,comp,c_where);
573  }
574 
575  // cols[g_map.template getProp<0>(kmap)*Sys_eqs::nvar + FD::var_id + comp] += coeff;
576  }
577 
578 
586  gtype & getGrid()
587  {
588  return o1.getGrid();
589  }
590 
598  const gtype & getGrid() const
599  {
600  return o1.getGrid();
601  }
602  };
603 
604 
605 
606 
607  template<unsigned int dir,unsigned int ord_d ,unsigned int ord, unsigned int impl>
608  class Derivative
609  {
610 
611  public:
612 
613  Derivative()
614  {
615  }
616 
617  template<typename operand_type>
618  grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<dir,ord_d,ord,impl>> operator()(operand_type arg)
619  {
620  return grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<dir,ord_d,ord,impl>>(arg);
621  }
622  };
623 
624  template<unsigned int dim, unsigned int impl>
625  class Laplacian
626  {
627 
628  public:
629 
630  Laplacian()
631  {
632  }
633 
634  template<typename operand_type>
635  FD::grid_dist_expression_op<
636  FD::grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<0,2,2,impl>>,
637  FD::grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<1,2,2,impl>>,
638  FD::sum>
639  operator()(operand_type arg)
640  {
641  grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<0,2,2,impl>> d_xx(arg);
642  grid_dist_expression_op<operand_type,void,GRID_DERIVATIVE<1,2,2,impl>> d_yy(arg);
643  return d_xx + d_yy;
644  }
645  };
646 
647  class L2Error
648  {
649  public:
650  L2Error() {}
651 
652  template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
653  auto operator()(const FD::grid_dist_expression<p1,g1,impl_p1> ga, const FD::grid_dist_expression<p2,g2,impl_p2> gb) const ->
654  typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type
655  {
656  typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type result = 0.0;
657  comb<FD::grid_dist_expression<p1,g1,impl_p1>::gtype::dims> s_pos;
658  auto diff_expr = ga + -1*gb;
659 
660  auto & v_cl = create_vcluster();
661 
662  auto it = ga.getGrid().getDomainIterator();
663  while (it.isNext())
664  {
665  auto key = it.get();
666  auto diff = diff_expr.value(key,s_pos);
667  result += diff*diff;
668  ++it;
669  }
670 
671  v_cl.sum(result);
672  v_cl.execute();
673 
674  return result;
675  }
676  };
677 
678 
679  class LInfError
680  {
681  public:
682  LInfError() {}
683 
684  template<unsigned int p1, unsigned int p2, typename g1, typename g2, unsigned int impl_p1, unsigned int impl_p2>
685  auto operator()(const FD::grid_dist_expression<p1,g1,impl_p1> ga, const FD::grid_dist_expression<p2,g2,impl_p2> gb) const ->
686  typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type
687  {
688  typename std::remove_reference<decltype(ga.getGrid().template getProp<p1>(ga.getGrid().getDomainIterator().get()))>::type result = 0.0;
689  comb<FD::grid_dist_expression<p1,g1,impl_p1>::gtype::dims> s_pos;
690  auto diff_expr = ga + -1*gb;
691 
692  auto & v_cl = create_vcluster();
693 
694  auto it = ga.getGrid().getDomainIterator();
695  while (it.isNext())
696  {
697  auto key = it.get();
698  auto diff = fabs(diff_expr.value(key,s_pos));
699 
700  if (diff > result) result = diff;
701  ++it;
702  }
703 
704  v_cl.max(result);
705  v_cl.execute();
706 
707  return result;
708  }
709  };
710 
711  typedef Derivative<0,1,2,CENTRAL_STAG> Derivative_x_stag;
712  typedef Derivative<1,1,2,CENTRAL_STAG> Derivative_y_stag;
713  typedef Derivative<0,1,2,CENTRAL> Derivative_x;
714  typedef Derivative<1,1,2,CENTRAL> Derivative_y;
715  typedef Derivative<0,2,2,CENTRAL> Derivative_xx;
716  typedef Derivative<1,2,2,CENTRAL> Derivative_yy;
717  typedef Laplacian<2,CENTRAL> Lap;
718 };
719 
720 
721 #endif
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::o1
const exp1 o1
expression 1
Definition: FD_op.hpp:496
grid_dist_key_dx::getKeyRef
base_key & getKeyRef()
Get the reference key.
Definition: grid_dist_key.hpp:91
comb
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition: comb.hpp:34
FD::grid_dist_expression_op
Definition: FD_expressions.hpp:19
FD::grid_dist_expression
Definition: FD_expressions.hpp:296
grid_dist_key_dx
Grid key for a distributed grid.
Definition: grid_dist_key.hpp:6
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::value
auto value(grid_dist_key_dx< gtype::dims > &key, comb< gtype::dims > &c_where) const -> typename std::remove_reference< decltype(o1.value(key, c_where))>::type
Evaluate the expression.
Definition: FD_op.hpp:524
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::getGrid
const gtype & getGrid() const
Return the vector on which is acting.
Definition: FD_op.hpp:598
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::init
void init() const
This function must be called before value.
Definition: FD_op.hpp:512
grid_sm< Sys_eqs::dims, void >
FD::Derivative_impl
Definition: FD_op.hpp:24
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::grid_dist_expression_op
grid_dist_expression_op(const exp1 &o1)
Costruct a FD expression out of two expressions.
Definition: FD_op.hpp:503
FD::grid_dist_expression_op< exp1, void, GRID_DERIVATIVE< dir, ord_d, ord, impl > >::getGrid
gtype & getGrid()
Return the vector on which is acting.
Definition: FD_op.hpp:586
FD::GRID_DERIVATIVE
Definition: FD_op.hpp:484
comb::c
char c[dim]
Array that store the combination.
Definition: comb.hpp:37