diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index b2e55433876fd5890357de1812bfc9524128bd68..00cac997777bd5d57efcfddb1e3986cf34638caf 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -134,6 +134,38 @@ public:
private:
+ // Encapsulation of the b term as constant
+ struct constant_b
+ {
+ typename Sys_eqs::stype scal;
+
+ constant_b(typename Sys_eqs::stype scal)
+ {
+ this->scal = scal;
+ }
+
+ typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
+ {
+ return scal;
+ }
+ };
+
+ // Encapsulation of the b term as grid
+ template<typename grid, unsigned int prp>
+ struct grid_b
+ {
+ grid & gr;
+
+ grid_b(grid & gr)
+ :gr(gr)
+ {}
+
+ typename Sys_eqs::stype get(grid_dist_key_dx<Sys_eqs::dims> & key)
+ {
+ return gr.template get<prp>(key);
+ }
+ };
+
//! Padding
Padding<Sys_eqs::dims> pd;
@@ -247,34 +279,36 @@ private:
nz_rows.resize(row_b);
for (size_t i = 0 ; i < trpl.size() ; i++)
- {
- nz_rows.get(trpl.get(i).row()) = true;
- }
+ nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
// Indicate all the non zero colums
- openfpm::vector<unsigned> nz_cols;
- nz_cols.resize(row_b);
+ // This check can be done only on single processor
- for (size_t i = 0 ; i < trpl.size() ; i++)
+ Vcluster & v_cl = create_vcluster();
+ if (v_cl.getProcessingUnits() == 1)
{
- nz_cols.get(trpl.get(i).col()) = true;
- }
+ openfpm::vector<unsigned> nz_cols;
+ nz_cols.resize(row_b);
- // all the rows must have a non zero element
- for (size_t i = 0 ; i < nz_rows.size() ; i++)
- {
- if (nz_rows.get(i) == false)
+ for (size_t i = 0 ; i < trpl.size() ; i++)
+ nz_cols.get(trpl.get(i).col()) = true;
+
+ // all the rows must have a non zero element
+ for (size_t i = 0 ; i < nz_rows.size() ; i++)
{
- key_and_eq ke = from_row_to_key(i);
- std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i << " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
+ if (nz_rows.get(i) == false)
+ {
+ key_and_eq ke = from_row_to_key(i);
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i << " is not filled, position " << ke.key.to_string() << " equation: " << ke.eq << "\n";
+ }
}
- }
- // all the colums must have a non zero element
- for (size_t i = 0 ; i < nz_cols.size() ; i++)
- {
- if (nz_cols.get(i) == false)
- std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
+ // all the colums must have a non zero element
+ for (size_t i = 0 ; i < nz_cols.size() ; i++)
+ {
+ if (nz_cols.get(i) == false)
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
+ }
}
}
@@ -328,6 +362,178 @@ private:
}
}
+ /*! \brief Copy a given solution vector in a normal grid
+ *
+ * \tparam Vct Vector type
+ * \tparam Grid_dst target grid
+ * \tparam pos set of property
+ *
+ * \param v Vector
+ * \param g_dst target normal grid
+ *
+ */
+ template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_normal(Vct & v, Grid_dst & g_dst)
+ {
+ // check that g_dst is staggered
+ if (g_dst.is_staggered() == true)
+ std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be normal " << std::endl;
+
+#ifdef SE_CLASS1
+
+ if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
+ std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid in size does not match the mapping grid" << std::endl;
+#endif
+
+ // sub-grid iterator over the grid map
+ auto g_map_it = g_map.getDomainIterator();
+
+ // Iterator over the destination grid
+ auto g_dst_it = g_dst.getDomainIterator();
+
+ while (g_map_it.isNext() == true)
+ {
+ typedef typename to_boost_vmpl<pos...>::type vid;
+ typedef boost::mpl::size<vid> v_size;
+
+ auto key_src = g_map_it.get();
+ size_t lin_id = g_map.template get<0>(key_src);
+
+ // destination point
+ auto key_dst = g_dst_it.get();
+
+ // Transform this id into an id for the vector
+
+ copy_ele<Sys_eqs_typ, Grid_dst,Vct> cp(key_dst,g_dst,v,lin_id,g_map.size());
+
+ boost::mpl::for_each_ref<boost::mpl::range_c<int,0,v_size::value>>(cp);
+
+ ++g_map_it;
+ ++g_dst_it;
+ }
+ }
+
+ /*! \brief Impose an operator
+ *
+ * This function impose an operator on a particular grid region to produce the system
+ *
+ * Ax = b
+ *
+ * ## Stokes equation 2D, lid driven cavity with one splipping wall
+ * \snippet eq_unit_test.hpp Copy the solution to grid
+ *
+ * \param op Operator to impose (A term)
+ * \param num right hand side of the term (b term)
+ * \param id Equation id in the system that we are imposing
+ * \param it_d iterator that define where you want to impose
+ *
+ */
+ template<typename T, typename bop, typename iterator> void impose(const T & op ,
+ bop num,
+ long int id ,
+ const iterator & it_d)
+ {
+ openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
+
+ auto it = it_d;
+ grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
+
+ std::unordered_map<long int,float> cols;
+
+ // iterate all the grid points
+ while (it.isNext())
+ {
+ // get the position
+ auto key = it.get();
+
+ // Calculate the non-zero colums
+ T::value(g_map,key,gs,spacing,cols,1.0);
+
+ // indicate if the diagonal has been set
+ bool is_diag = false;
+
+ // create the triplet
+ for ( auto it = cols.begin(); it != cols.end(); ++it )
+ {
+ trpl.add();
+ trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+ trpl.last().col() = it->first;
+ trpl.last().value() = it->second;
+
+ if (trpl.last().row() == trpl.last().col())
+ is_diag = true;
+
+// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
+ }
+
+ // If does not have a diagonal entry put it to zero
+ if (is_diag == false)
+ {
+ trpl.add();
+ trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+ trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
+ trpl.last().value() = 0.0;
+ }
+
+ b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num.get(key);
+
+ cols.clear();
+// std::cout << "\n";
+
+ // if SE_CLASS1 is defined check the position
+#ifdef SE_CLASS1
+// T::position(key,gs,s_pos);
+#endif
+
+ ++row;
+ ++row_b;
+ ++it;
+ }
+ }
+
+ /*! \brief Construct the gmap structure
+ *
+ */
+ void construct_gmap()
+ {
+ Vcluster & v_cl = create_vcluster();
+
+ // Calculate the size of the local domain
+ size_t sz = g_map.getLocalDomainSize();
+
+ // Get the total size of the local grids on each processors
+ v_cl.allGather(sz,pnt);
+ v_cl.execute();
+ s_pnt = 0;
+
+ // calculate the starting point for this processor
+ for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
+ s_pnt += pnt.get(i);
+
+ // resize b if needed
+ b.resize(Sys_eqs::nvar * g_map.size(),Sys_eqs::nvar * sz);
+
+ // Calculate the starting point
+
+ // Counter
+ size_t cnt = 0;
+
+ // Create the re-mapping grid
+ auto it = g_map.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ g_map.template get<0>(key) = cnt + s_pnt;
+
+ ++cnt;
+ ++it;
+ }
+
+ // sync the ghost
+ g_map.template ghost_get<0>();
+ }
+
public:
/*! \brief set the staggered position for each property
@@ -382,12 +588,14 @@ public:
}
/*! \brief Constructor
+ *
+ * The periodicity is given by the grid b_g
*
* \param pd Padding, how many points out of boundary are present
- * \param stencil how many ghost points are required for this calculation
+ * \param stencil maximum extension of the stencil on each directions
* \param domain the domain
* \param gs grid infos where Finite differences work
- * \param b_g maximum extension of the stencil on each directions
+ * \param b_g Finite difference grid (it is not important what it contain only its decomposition is used)
*
*/
FDScheme(Padding<Sys_eqs::dims> & pd,
@@ -397,49 +605,18 @@ public:
const typename Sys_eqs::b_grid & b_g)
:pd(pd),gs(gs),g_map(b_g,stencil,pd),row(0),row_b(0)
{
- Vcluster & v_cl = b_g.getDecomposition().getVC();
-
- // Calculate the size of the local domain
- size_t sz = g_map.getLocalDomainSize();
-
- // Get the total size of the local grids on each processors
- v_cl.allGather(sz,pnt);
- v_cl.execute();
- s_pnt = 0;
-
- // calculate the starting point for this processor
- for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
- s_pnt += pnt.get(i);
-
- // resize b if needed
- b.resize(Sys_eqs::nvar * g_map.size(),Sys_eqs::nvar * sz);
-
- // Calculate the starting point
-
- // Counter
- size_t cnt = 0;
-
- // Create the re-mapping grid
- auto it = g_map.getDomainIterator();
-
- while (it.isNext())
- {
- auto key = it.get();
-
- g_map.template get<0>(key) = cnt + s_pnt;
-
- ++cnt;
- ++it;
- }
-
- // sync the ghost
- g_map.template ghost_get<0>();
+ construct_gmap();
// Create a CellDecomposer and calculate the spacing
size_t sz_g[Sys_eqs::dims];
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
- sz_g[i] = gs.getSize()[i] - 1;
+ {
+ if (Sys_eqs::boundary[i] == NON_PERIODIC)
+ sz_g[i] = gs.getSize()[i] - 1;
+ else
+ sz_g[i] = gs.getSize()[i];
+ }
CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);
@@ -495,7 +672,9 @@ public:
if (increment == true)
++it;
- impose(op,num,id,it);
+ constant_b b(num);
+
+ impose(op,b,id,it);
}
@@ -514,64 +693,39 @@ public:
* \param it_d iterator that define where you want to impose
*
*/
- template<typename T> void impose(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)
+ template<typename T> void impose(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)
{
- openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
-
- auto it = it_d;
- grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
-
- std::unordered_map<long int,float> cols;
-
- // iterate all the grid points
- while (it.isNext())
- {
- // get the position
- auto key = it.get();
-
- // Calculate the non-zero colums
- T::value(g_map,key,gs,spacing,cols,1.0);
-
- // indicate if the diagonal has been set
- bool is_diag = false;
-
- // create the triplet
- for ( auto it = cols.begin(); it != cols.end(); ++it )
- {
- trpl.add();
- trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
- trpl.last().col() = it->first;
- trpl.last().value() = it->second;
-
- if (trpl.last().row() == trpl.last().col())
- is_diag = true;
-
-// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
- }
-
- // If does not have a diagonal entry put it to zero
- if (is_diag == false)
- {
- trpl.add();
- trpl.last().row() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
- trpl.last().col() = g_map.template get<0>(key)*Sys_eqs::nvar + id;
- trpl.last().value() = 0.0;
- }
-
- b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
+ constant_b b(num);
- cols.clear();
-// std::cout << "\n";
+ impose(op,b,id,it_d);
+ }
- // if SE_CLASS1 is defined check the position
-#ifdef SE_CLASS1
-// T::position(key,gs,s_pos);
-#endif
+ /*! \brief Impose an operator
+ *
+ * This function impose an operator on a particular grid region to produce the system
+ *
+ * Ax = b
+ *
+ * ## Stokes equation 2D, lid driven cavity with one splipping wall
+ * \snippet eq_unit_test.hpp Copy the solution to grid
+ *
+ * \param op Operator to impose (A term)
+ * \param num right hand side of the term (b term)
+ * \param id Equation id in the system that we are imposing
+ * \param it_d iterator that define where you want to impose
+ *
+ */
+ template<unsigned int prp, typename T, typename b_term, typename iterator> void impose(const T & op ,
+ b_term & b_t,
+ long int id ,
+ const iterator & it_d)
+ {
+ grid_b<b_term,prp> b(b_t);
- ++row;
- ++row_b;
- ++it;
- }
+ impose(op,b,id,it_d);
}
//! type of the sparse matrix
@@ -587,7 +741,9 @@ public:
#ifdef SE_CLASS1
consistency();
#endif
- A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar,g_map.getLocalDomainSize()*Sys_eqs::nvar,g_map.getLocalDomainSize()*Sys_eqs::nvar);
+ A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar,
+ g_map.getLocalDomainSize()*Sys_eqs::nvar,
+ g_map.getLocalDomainSize()*Sys_eqs::nvar);
return A;
@@ -632,7 +788,17 @@ public:
{
// Create a temporal staggered grid and copy the data there
auto & g_map = this->getMap();
- 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_map.getDecomposition().getGhost(),this->getPadding());
+
+ // Convert the ghost in grid units
+
+ Ghost<Grid_dst::dims,long int> g_int;
+ for (size_t i = 0 ; i < Grid_dst::dims ; i++)
+ {
+ g_int.setLow(i,g_map.getDecomposition().getGhost().getLow(i) / g_map.spacing(i));
+ g_int.setHigh(i,g_map.getDecomposition().getGhost().getHigh(i) / g_map.spacing(i));
+ }
+
+ 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());
stg.setDefaultStagPosition();
copy_staggered<Vct,decltype(stg),pos...>(v,stg);
@@ -641,6 +807,10 @@ public:
stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
}
}
+ else
+ {
+ copy_normal<Vct,Grid_dst,pos...>(v,g_dst);
+ }
}
};
diff --git a/src/Makefile.am b/src/Makefile.am
index d802738e4a88986d8fed6e1501633cc39a4ccaae..bb2f3e1442f77bbb79de0a7160310502071c7b40 100755
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -12,7 +12,8 @@ Solvers/umfpack_solver.hpp Solvers/petsc_solver.hpp \
util/petsc_util.hpp util/linalgebra_lib.hpp util/util_num.hpp util/grid_dist_testing.hpp \
FiniteDifference/Average.hpp FiniteDifference/Derivative.hpp FiniteDifference/eq.hpp FiniteDifference/FDScheme.hpp FiniteDifference/Laplacian.hpp FiniteDifference/mul.hpp FiniteDifference/sum.hpp FiniteDifference/util/common.hpp \
PSE/Kernels.hpp PSE/Kernels_test_util.hpp Operators/Vector/vector_dist_operators_extensions.hpp Operators/Vector/vector_dist_operators.hpp Operators/Vector/vector_dist_operators_apply_kernel.hpp Operators/Vector/vector_dist_operators_functions.hpp Operators/Vector/vector_dist_operator_assign.hpp \
-Draw/DrawParticles.hpp Draw/PointIterator.hpp Draw/PointIteratorSkin.hpp
+Draw/DrawParticles.hpp Draw/PointIterator.hpp Draw/PointIteratorSkin.hpp \
+interpolation/interpolation.hpp interpolation/mp4_kernel.hpp interpolation/z_spline.hpp
.cu.o :
$(NVCC) $(NVCCFLAGS) -o $@ -c $<
diff --git a/src/Matrix/SparseMatrix_petsc.hpp b/src/Matrix/SparseMatrix_petsc.hpp
index 78d352d1ab1d7d4a3921231173b26324edec3724..b9c2191407d21e12ec98456b8793235aa763ecd1 100644
--- a/src/Matrix/SparseMatrix_petsc.hpp
+++ b/src/Matrix/SparseMatrix_petsc.hpp
@@ -103,31 +103,36 @@ public:
private:
- // Size of the matrix
+ //! Number of matrix row (global)
size_t g_row;
+ //! Number of matrix colums (global)
size_t g_col;
- // Local size of the matrix
+ //! Number of matrix row (local)
size_t l_row;
+ //! Number of matrix colums (local)
size_t l_col;
- // starting row for this processor
+ //! starting row for this processor
size_t start_row;
- // indicate if the matrix has been created
+ //! indicate if the matrix has been created
bool m_created = false;
- // PETSC Matrix
+ //! PETSC Matrix
Mat mat;
+
+ //! Triplets of the matrix
openfpm::vector<triplet_type> trpl;
- openfpm::vector<triplet_type> trpl_recv;
- // temporary list of values
+
+ //! temporary list of values
mutable openfpm::vector<PetscScalar> vals;
- // temporary list of colums
+ //! temporary list of colums
mutable openfpm::vector<PetscInt> cols;
- // PETSC d_nnz and o_nnz
+ //! PETSC d_nnz
mutable openfpm::vector<PetscInt> d_nnz;
+ //! PETSC o_nnz
mutable openfpm::vector<PetscInt> o_nnz;
/*! \brief Fill the petsc Matrix
@@ -152,7 +157,7 @@ private:
{
PetscInt row = trpl.get(i).row();
- while(row == trpl.get(i).row() && i < trpl.size())
+ while(i < trpl.size() && row == trpl.get(i).row())
{
if ((size_t)trpl.get(i).col() >= start_row && (size_t)trpl.get(i).col() < start_row + l_row)
d_nnz.get(row - start_row)++;
@@ -176,7 +181,7 @@ private:
PetscInt row = trpl.get(i).row();
- while(row == trpl.get(i).row() && i < trpl.size())
+ while(i < trpl.size() && row == trpl.get(i).row())
{
vals.add(trpl.get(i).value());
cols.add(trpl.get(i).col());
@@ -237,7 +242,7 @@ public:
{PETSC_SAFE_CALL(MatDestroy(&mat));}
}
- /*! \brief Get the Matrix triplets bugger
+ /*! \brief Get the Matrix triplets buffer
*
* It return a buffer that can be filled with triplets
*
diff --git a/src/Solvers/petsc_solver.hpp b/src/Solvers/petsc_solver.hpp
index d7ec9a30db39d9728c1d55402271cf01315e1b9d..795e0f34cf8f603d38b4353f1d257a54a7dfa60f 100644
--- a/src/Solvers/petsc_solver.hpp
+++ b/src/Solvers/petsc_solver.hpp
@@ -16,6 +16,8 @@
#include <petscksp.h>
#include <petsctime.h>
#include "Plot/GoogleChart.hpp"
+#include "Matrix/SparseMatrix.hpp"
+#include "Vector/Vector.hpp"
#define UMFPACK_NONE 0
diff --git a/src/main.cpp b/src/main.cpp
index 6f3c9131de95cec2a6ff080c2bb25f425ee81bb2..fdb08845570c9f532770a5cf3b80a773ec90cecb 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -16,3 +16,4 @@
#include "PSE/Kernels_unit_tests.hpp"
#include "Operators/Vector/vector_dist_operators_unit_tests.hpp"
#include "Draw/DrawParticles_unit_tests.hpp"
+#include "interpolation/interpolation_unit_tests.hpp"