8 #ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
9 #define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_FDSCHEME_HPP_
11 #include "../Matrix/SparseMatrix.hpp"
12 #include "Grid/grid_dist_id.hpp"
13 #include "Grid/Iterators/grid_dist_id_iterator_sub.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"
123 template<
typename Sys_eqs>
168 template<
typename gr
id,
unsigned int prp>
192 return gr.template get<prp>(
key);
200 typedef typename Sys_eqs::SparseMatrix_type::triplet_type
triplet;
203 typename Sys_eqs::Vector_type
b;
209 typename Sys_eqs::stype
spacing[Sys_eqs::dims];
258 size_t row_low =
g_map.template get<0>(it.get());
260 if (row >= row_low * Sys_eqs::nvar && row < row_low * Sys_eqs::nvar + Sys_eqs::nvar)
262 ke.
eq = row - row_low * Sys_eqs::nvar;
269 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" the row does not map to any position" <<
"\n";
284 std::vector<size_t> g_sz_pad(Sys_eqs::dims);
286 for (
size_t i = 0 ; i < Sys_eqs::dims ; i++)
301 std::cerr <<
"Error " << __FILE__ <<
":" << __LINE__ <<
"the term B and the Matrix A for Ax=B must contain the same number of rows\n";
305 nz_rows.resize(
row_b);
307 for (
size_t i = 0 ; i < trpl.
size() ; i++)
308 nz_rows.get(trpl.get(i).row() -
s_pnt*Sys_eqs::nvar) =
true;
313 Vcluster & v_cl = create_vcluster();
317 nz_cols.resize(
row_b);
319 for (
size_t i = 0 ; i < trpl.
size() ; i++)
320 nz_cols.get(trpl.get(i).col()) =
true;
323 for (
size_t i = 0 ; i < nz_rows.
size() ; i++)
325 if (nz_rows.get(i) ==
false)
328 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" Ill posed matrix row " << i <<
" is not filled, position " << ke.
key.
to_string() <<
" equation: " << ke.
eq <<
"\n";
333 for (
size_t i = 0 ; i < nz_cols.
size() ; i++)
335 if (nz_cols.get(i) ==
false)
336 std::cerr <<
"Error: " << __FILE__ <<
":" << __LINE__ <<
" Ill posed matrix colum " << i <<
" is not filled\n";
351 template<
typename Vct,
typename Grid_dst ,
unsigned int ... pos>
void copy_staggered(Vct & v, Grid_dst & g_dst)
354 if (g_dst.is_staggered() ==
false)
355 std::cerr << __FILE__ <<
":" << __LINE__ <<
" The destination grid must be staggered " << std::endl;
360 std::cerr << __FILE__ <<
":" << __LINE__ <<
" The staggered and destination grid in size does not match " << std::endl;
367 auto g_dst_it = g_dst.getDomainIterator();
369 while (g_map_it.isNext() ==
true)
372 typedef boost::mpl::size<vid> v_size;
374 auto key_src = g_map_it.get();
375 size_t lin_id =
g_map.template get<0>(key_src);
378 auto key_dst = g_dst_it.
get();
384 boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
401 template<
typename Vct,
typename Grid_dst ,
unsigned int ... pos>
void copy_normal(Vct & v, Grid_dst & g_dst)
404 if (g_dst.is_staggered() ==
true)
405 std::cerr << __FILE__ <<
":" << __LINE__ <<
" The destination grid must be normal " << std::endl;
410 for (
size_t i = 0 ; i < Grid_dst::dims ; i++)
420 auto g_dst_it = g_dst.getDomainIterator();
422 while (g_dst_it.isNext() ==
true)
425 typedef boost::mpl::size<vid> v_size;
427 auto key_src = g_map_it.get();
428 size_t lin_id =
g_map.template get<0>(key_src);
431 auto key_dst = g_dst_it.
get();
437 boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
457 template<
typename bop,
typename iterator>
460 const iterator & it_d)
472 b(
g_map.template get<0>(
key)*Sys_eqs::nvar +
id) = num.get(
key);
499 template<
typename T,
typename bop,
typename iterator>
void impose_dit(
const T & op,
502 const iterator & it_d)
509 std::unordered_map<long int,typename Sys_eqs::stype> cols;
518 for (
size_t i = 0 ; i < Sys_eqs::dims ; i++)
525 bool is_diag =
false;
528 for (
auto it = cols.begin(); it != cols.end(); ++it )
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;
535 if (trpl.last().row() == trpl.last().col())
541 if (is_diag ==
false)
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;
549 b(
g_map.template get<0>(
key)*Sys_eqs::nvar +
id) = num.get(
key);
579 template<
typename T,
typename bop,
typename iterator>
void impose_git(
const T & op ,
582 const iterator & it_d)
589 std::unordered_map<long int,float> cols;
601 bool is_diag =
false;
604 for (
auto it = cols.begin(); it != cols.end(); ++it )
607 trpl.last().row() =
g_map.template get<0>(
key)*Sys_eqs::nvar +
id;
608 trpl.last().col() = it->first;
609 trpl.last().value() = it->second;
611 if (trpl.last().row() == trpl.last().col())
618 if (is_diag ==
false)
621 trpl.last().row() =
g_map.template get<0>(
key)*Sys_eqs::nvar +
id;
622 trpl.last().col() =
g_map.template get<0>(
key)*Sys_eqs::nvar +
id;
623 trpl.last().value() = 0.0;
626 b(
g_map.template get<0>(
key)*Sys_eqs::nvar +
id) = num.get(
key);
646 Vcluster & v_cl = create_vcluster();
661 b.resize(Sys_eqs::nvar *
g_map.
size(),Sys_eqs::nvar * sz);
682 g_map.template ghost_get<0>();
696 size_t sz_g[Sys_eqs::dims];
697 for (
size_t i = 0 ; i < Sys_eqs::dims ; i++)
699 if (Sys_eqs::boundary[i] == NON_PERIODIC)
705 CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);
707 for (
size_t i = 0 ; i < Sys_eqs::dims ; i++)
708 spacing[i] = cd.getCellBox().getHigh(i);
720 for (
size_t i = 0 ; i < Sys_eqs::nvar ; i++)
733 typedef typename Sys_eqs::b_grid::value_type prp_type;
739 boost::mpl::for_each_ref< boost::mpl::range_c<int,0,prp_type::max_prop> >(ssp);
775 const typename Sys_eqs::b_grid & b_g)
776 :
pd({0,0,0},{0,0,0}),
gs(b_g.getGridInfoVoid()),
g_map(b_g,stencil,
pd),
row(0),
row_b(0)
794 const typename Sys_eqs::b_grid & b_g)
795 :pd(pd),
gs(b_g.getGridInfoVoid()),
g_map(b_g,stencil,pd),
row(0),
row_b(0)
817 template<
typename T>
void impose(
const T & op,
818 typename Sys_eqs::stype num,
820 const long int (& start)[Sys_eqs::dims],
821 const long int (& stop)[Sys_eqs::dims],
822 bool skip_first =
false)
827 bool increment =
false;
828 if (skip_first ==
true)
835 if (it.isNext() ==
true)
845 if (increment ==
true)
882 template<
unsigned int prp,
typename b_term,
typename iterator>
884 const iterator & it_d,
908 typename Sys_eqs::stype num,
932 template<
unsigned int prp,
typename T,
typename b_term,
typename iterator>
void impose_dit(
const T & op ,
934 const iterator & it_d,
943 typename Sys_eqs::SparseMatrix_type
A;
950 typename Sys_eqs::SparseMatrix_type &
getA()
968 typename Sys_eqs::Vector_type &
getB()
992 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)
996 if (g_dst.is_staggered() ==
true)
997 copy_staggered<Vct,Grid_dst,pos...>(v,g_dst);
1006 for (
size_t i = 0 ; i < Grid_dst::dims ; i++)
1017 stg.template ghost_get<pos...>();
1018 stg.template to_normal<Grid_dst,pos...>(g_dst,this->
getPadding(),start,stop);
1040 template<
unsigned int ... pos,
typename Vct,
typename Grid_dst>
void copy(Vct & v, Grid_dst & g_dst)
1042 long int start[Sys_eqs_typ::dims];
1043 long int stop[Sys_eqs_typ::dims];
1045 for (
size_t i = 0 ; i < Sys_eqs_typ::dims ; i++)
1048 stop[i] = g_dst.size(i);
1051 this->
copy<pos...>(v,start,stop,g_dst);
1055 #define EQS_FIELDS 0
Padding< Sys_eqs::dims > pd
Padding.
void setDefaultStagPosition()
It set all the properties defined to be staggered on the default location.
grid_key_dx< Sys_eqs::dims > key
space position
void copy_staggered(Vct &v, Grid_dst &g_dst)
Copy a given solution vector in a staggered grid.
size_t row
row of the matrix
void computeStag()
compute the staggered position for each property
auto get(const grid_dist_key_dx< dim > &v1) const -> typename std::add_lvalue_reference< decltype(loc_grid.get(v1.getSub()).template get< p >(v1.getKey()))>::type
Get the reference of the selected element.
grid_dist_iterator< dim, device_grid, FREE > getDomainIterator() const
It return an iterator that span the full grid domain (each processor span its local domain) ...
Sys_eqs::Vector_type b
Vector b.
FDScheme(const Ghost< Sys_eqs::dims, long int > &stencil, const Box< Sys_eqs::dims, typename Sys_eqs::stype > &domain, const typename Sys_eqs::b_grid &b_g)
Constructor.
g_map_type g_map
mapping grid
T getLow(int i) const
get the i-coordinate of the low bound interval of the box
size_t getProcessUnitID()
Get the process unit id.
FDScheme(Padding< Sys_eqs::dims > &pd, const Ghost< Sys_eqs::dims, long int > &stencil, const Box< Sys_eqs::dims, typename Sys_eqs::stype > &domain, const typename Sys_eqs::b_grid &b_g)
Constructor.
const grid_sm< dim, void > & getGridInfoVoid() const
Get an object containing the grid informations without type.
void execute()
Execute all the requests.
comb< Sys_eqs::dims > s_pos[Sys_eqs::nvar]
Staggered position for each property.
Sys_eqs::stype spacing[Sys_eqs::dims]
Get the grid spacing.
grid_b(grid &gr)
gr grid that encapsulate
const grid_sm< Sys_eqs::dims, void > & gs
Domain Grid informations.
void impose_dit(const T &op, typename Sys_eqs::stype num, long int id, grid_dist_iterator_sub< Sys_eqs::dims, typename g_map_type::d_grid > it_d)
Impose an operator.
void copy_normal(Vct &v, Grid_dst &g_dst)
Copy a given solution vector in a normal grid.
T getHigh(int i) const
get the high interval of the box
const Padding< Sys_eqs::dims > & getPadding()
Get the specified padding.
size_t size() const
Return the total number of points in the grid.
void setHigh(int i, T val)
set the high interval of the box
void impose(const T &op, typename Sys_eqs::stype num, long int id, const long int(&start)[Sys_eqs::dims], const long int(&stop)[Sys_eqs::dims], bool skip_first=false)
Impose an operator.
std::string to_string()
convert the information into a string
Grid key for a distributed grid.
Sys_eqs::stype scal
scalar
void new_A()
In case we want to impose a new A re-using FDScheme we have to call This function.
grid_dist_iterator_sub< dim, device_grid > getSubDomainIterator(const grid_key_dx< dim > &start, const grid_key_dx< dim > &stop) const
It return an iterator that span the grid domain only in the specified part.
Sys_eqs::SparseMatrix_type & getA()
produce the Matrix
Encapsulation of the b term as constant.
void setStagPos(comb< Sys_eqs::dims >(&sp)[Sys_eqs::nvar])
set the staggered position for each property
Implementation of VCluster class.
void copy(Vct &v, Grid_dst &g_dst)
Copy the vector into the grid.
void construct_gmap()
Construct the gmap structure.
grid_key_dx< dim > getGKey(const grid_dist_key_dx< dim > &k)
Convert a g_dist_key_dx into a global key.
This is a distributed grid.
void new_b()
In case we want to impose a new b re-using FDScheme we have to call This function.
Implementation of the staggered grid.
constant_b(typename Sys_eqs::stype scal)
Constrictor from a scalar.
this class is a functor for "for_each" algorithm
const g_map_type & getMap()
Return the map between the grid index position and the position in the distributed vector...
const size_t(& getSize() const)[N]
Return the size of the grid as an array.
is_grid_staggered analyse T if it has a property grid_type defined and indicate that the grid is stag...
key_and_eq from_row_to_key(size_t row)
From the row Matrix position to the spatial position.
Distributed grid iterator.
void setLow(int i, T val)
set the low interval of the box
void impose_dit_b(bop num, long int id, const iterator &it_d)
Impose an operator.
void impose_dit(const T &op, b_term &b_t, const iterator &it_d, long int id=0)
Impose an operator.
Encapsulation of the b term as grid.
This class represent an N-dimensional box.
Sys_eqs Sys_eqs_typ
Type that specify the properties of the system of equations.
This class is a trick to indicate the compiler a specific specialization pattern. ...
void consistency()
Check if the Matrix is consistent.
void impose_dit(const T &op, bop num, long int id, const iterator &it_d)
Impose an operator.
openfpm::vector< size_t > pnt
Grid points that has each processor.
const std::vector< size_t > padding(const size_t(&sz)[Sys_eqs::dims], Padding< Sys_eqs::dims > &pd)
calculate the mapping grid size with padding
Sys_eqs::SparseMatrix_type::triplet_type triplet
Sparse matrix triplet type.
Sys_eqs::Vector_type & getB()
produce the B vector
size_t getLocalDomainSize() const
Get the total number of grid points for the calling processor.
void Initialize(const Box< Sys_eqs::dims, typename Sys_eqs::stype > &domain)
Equation id + space position.
void set_d(size_t i, mem_id id)
Set the i index.
void impose_git(const T &op, bop num, long int id, const iterator &it_d)
Impose an operator.
grid_dist_id< Sys_eqs::dims, typename Sys_eqs::stype, aggregate< size_t >, typename Sys_eqs::b_grid::decomposition::extended_type > g_map_type
Distributed grid map.
this class is a functor for "for_each" algorithm
void copy(Vct &v, const long int(&start)[Sys_eqs_typ::dims], const long int(&stop)[Sys_eqs_typ::dims], Grid_dst &g_dst)
Copy the vector into the grid.
size_t getProcessingUnits()
Get the total number of processors.
Decomposition & getDecomposition()
Get the object that store the information about the decomposition.
bool allGather(T &send, openfpm::vector< T, Mem, gr > &v)
Gather the data from all processors.
void impose_dit_b(b_term &b_t, const iterator &it_d, long int id=0)
Impose an operator.
Sys_eqs::SparseMatrix_type A
type of the sparse matrix
grid & gr
b term fo the grid
St spacing(size_t i) const
Get the spacing of the grid in direction i.