From ed634ca877baf88c26b76f9f7f0823866fcf3246 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <i-bird@private-incardon-3.mpi-cbg.de>
Date: Wed, 13 Apr 2016 09:55:27 +0200
Subject: [PATCH] Some small refactorization
---
src/FiniteDifference/FDScheme.hpp | 100 +++++++++++++++++++++--
src/FiniteDifference/eq_unit_test.hpp | 26 ++++--
src/FiniteDifference/eq_unit_test_3d.hpp | 5 +-
src/Vector/Vector_eigen.hpp | 83 -------------------
4 files changed, 116 insertions(+), 98 deletions(-)
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 5e5119ec..32b12e3c 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -18,13 +18,13 @@
/*! \brief Finite Differences
*
- * This class is able to discreatize on a Matrix any system of equations producing a linear system of type \f$Ax=B\f$. In order to create a consistent
- * Matrix it is required that each processor must contain a contiguos range on grid points without
- * holes. In order to ensure this, each processor produce a contiguos local labelling of its local
+ * 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
+ * Matrix it is required that each processor must contain a contiguous range on grid points without
+ * holes. In order to ensure this, each processor produce a contiguous local labeling of its local
* points. Each processor also add an offset equal to the number of local
* points of the processors with id smaller than him, to produce a global and non overlapping
- * labelling. An example is shown in the figures down, here we have
- * a grid 8x6 divided across two processor each processor label locally its grid points
+ * labeling. An example is shown in the figures down, here we have
+ * a grid 8x6 divided across four processors each processor label locally its grid points
*
* \verbatim
*
@@ -266,6 +266,56 @@ private:
}
}
+ /*! \brief Copy a given solution vector in a staggered grid
+ *
+ * \tparam Vct Vector type
+ * \tparam Grid_dst target grid
+ * \tparam pos set of properties
+ *
+ * \param v Vector
+ * \param g_dst target staggered grid
+ *
+ */
+ template<typename Vct, typename Grid_dst ,unsigned int ... pos> void copy_staggered(Vct & v, Grid_dst & g_dst)
+ {
+ // check that g_dst is staggered
+ if (g_dst.is_staggered() == false)
+ std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;
+
+#ifdef SE_CLASS1
+
+ if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
+ std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << 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 Eigen 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;
+ }
+ }
+
public:
/*! \brief set the staggered position for each property
@@ -387,6 +437,7 @@ public:
* Ax = b
*
* ## Stokes equation, lid driven cavity with one splipping wall
+ * \snippet eq_unit_test.hpp lid-driven cavity 2D
*
* \param op Operator to impose (A term)
* \param num right hand side of the term (b term)
@@ -412,7 +463,8 @@ public:
*
* Ax = b
*
- * ## Stokes equation, lid driven cavity with one splipping wall
+ * ## 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)
@@ -513,6 +565,42 @@ public:
return b;
}
+
+ /*! \brief Copy the vector into the grid
+ *
+ * ## Copy the solution into the grid
+ * \snippet eq_unit_test.hpp Copy the solution to grid
+ *
+ * \tparam scheme Discretization scheme
+ * \tparam Grid_dst type of the target grid
+ * \tparam pos target properties
+ *
+ * \param scheme Discretization scheme
+ * \param start point
+ * \param stop point
+ * \param g_dst Destination grid
+ *
+ */
+ 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)
+ {
+ if (is_grid_staggered<Sys_eqs>::value())
+ {
+ if (g_dst.is_staggered() == true)
+ copy_staggered<Vct,Grid_dst,pos...>(v,g_dst);
+ else
+ {
+ // 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());
+ stg.setDefaultStagPosition();
+ copy_staggered<Vct,decltype(stg),pos...>(v,stg);
+
+ // sync the ghost and interpolate to the normal grid
+ stg.template ghost_get<pos...>();
+ stg.template to_normal<Grid_dst,pos...>(g_dst,this->getPadding(),start,stop);
+ }
+ }
+ }
};
#define EQS_FIELDS 0
diff --git a/src/FiniteDifference/eq_unit_test.hpp b/src/FiniteDifference/eq_unit_test.hpp
index 8aff09c9..511e5d01 100644
--- a/src/FiniteDifference/eq_unit_test.hpp
+++ b/src/FiniteDifference/eq_unit_test.hpp
@@ -39,18 +39,18 @@ struct lid_nn
typedef float stype;
// type of base grid, it is the distributed grid that will store the result
- // Note the first property is a 2D vector (velocity), the second is a scalar
+ // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
typedef grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
- // that you are using
+ // that you are using, in case of umfpack it is the only possible choice
typedef SparseMatrix<double,int> SparseMatrix_type;
// type of Vector for the linear system, this parameter is bounded by the solver
- // that you are using
+ // that you are using, in case of umfpack it is the only possible choice
typedef Vector<double> Vector_type;
- // Define that the underline grid where we discretize the operators is staggered
+ // Define that the underline grid where we discretize the system of equation is staggered
static const int grid_type = STAGGERED_GRID;
};
@@ -142,8 +142,15 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
Vcluster & v_cl = *global_v_cluster;
+ if (v_cl.getProcessingUnits() > 3)
+ return;
+
//! [lid-driven cavity 2D]
+ // velocity in the grid is the property 0, pressure is the property 1
+ constexpr int velocity = 0;
+ constexpr int pressure = 1;
+
// Domain, a rectangle
Box<2,float> domain({0.0,0.0},{3.0,1.0});
@@ -168,7 +175,7 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Distributed grid that store the solution
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
- // Ghost stencil
+ // It is the maximum extension of the stencil
Ghost<2,long int> stencil_max(1);
// Finite difference scheme
@@ -219,11 +226,14 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
- // Copy the solution to grid
- x.copy<FDScheme<lid_nn>,decltype(g_dist),0,1>(fd,{0,0},{sz[0]-1,sz[1]-1},g_dist);
-
//! [lid-driven cavity 2D]
+ //! [Copy the solution to grid]
+
+ fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
+
+ //! [Copy the solution to grid]
+
g_dist.write("lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()));
if (v_cl.getProcessUnitID() == 0)
diff --git a/src/FiniteDifference/eq_unit_test_3d.hpp b/src/FiniteDifference/eq_unit_test_3d.hpp
index 2bae893d..cb6defb1 100644
--- a/src/FiniteDifference/eq_unit_test_3d.hpp
+++ b/src/FiniteDifference/eq_unit_test_3d.hpp
@@ -132,6 +132,9 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
Vcluster & v_cl = *global_v_cluster;
+ if (v_cl.getProcessingUnits() > 3)
+ return;
+
// Domain
Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
@@ -231,7 +234,7 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
// Bring the solution to grid
- x.copy<FDScheme<lid_nn_3d>,decltype(g_dist),velocity,pressure>(fd,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
+ fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
g_dist.write("lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()));
diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp
index 57e1c690..b0d161e7 100644
--- a/src/Vector/Vector_eigen.hpp
+++ b/src/Vector/Vector_eigen.hpp
@@ -88,53 +88,6 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
//size of each chunk
mutable openfpm::vector<size_t> sz;
- /*! \brief Copy in a staggered grid
- *
- *
- */
- template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy_staggered_to_staggered(scheme & sc, Grid_dst & g_dst)
- {
- // get the map
- const auto & g_map = sc.getMap();
-
- // check that g_dst is staggered
- if (is_grid_staggered<Grid_dst>::value() == false)
- std::cerr << __FILE__ << ":" << __LINE__ << " The destination grid must be staggered " << std::endl;
-
-#ifdef SE_CLASS1
-
- if (g_map.getLocalDomainSize() != g_dst.getLocalDomainSize())
- std::cerr << __FILE__ << ":" << __LINE__ << " The staggered and destination grid in size does not match " << 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 Eigen vector
-
- copy_ele<typename scheme::Sys_eqs_typ,Grid_dst,typename std::remove_reference<decltype(*this)>::type> cp(key_dst,g_dst,*this,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 Here we collect the full matrix on master
*
@@ -350,42 +303,6 @@ public:
return v;
}
- /*! \brief Copy the vector into the grid
- *
- * ## Copy from the vector to the destination grid
- * \snippet eq_unit_tests.hpp
- *
- * \tparam scheme Discretization scheme
- * \tparam Grid_dst type of the target grid
- * \tparam pos target properties
- *
- * \param scheme Discretization scheme
- * \param start point
- * \param stop point
- * \param g_dst Destination grid
- *
- */
- template<typename scheme, typename Grid_dst ,unsigned int ... pos> void copy(scheme & sc,const long int (& start)[scheme::Sys_eqs_typ::dims], const long int (& stop)[scheme::Sys_eqs_typ::dims], Grid_dst & g_dst)
- {
- if (is_grid_staggered<typename scheme::Sys_eqs_typ>::value())
- {
- if (g_dst.is_staggered() == true)
- copy_staggered_to_staggered<scheme,Grid_dst,pos...>(sc,g_dst);
- else
- {
- // Create a temporal staggered grid and copy the data there
- auto & g_map = sc.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(),sc.getPadding());
- stg.setDefaultStagPosition();
- copy_staggered_to_staggered<scheme,decltype(stg),pos...>(sc,stg);
-
- // sync the ghost and interpolate to the normal grid
- stg.template ghost_get<pos...>();
- stg.template to_normal<Grid_dst,pos...>(g_dst,sc.getPadding(),start,stop);
- }
- }
- }
-
/*! \brief Scatter the vector information to the other processors
*
* Eigen does not have a real parallel vector, so in order to work we have to scatter
--
GitLab