1/*
2 * FiniteDifferences.hpp
3 *
4 * Created on: Sep 17, 2015
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
9#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
10
11#include "../Matrix/SparseMatrix.hpp"
12#include "Grid/grid_dist_id.hpp"
13#include "Grid/Iterators/grid_dist_id_iterator_sub.hpp"
14#include "eq.hpp"
15#include "NN/CellList/CellDecomposer.hpp"
16#include "Grid/staggered_dist_grid_util.hpp"
17#include "Grid/grid_dist_id.hpp"
18#include "Vector/Vector_util.hpp"
19#include "Grid/staggered_dist_grid.hpp"
20
21/*! \brief Finite Differences
22 *
23 * This class is able to discretize on a Matrix any system of equations producing a linear system of type \f$Ax=b\f$. In order to create a consistent
24 * Matrix it is required that each processor must contain a contiguous range on grid points without
25 * holes. In order to ensure this, each processor produce a contiguous local labeling of its local
26 * points. Each processor also add an offset equal to the number of local
27 * points of the processors with id smaller than him, to produce a global and non overlapping
28 * labeling. An example is shown in the figures down, here we have
29 * a grid 8x6 divided across four processors each processor label locally its grid points
30 *
31 * \verbatim
32 *
33+--------------------------+
34| 1 2 3 4| 1 2 3 4|
35| | |
36| 5 6 7 8| 5 6 7 8|
37| | |
38| 9 10 11 12| 9 10 11 12|
39+--------------------------+
40|13 14 15| 13 14 15 16 17|
41| | |
42|16 17 18| 18 19 20 21 22|
43| | |
44|19 20 21| 23 24 25 26 27|
45+--------------------------+
46
47 *
48 *
49 * \endverbatim
50 *
51 * To the local relabelling is added an offset to make the local id global and non overlapping
52 *
53 *
54 * \verbatim
55 *
56+--------------------------+
57| 1 2 3 4|23 24 25 26|
58| | |
59| 5 6 7 8|27 28 29 30|
60| | |
61| 9 10 12 13|31 32 33 34|
62+--------------------------+
63|14 15 16| 35 36 37 38 39|
64| | |
65|17 18 19| 40 41 42 43 44|
66| | |
67|20 21 22| 45 46 47 48 49|
68+--------------------------+
69 *
70 *
71 * \endverbatim
72 *
73 * \tparam Sys_eqs Definition of the system of equations
74 *
75 * # Examples
76 *
77 * ## Solve lid-driven cavity 2D for incompressible fluid (inertia=0 --> Re=0)
78 *
79 * In this case the system of equation to solve is
80 *
81 * \f$
82\left\{
83\begin{array}{c}
84\eta\nabla v_x + \partial_x P = 0 \quad Eq1 \\
85\eta\nabla v_y + \partial_y P = 0 \quad Eq2 \\
86\partial_x v_x + \partial_y v_y = 0 \quad Eq3
87\end{array}
88\right. \f$
89
90and boundary conditions
91
92 * \f$
93\left\{
94\begin{array}{c}
95v_x = 0, v_y = 0 \quad x = 0 \quad B1\\
96v_x = 0, v_y = 1.0 \quad x = L \quad B2\\
97v_x = 0, v_y = 0 \quad y = 0 \quad B3\\
98v_x = 0, v_y = 0 \quad y = L \quad B4\\
99\end{array}
100\right. \f$
101
102 *
103 * with \f$v_x\f$ and \f$v_y\f$ the velocity in x and y and \f$P\f$ Pressure
104 *
105 * In order to solve such system first we define the general properties of the system
106 *
107 * \snippet eq_unit_test.hpp Definition of the system
108 *
109 * ## Define the equations of the system
110 *
111 * \snippet eq_unit_test.hpp Definition of the equation of the system in the bulk and at the boundary
112 *
113 * ## Define the domain and impose the equations
114 *
115 * \snippet eq_unit_test.hpp lid-driven cavity 2D
116 *
117 * # 3D
118 *
119 * A 3D case is given in the examples
120 *
121 */
122
123template<typename Sys_eqs>
124class FDScheme
125{
126public:
127
128 //! Distributed grid map
129 typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,aggregate<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
130
131 //! Type that specify the properties of the system of equations
132 typedef Sys_eqs Sys_eqs_typ;
133
134 //! Encapsulation of the b term as constant
135 struct constant_b
136 {
137 //! scalar
138 typename Sys_eqs::stype scal;
139
140 /*! \brief Constrictor from a scalar
141 *
142 * \param scal scalar
143 *
144 */
145 constant_b(typename Sys_eqs::stype scal)
146 {
147 this->scal = scal;
148 }
149
150 /*! \brief Get the b term on a grid point
151 *
152 * \note It does not matter the grid point it is a scalar
153 *
154 * \param key grid position (unused because it is a constant)
155 *
156 * \return the scalar
157 *
158 */
159 typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
160 {
161 return scal;
162 }
163 };
164
165 //! Encapsulation of the b term as grid
166 template<typename grid, unsigned int prp>
167 struct grid_b
168 {
169 //! b term fo the grid
170 grid & gr;
171
172 /*! \brief gr grid that encapsulate
173 *
174 * \param gr grid
175 *
176 */
177 grid_b(grid & gr)
178 :gr(gr)
179 {}
180
181 /*! \brief Get the value of the b term on a grid point
182 *
183 * \param key grid point
184 *
185 * \return the value
186 *
187 */
188 auto get(grid_dist_key_dx<Sys_eqs::dims> & key) -> decltype(gr.template get<prp>(key))
189 {
190 return gr.template get<prp>(key);
191 }
192 };
193
194private:
195
196 //! Padding
197 Padding<Sys_eqs::dims> pd;
198
199 //! Sparse matrix triplet type
200 typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;
201
202 //! Vector b
203 typename Sys_eqs::Vector_type b;
204
205 //! Domain Grid informations
206 const grid_sm<Sys_eqs::dims,void> & gs;
207
208 //! Get the grid spacing
209 typename Sys_eqs::stype spacing[Sys_eqs::dims];
210
211 //! mapping grid
212 g_map_type g_map;
213
214 //! row of the matrix
215 size_t row;
216
217 //! row on b
218 size_t row_b;
219
220 //! Grid points that has each processor
221 openfpm::vector<size_t> pnt;
222
223 //! Staggered position for each property
224 comb<Sys_eqs::dims> s_pos[Sys_eqs::nvar];
225
226 //! Each point in the grid has a global id, to decompose correctly the Matrix each processor contain a
227 //! contiguos range of global id, example processor 0 can have from 0 to 234 and processor 1 from 235 to 512
228 //! no processors can have holes in the sequence, this number indicate where the sequence start for this
229 //! processor
230 size_t s_pnt;
231
232 /*! \brief Equation id + space position
233 *
234 */
235 struct key_and_eq
236 {
237 //! space position
238 grid_key_dx<Sys_eqs::dims> key;
239
240 //! equation id
241 size_t eq;
242 };
243
244 /*! \brief From the row Matrix position to the spatial position
245 *
246 * \param row Matrix row
247 *
248 * \return spatial position + equation id
249 *
250 */
251 inline key_and_eq from_row_to_key(size_t row)
252 {
253 key_and_eq ke;
254 auto it = g_map.getDomainIterator();
255
256 while (it.isNext())
257 {
258 size_t row_low = g_map.template get<0>(it.get());
259
260 if (row >= row_low * Sys_eqs::nvar && row < row_low * Sys_eqs::nvar + Sys_eqs::nvar)
261 {
262 ke.eq = row - row_low * Sys_eqs::nvar;
263 ke.key = g_map.getGKey(it.get());
264 return ke;
265 }
266
267 ++it;
268 }
269 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " the row does not map to any position" << "\n";
270
271 return ke;
272 }
273
274 /*! \brief calculate the mapping grid size with padding
275 *
276 * \param sz original grid size
277 * \param pd padding
278 *
279 * \return padded grid size
280 *
281 */
282 inline const std::vector<size_t> padding( const size_t (& sz)[Sys_eqs::dims], Padding<Sys_eqs::dims> & pd)
283 {
284 std::vector<size_t> g_sz_pad(Sys_eqs::dims);
285
286 for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
287 g_sz_pad[i] = sz[i] + pd.getLow(i) + pd.getHigh(i);
288
289 return g_sz_pad;
290 }
291
292 /*! \brief Check if the Matrix is consistent
293 *
294 */
295 void consistency()
296 {
297 openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
298
299 // A and B must have the same rows
300 if (row != row_b)
301 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << "the term B and the Matrix A for Ax=B must contain the same number of rows\n";
302
303 // Indicate all the non zero rows
304 openfpm::vector<unsigned char> nz_rows;
305 nz_rows.resize(row_b);
306
307 for (size_t i = 0 ; i < trpl.size() ; i++)
308 nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
309
310 // Indicate all the non zero colums
311 // This check can be done only on single processor
312
313 Vcluster<> & v_cl = create_vcluster();
314 if (v_cl.getProcessingUnits() == 1)
315 {
316 openfpm::vector<unsigned> nz_cols;
317 nz_cols.resize(row_b);
318
319 for (size_t i = 0 ; i < trpl.size() ; i++)
320 nz_cols.get(trpl.get(i).col()) = true;
321
322 // all the rows must have a non zero element
323 for (size_t i = 0 ; i < nz_rows.size() ; i++)
324 {
325 if (nz_rows.get(i) == false)
326 {
327 key_and_eq ke = from_row_to_key(i);
328 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i << " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
329 }
330 }
331
332 // all the colums must have a non zero element
333 for (size_t i = 0 ; i < nz_cols.size() ; i++)
334 {
335 if (nz_cols.get(i) == false)
336 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
337 }
338 }
339 }
340
341 /*! \brief Copy a given solution vector in a staggered grid
342 *
343 * \tparam Vct Vector type
344 * \tparam Grid_dst target grid
345 * \tparam pos set of properties
346 *
347 * \param v Vector
348 * \param g_dst target staggered grid
349 *
350 */
351 template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_staggered(Vct & v, Grid_dst & g_dst)
352 {
353 // check that g_dst is staggered
354 if (g_dst.is_staggered() == false)
355 std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;
356
357#ifdef SE_CLASS1
358
359 if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
360 std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << std::endl;
361#endif
362
363 // sub-grid iterator over the grid map
364 auto g_map_it = g_map.getDomainIterator();
365
366 // Iterator over the destination grid
367 auto g_dst_it = g_dst.getDomainIterator();
368
369 while (g_map_it.isNext() == true)
370 {
371 typedef typename to_boost_vmpl<pos...>::type vid;
372 typedef boost::mpl::size<vid> v_size;
373
374 auto key_src = g_map_it.get();
375 size_t lin_id = g_map.template get<0>(key_src);
376
377 // destination point
378 auto key_dst = g_dst_it.get();
379
380 // Transform this id into an id for the Eigen vector
381
382 copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());
383
384 boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
385
386 ++g_map_it;
387 ++g_dst_it;
388 }
389 }
390
391 /*! \brief Copy a given solution vector in a normal grid
392 *
393 * \tparam Vct Vector type
394 * \tparam Grid_dst target grid
395 * \tparam pos set of property
396 *
397 * \param v Vector
398 * \param g_dst target normal grid
399 *
400 */
401 template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_normal(Vct & v, Grid_dst & g_dst)
402 {
403 // check that g_dst is staggered
404 if (g_dst.is_staggered() == true)
405 std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be normal " << std::endl;
406
407 grid_key_dx<Grid_dst::dims> start;
408 grid_key_dx<Grid_dst::dims> stop;
409
410 for (size_t i = 0 ; i < Grid_dst::dims ; i++)
411 {
412 start.set_d(i,pd.getLow(i));
413 stop.set_d(i,g_map.size(i) - pd.getHigh(i));
414 }
415
416 // sub-grid iterator over the grid map
417 auto g_map_it = g_map.getSubDomainIterator(start,stop);
418
419 // Iterator over the destination grid
420 auto g_dst_it = g_dst.getDomainIterator();
421
422 while (g_dst_it.isNext() == true)
423 {
424 typedef typename to_boost_vmpl<pos...>::type vid;
425 typedef boost::mpl::size<vid> v_size;
426
427 auto key_src = g_map_it.get();
428 size_t lin_id = g_map.template get<0>(key_src);
429
430 // destination point
431 auto key_dst = g_dst_it.get();
432
433 // Transform this id into an id for the vector
434
435 copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());
436
437 boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
438
439 ++g_map_it;
440 ++g_dst_it;
441 }
442 }
443
444 /*! \brief Impose an operator
445 *
446 * This function the RHS no matrix is imposed. This
447 * function is usefull if the Matrix has been constructed and only
448 * the right hand side b must be changed
449 *
450 * Ax = b
451 *
452 * \param num right hand side of the term (b term)
453 * \param id Equation id in the system that we are imposing
454 * \param it_d iterator that define where you want to impose
455 *
456 */
457 template<typename bop, typename iterator>
458 void impose_dit_b(bop num,
459 long int id ,
460 const iterator & it_d)
461 {
462
463 auto it = it_d;
464 grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
465
466 // iterate all the grid points
467 while (it.isNext())
468 {
469 // get the position
470 auto key = it.get();
471
472 b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
473
474 // if SE_CLASS1 is defined check the position
475#ifdef SE_CLASS1
476// T::position(key,gs,s_pos);
477#endif
478
479 ++row_b;
480 ++it;
481 }
482 }
483
484 /*! \brief Impose an operator
485 *
486 * This function impose an operator on a particular grid region to produce the system
487 *
488 * Ax = b
489 *
490 * ## Stokes equation 2D, lid driven cavity with one splipping wall
491 * \snippet eq_unit_test.hpp Copy the solution to grid
492 *
493 * \param op Operator to impose (A term)
494 * \param num right hand side of the term (b term)
495 * \param id Equation id in the system that we are imposing
496 * \param it_d iterator that define where you want to impose
497 *
498 */
499 template<typename T, typename bop, typename iterator> void impose_dit(const T & op,
500 bop num,
501 long int id ,
502 const iterator & it_d)
503 {
504 openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
505
506 auto it = it_d;
507 grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
508
509 std::unordered_map<long int,typename Sys_eqs::stype> cols;
510
511 // iterate all the grid points
512 while (it.isNext())
513 {
514 // get the position
515 auto key = it.get();
516
517 // Add padding
518 for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
519 key.getKeyRef().set_d(i,key.getKeyRef().get(i) + pd.getLow(i));
520
521 // Calculate the non-zero colums
522 T::value(g_map,key,gs,spacing,cols,1.0);
523
524 // indicate if the diagonal has been set
525 bool is_diag = false;
526
527 // create the triplet
528 for ( auto it = cols.begin(); it != cols.end(); ++it )
529 {
530 trpl.add();
531 trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
532 trpl.last().col() = it->first;
533 trpl.last().value() = it->second;
534
535 if (trpl.last().row() == trpl.last().col())
536 {is_diag = true;}
537
538 }
539
540 // If does not have a diagonal entry put it to zero
541 if (is_diag == false)
542 {
543 trpl.add();
544 trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
545 trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
546 trpl.last().value() = 0.0;
547 }
548
549 b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
550
551 cols.clear();
552
553 // if SE_CLASS1 is defined check the position
554#ifdef SE_CLASS1
555// T::position(key,gs,s_pos);
556#endif
557
558 ++row;
559 ++row_b;
560 ++it;
561 }
562 }
563
564
565
566 /*! \brief Construct the gmap structure
567 *
568 */
569 void construct_gmap()
570 {
571 Vcluster<> & v_cl = create_vcluster();
572
573 // Calculate the size of the local domain
574 size_t sz = g_map.getLocalDomainSize();
575
576 // Get the total size of the local grids on each processors
577 v_cl.allGather(sz,pnt);
578 v_cl.execute();
579 s_pnt = 0;
580
581 // calculate the starting point for this processor
582 for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
583 s_pnt += pnt.get(i);
584
585 // resize b if needed
586 b.resize(Sys_eqs::nvar * g_map.size(),Sys_eqs::nvar * sz);
587
588 // Calculate the starting point
589
590 // Counter
591 size_t cnt = 0;
592
593 // Create the re-mapping grid
594 auto it = g_map.getDomainIterator();
595
596 while (it.isNext())
597 {
598 auto key = it.get();
599
600 g_map.template get<0>(key) = cnt + s_pnt;
601
602 ++cnt;
603 ++it;
604 }
605
606 // sync the ghost
607 g_map.template ghost_get<0>();
608 }
609
610 /*! \initialize the object FDScheme
611 *
612 * \param domain simulation domain
613 *
614 */
615 void Initialize(const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain)
616 {
617 construct_gmap();
618
619 // Create a CellDecomposer and calculate the spacing
620
621 size_t sz_g[Sys_eqs::dims];
622 for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
623 {
624 if (Sys_eqs::boundary[i] == NON_PERIODIC)
625 sz_g[i] = gs.getSize()[i] - 1;
626 else
627 sz_g[i] = gs.getSize()[i];
628 }
629
630 CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);
631
632 for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
633 spacing[i] = cd.getCellBox().getHigh(i);
634 }
635
636public:
637
638 /*! \brief set the staggered position for each property
639 *
640 * \param sp vector containing the staggered position for each property
641 *
642 */
643 void setStagPos(comb<Sys_eqs::dims> (& sp)[Sys_eqs::nvar])
644 {
645 for (size_t i = 0 ; i < Sys_eqs::nvar ; i++)
646 s_pos[i] = sp[i];
647 }
648
649 /*! \brief compute the staggered position for each property
650 *
651 * This is compute from the value_type stored by Sys_eqs::b_grid::value_type
652 * the position of the staggered properties
653 *
654 *
655 */
656 void computeStag()
657 {
658 typedef typename Sys_eqs::b_grid::value_type prp_type;
659
660 openfpm::vector<comb<Sys_eqs::dims>> c_prp[prp_type::max_prop];
661
662 stag_set_position<Sys_eqs::dims,prp_type> ssp(c_prp);
663
664 boost::mpl::for_each_ref< boost::mpl::range_c<int,0,prp_type::max_prop> >(ssp);
665 }
666
667 /*! \brief Get the specified padding
668 *
669 * \return the padding specified
670 *
671 */
672 const Padding<Sys_eqs::dims> & getPadding()
673 {
674 return pd;
675 }
676
677 /*! \brief Return the map between the grid index position and the position in the distributed vector
678 *
679 * It is the map explained in the intro of the FDScheme
680 *
681 * \return the map
682 *
683 */
684 const g_map_type & getMap()
685 {
686 return g_map;
687 }
688
689 /*! \brief Constructor
690 *
691 * The periodicity is given by the grid b_g
692 *
693 * \param stencil maximum extension of the stencil on each directions
694 * \param domain the domain
695 * \param b_g object grid that will store the solution
696 *
697 */
698 FDScheme(const Ghost<Sys_eqs::dims,long int> & stencil,
699 const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
700 const typename Sys_eqs::b_grid & b_g)
701 :pd({0,0,0},{0,0,0}),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
702 {
703 Initialize(domain);
704 }
705
706 /*! \brief Constructor
707 *
708 * The periodicity is given by the grid b_g
709 *
710 * \param pd Padding, how many points out of boundary are present
711 * \param stencil maximum extension of the stencil on each directions
712 * \param domain the domain
713 * \param b_g object grid that will store the solution
714 *
715 */
716 FDScheme(Padding<Sys_eqs::dims> & pd,
717 const Ghost<Sys_eqs::dims,long int> & stencil,
718 const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
719 const typename Sys_eqs::b_grid & b_g)
720 :pd(pd),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
721 {
722 Initialize(domain);
723 }
724
725 /*! \brief Impose an operator
726 *
727 * This function impose an operator on a box region to produce the system
728 *
729 * Ax = b
730 *
731 * ## Stokes equation, lid driven cavity with one splipping wall
732 * \snippet eq_unit_test.hpp lid-driven cavity 2D
733 *
734 * \param op Operator to impose (A term)
735 * \param num right hand side of the term (b term)
736 * \param id Equation id in the system that we are imposing
737 * \param start starting point of the box
738 * \param stop stop point of the box
739 * \param skip_first skip the first point
740 *
741 */
742 template<typename T> void impose(const T & op,
743 typename Sys_eqs::stype num,
744 long int id,
745 const long int (& start)[Sys_eqs::dims],
746 const long int (& stop)[Sys_eqs::dims],
747 bool skip_first = false)
748 {
749 grid_key_dx<Sys_eqs::dims> start_k;
750 grid_key_dx<Sys_eqs::dims> stop_k;
751
752 bool increment = false;
753 if (skip_first == true)
754 {
755 start_k = grid_key_dx<Sys_eqs::dims>(start);
756 stop_k = grid_key_dx<Sys_eqs::dims>(start);
757
758 auto it = g_map.getSubDomainIterator(start_k,stop_k);
759
760 if (it.isNext() == true)
761 increment = true;
762 }
763
764 // add padding to start and stop
765 start_k = grid_key_dx<Sys_eqs::dims>(start);
766 stop_k = grid_key_dx<Sys_eqs::dims>(stop);
767
768
769 auto it = g_map.getSubDomainIterator(start_k,stop_k);
770
771 if (increment == true)
772 ++it;
773
774 constant_b b(num);
775
776 impose_git_gmap(op,b,id,it);
777
778 }
779
780 /*! \brief In case we want to impose a new b re-using FDScheme we have to call
781 * This function
782 */
783 void new_b()
784 {row_b = 0;}
785
786 /*! \brief In case we want to impose a new A re-using FDScheme we have to call
787 * This function
788 *
789 */
790 void new_A()
791 {row = 0;}
792
793 /*! \brief Impose an operator
794 *
795 * This function impose an operator on a particular grid region to produce the system
796 *
797 * Ax = b
798 *
799 * ## Stokes equation 2D, lid driven cavity with one splipping wall
800 * \snippet eq_unit_test.hpp Copy the solution to grid
801 *
802 * \param op Operator to impose (A term)
803 * \param num right hand side of the term (b term)
804 * \param id Equation id in the system that we are imposing
805 * \param it_d iterator that define where you want to impose
806 *
807 */
808 template<unsigned int prp, typename b_term, typename iterator>
809 void impose_dit_b(b_term & b_t,
810 const iterator & it_d,
811 long int id = 0)
812 {
813 grid_b<b_term,prp> b(b_t);
814
815 impose_dit_b(b,id,it_d);
816 }
817
818 /*! \brief Impose an operator
819 *
820 * This function impose an operator on a particular grid region to produce the system
821 *
822 * Ax = b
823 *
824 * ## Stokes equation 2D, lid driven cavity with one splipping wall
825 * \snippet eq_unit_test.hpp Copy the solution to grid
826 *
827 * \param op Operator to impose (A term)
828 * \param num right hand side of the term (b term)
829 * \param id Equation id in the system that we are imposing
830 * \param it_d iterator that define where you want to impose
831 *
832 */
833 template<typename T> void impose_dit(const T & op ,
834 typename Sys_eqs::stype num,
835 long int id ,
836 grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d)
837 {
838 constant_b b(num);
839
840 impose_dit(op,b,id,it_d);
841 }
842
843 /*! \brief Impose an operator
844 *
845 * This function impose an operator on a particular grid region to produce the system
846 *
847 * Ax = b
848 *
849 * ## Stokes equation 2D, lid driven cavity with one splipping wall
850 * \snippet eq_unit_test.hpp Copy the solution to grid
851 *
852 * \param op Operator to impose (A term)
853 * \param num right hand side of the term (b term)
854 * \param id Equation id in the system that we are imposing
855 * \param it_d iterator that define where you want to impose
856 *
857 */
858 template<typename T, typename bop, typename iterator> void impose_git_gmap(const T & op ,
859 bop num,
860 long int id ,
861 const iterator & it_d)
862 {
863 openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
864
865 auto it = it_d;
866 grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
867
868 std::unordered_map<long int,float> cols;
869
870 // iterate all the grid points
871 while (it.isNext())
872 {
873 // get the position
874 auto key = it.get();
875
876 // Calculate the non-zero colums
877 T::value(g_map,key,gs,spacing,cols,1.0);
878
879 // indicate if the diagonal has been set
880 bool is_diag = false;
881
882 // create the triplet
883 for ( auto it = cols.begin(); it != cols.end(); ++it )
884 {
885 trpl.add();
886 trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
887 trpl.last().col() = it->first;
888 trpl.last().value() = it->second;
889
890 if (trpl.last().row() == trpl.last().col())
891 is_diag = true;
892
893// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
894 }
895
896 // If does not have a diagonal entry put it to zero
897 if (is_diag == false)
898 {
899 trpl.add();
900 trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
901 trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
902 trpl.last().value() = 0.0;
903 }
904
905 b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
906
907 cols.clear();
908
909 // if SE_CLASS1 is defined check the position
910#ifdef SE_CLASS1
911// T::position(key,gs,s_pos);
912#endif
913
914 ++row;
915 ++row_b;
916 ++it;
917 }
918 }
919
920 /*! \brief Impose an operator
921 *
922 * This function impose an operator on a particular grid region to produce the system
923 *
924 * Ax = b
925 *
926 * ## Stokes equation 2D, lid driven cavity with one splipping wall
927 * \snippet eq_unit_test.hpp Copy the solution to grid
928 *
929 * \param op Operator to impose (A term)
930 * \param num right hand side of the term (b term)
931 * \param id Equation id in the system that we are imposing
932 * \param it_d iterator that define where you want to impose
933 *
934 */
935 template<typename T, typename bop, typename iterator> void impose_git(const T & op ,
936 bop num,
937 long int id ,
938 const iterator & it_d)
939 {
940 openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
941
942 auto start = it_d.getStart();
943 auto stop = it_d.getStop();
944
945 auto itg = g_map.getSubDomainIterator(start,stop);
946
947 auto it = it_d;
948 grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
949
950 std::unordered_map<long int,float> cols;
951
952 // iterate all the grid points
953 while (it.isNext())
954 {
955 // get the position
956 auto key = it.get();
957 auto keyg = itg.get();
958
959 // Calculate the non-zero colums
960 T::value(g_map,keyg,gs,spacing,cols,1.0);
961
962 // indicate if the diagonal has been set
963 bool is_diag = false;
964
965 // create the triplet
966 for ( auto it = cols.begin(); it != cols.end(); ++it )
967 {
968 trpl.add();
969 trpl.last().row() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
970 trpl.last().col() = it->first;
971 trpl.last().value() = it->second;
972
973 if (trpl.last().row() == trpl.last().col())
974 is_diag = true;
975
976// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
977 }
978
979 // If does not have a diagonal entry put it to zero
980 if (is_diag == false)
981 {
982 trpl.add();
983 trpl.last().row() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
984 trpl.last().col() = g_map.template get<0>(keyg)*Sys_eqs::nvar + id;
985 trpl.last().value() = 0.0;
986 }
987
988 b(g_map.template get<0>(keyg)*Sys_eqs::nvar + id) = num.get(key);
989
990 cols.clear();
991
992 // if SE_CLASS1 is defined check the position
993#ifdef SE_CLASS1
994// T::position(key,gs,s_pos);
995#endif
996
997 ++row;
998 ++row_b;
999 ++it;
1000 ++itg;
1001 }
1002 }
1003
1004 /*! \brief Impose an operator
1005 *
1006 * This function impose an operator on a particular grid region to produce the system
1007 *
1008 * Ax = b
1009 *
1010 * ## Stokes equation 2D, lid driven cavity with one splipping wall
1011 * \snippet eq_unit_test.hpp Copy the solution to grid
1012 *
1013 * \param op Operator to impose (A term)
1014 * \param num right hand side of the term (b term)
1015 * \param id Equation id in the system that we are imposing
1016 * \param it_d iterator that define where you want to impose
1017 *
1018 */
1019 template<unsigned int prp, typename T, typename b_term, typename iterator> void impose_dit(const T & op ,
1020 b_term & b_t,
1021 const iterator & it_d,
1022 long int id = 0)
1023 {
1024 grid_b<b_term,prp> b(b_t);
1025
1026 impose_dit(op,b,id,it_d);
1027 }
1028
1029 //! type of the sparse matrix
1030 typename Sys_eqs::SparseMatrix_type A;
1031
1032 /*! \brief produce the Matrix
1033 *
1034 * \return the Sparse matrix produced
1035 *
1036 */
1037 typename Sys_eqs::SparseMatrix_type & getA()
1038 {
1039#ifdef SE_CLASS1
1040 consistency();
1041#endif
1042 A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar,
1043 g_map.getLocalDomainSize()*Sys_eqs::nvar,
1044 g_map.getLocalDomainSize()*Sys_eqs::nvar);
1045
1046 return A;
1047
1048 }
1049
1050 /*! \brief produce the B vector
1051 *
1052 * \return the vector produced
1053 *
1054 */
1055 typename Sys_eqs::Vector_type & getB()
1056 {
1057#ifdef SE_CLASS1
1058 consistency();
1059#endif
1060
1061 return b;
1062 }
1063
1064 /*! \brief Copy the vector into the grid
1065 *
1066 * ## Copy the solution into the grid
1067 * \snippet eq_unit_test.hpp Copy the solution to grid
1068 *
1069 * \tparam scheme Discretization scheme
1070 * \tparam Grid_dst type of the target grid
1071 * \tparam pos target properties
1072 *
1073 * \param v Vector that contain the solution of the system
1074 * \param start point
1075 * \param stop point
1076 * \param g_dst Destination grid
1077 *
1078 */
1079 template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v,const long int (& start)[Sys_eqs_typ::dims], const long int (& stop)[Sys_eqs_typ::dims], Grid_dst & g_dst)
1080 {
1081 if (is_grid_staggered<Sys_eqs>::value())
1082 {
1083 if (g_dst.is_staggered() == true)
1084 copy_staggered<Vct,Grid_dst,pos...>(v,g_dst);
1085 else
1086 {
1087 // Create a temporal staggered grid and copy the data there
1088 auto & g_map = this->getMap();
1089
1090 // Convert the ghost in grid units
1091
1092 Ghost<Grid_dst::dims,long int> g_int;
1093 for (size_t i = 0 ; i < Grid_dst::dims ; i++)
1094 {
1095 g_int.setLow(i,g_map.getDecomposition().getGhost().getLow(i) / g_map.spacing(i));
1096 g_int.setHigh(i,g_map.getDecomposition().getGhost().getHigh(i) / g_map.spacing(i));
1097 }
1098
1099 staggered_grid_dist<Grid_dst::dims,typename Grid_dst::stype,typename Grid_dst::value_type,typename Grid_dst::decomposition::extended_type, typename Grid_dst::memory_type, typename Grid_dst::device_grid_type> stg(g_dst,g_int,this->getPadding());
1100 stg.setDefaultStagPosition();
1101 copy_staggered<Vct,decltype(stg),pos...>(v,stg);
1102
1103 // sync the ghost and interpolate to the normal grid
1104 stg.template ghost_get<pos...>();
1105 stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
1106 }
1107 }
1108 else
1109 {
1110 copy_normal<Vct,Grid_dst,pos...>(v,g_dst);
1111 }
1112 }
1113
1114 /*! \brief Copy the vector into the grid
1115 *
1116 * ## Copy the solution into the grid
1117 * \snippet eq_unit_test.hpp Copy the solution to grid
1118 *
1119 * \tparam scheme Discretization scheme
1120 * \tparam Grid_dst type of the target grid
1121 * \tparam pos target properties
1122 *
1123 * \param v Vector that contain the solution of the system
1124 * \param g_dst Destination grid
1125 *
1126 */
1127 template<unsigned int ... pos, typename Vct, typename Grid_dst> void copy(Vct & v, Grid_dst & g_dst)
1128 {
1129 long int start[Sys_eqs_typ::dims];
1130 long int stop[Sys_eqs_typ::dims];
1131
1132 for (size_t i = 0 ; i < Sys_eqs_typ::dims ; i++)
1133 {
1134 start[i] = 0;
1135 stop[i] = g_dst.size(i);
1136 }
1137
1138 this->copy<pos...>(v,start,stop,g_dst);
1139 }
1140};
1141
1142#define EQS_FIELDS 0
1143#define EQS_SPACE 1
1144
1145
1146#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_ */
1147