diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index 06dbcf55bae762ec9895d917634450a99d703a48..30b0877fbb557627b0f1f8950a0b676f324296db 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -400,6 +400,9 @@ public:
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
@@ -486,9 +489,6 @@ public:
std::unordered_map<long int,float> cols;
- // resize b if needed
- b.resize(Sys_eqs::nvar * g_map.size());
-
bool is_first = skip_first;
// iterate all the grid points
@@ -509,8 +509,10 @@ public:
// Calculate the non-zero colums
T::value(g_map,key,gs,spacing,cols,1.0);
- // create the triplet
+ // 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();
@@ -518,9 +520,21 @@ public:
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;
cols.clear();
@@ -549,7 +563,7 @@ public:
#ifdef SE_CLASS1
consistency();
#endif
- A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*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;
diff --git a/src/FiniteDifference/eq_unit_test_3d.hpp b/src/FiniteDifference/eq_unit_test_3d.hpp
index 2acac414925ba41531882f250e72f15621ce1993..0ca6d6cb5e0d6457775d179676c685fd214ff5c0 100644
--- a/src/FiniteDifference/eq_unit_test_3d.hpp
+++ b/src/FiniteDifference/eq_unit_test_3d.hpp
@@ -24,7 +24,7 @@ BOOST_AUTO_TEST_SUITE( eq_test_suite_3d )
//!
-struct lid_nn_3d
+struct lid_nn_3d_eigen
{
// dimensionaly of the equation ( 3D problem ...)
static const unsigned int dims = 3;
@@ -50,90 +50,56 @@ struct lid_nn_3d
static const int grid_type = STAGGERED_GRID;
};
-const bool lid_nn_3d::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
-
-// Constant Field
-struct eta
+struct lid_nn_3d_petsc
{
- typedef void const_field;
-
- static float val() {return 1.0;}
-};
-
-// Model the equations
-
-constexpr unsigned int v[] = {0,1,2};
-constexpr unsigned int P = 3;
-constexpr unsigned int ic = 3;
-
-typedef Field<v[x],lid_nn_3d> v_x;
-typedef Field<v[y],lid_nn_3d> v_y;
-typedef Field<v[z],lid_nn_3d> v_z;
-typedef Field<P,lid_nn_3d> Prs;
-
-// Eq1 V_x
-
-typedef mul<eta,Lap<v_x,lid_nn_3d>,lid_nn_3d> eta_lap_vx;
-typedef D<x,Prs,lid_nn_3d> p_x;
-typedef minus<p_x,lid_nn_3d> m_p_x;
-typedef sum<eta_lap_vx,m_p_x,lid_nn_3d> vx_eq;
-
-// Eq2 V_y
-
-typedef mul<eta,Lap<v_y,lid_nn_3d>,lid_nn_3d> eta_lap_vy;
-typedef D<y,Prs,lid_nn_3d> p_y;
-typedef minus<p_y,lid_nn_3d> m_p_y;
-typedef sum<eta_lap_vy,m_p_y,lid_nn_3d> vy_eq;
-
-// Eq3 V_z
-
-typedef mul<eta,Lap<v_z,lid_nn_3d>,lid_nn_3d> eta_lap_vz;
-typedef D<z,Prs,lid_nn_3d> p_z;
-typedef minus<p_z,lid_nn_3d> m_p_z;
-typedef sum<eta_lap_vz,m_p_z,lid_nn_3d> vz_eq;
-
-// Eq4 Incompressibility
+ // dimensionaly of the equation ( 3D problem ...)
+ static const unsigned int dims = 3;
+ // number of fields in the system
+ static const unsigned int nvar = 4;
-typedef D<x,v_x,lid_nn_3d,FORWARD> dx_vx;
-typedef D<y,v_y,lid_nn_3d,FORWARD> dy_vy;
-typedef D<z,v_z,lid_nn_3d,FORWARD> dz_vz;
-typedef sum<dx_vx,dy_vy,dz_vz,lid_nn_3d> ic_eq;
+ // boundary at X and Y
+ static const bool boundary[];
+ // type of space float, double, ...
+ typedef float stype;
-// Directional Avg
-typedef Avg<x,v_y,lid_nn_3d> avg_x_vy;
-typedef Avg<z,v_y,lid_nn_3d> avg_z_vy;
+ // type of base grid
+ typedef grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> b_grid;
-typedef Avg<y,v_x,lid_nn_3d> avg_y_vx;
-typedef Avg<z,v_x,lid_nn_3d> avg_z_vx;
+ // type of SparseMatrix for the linear solver
+ typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
-typedef Avg<y,v_z,lid_nn_3d> avg_y_vz;
-typedef Avg<x,v_z,lid_nn_3d> avg_x_vz;
+ // type of Vector for the linear solver
+ typedef Vector<double,PETSC_BASE> Vector_type;
-// Directional Avg
+ // Define the underline grid is staggered
+ static const int grid_type = STAGGERED_GRID;
+};
-typedef Avg<x,v_y,lid_nn_3d,FORWARD> avg_x_vy_f;
-typedef Avg<z,v_y,lid_nn_3d,FORWARD> avg_z_vy_f;
+//typedef lid_nn_3d_eigen lid_nn_3d;
-typedef Avg<y,v_x,lid_nn_3d,FORWARD> avg_y_vx_f;
-typedef Avg<z,v_x,lid_nn_3d,FORWARD> avg_z_vx_f;
+const bool lid_nn_3d_eigen::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+const bool lid_nn_3d_petsc::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
-typedef Avg<y,v_z,lid_nn_3d,FORWARD> avg_y_vz_f;
-typedef Avg<x,v_z,lid_nn_3d,FORWARD> avg_x_vz_f;
+// Constant Field
+struct eta
+{
+ typedef void const_field;
-#define EQ_1 0
-#define EQ_2 1
-#define EQ_3 2
-#define EQ_4 3
+ static float val() {return 1.0;}
+};
-// Lid driven cavity, uncompressible fluid
+//#define SYSEQ_TYPE lid_nn_3d_eigen;
+//#include "Equations/stoke_flow_eq_3d.hpp"
-BOOST_AUTO_TEST_CASE(lid_driven_cavity)
+template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
{
+ #include "Equations/stoke_flow_eq_3d.hpp"
+
Vcluster & v_cl = create_vcluster();
- if (v_cl.getProcessingUnits() > 3)
- return;
+// if (v_cl.getProcessingUnits() > 3)
+// return;
// Domain
Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
@@ -141,7 +107,7 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Ghost
Ghost<3,float> g(0.01);
- long int sz[] = {32,16,16};
+ long int sz[] = {96,64,64};
size_t szu[3];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
@@ -228,14 +194,20 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
- auto x = umfpack_solver<double>::solve(fd.getA(),fd.getB());
+ std::cout << "Solving " << "\n";
+
+ solver_type solver;
+ solver.best_solve();
+ auto x = solver.solve(fd.getA(),fd.getB());
+
+ std::cout << "Solved" << "\n";
// Bring the solution to grid
- fd.copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
+ fd.template 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()));
- if (v_cl.getProcessUnitID() == 0)
+/* if (v_cl.getProcessUnitID() == 0)
{
if (v_cl.getProcessingUnits() == 1)
{
@@ -261,7 +233,15 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
- }
+ }*/
+}
+
+// Lid driven cavity, uncompressible fluid
+
+BOOST_AUTO_TEST_CASE(lid_driven_cavity)
+{
+// lid_driven_cavity_3d<umfpack_solver<double>,lid_nn_3d_eigen>();
+ lid_driven_cavity_3d<petsc_solver<double>,lid_nn_3d_petsc>();
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Matrix/SparseMatrix.hpp b/src/Matrix/SparseMatrix.hpp
index 6afb15a4599920f1b5239ecfd0e01bf0949ee20e..28ff0ef1c08fced384ac918034e0d9fb90d24bf8 100644
--- a/src/Matrix/SparseMatrix.hpp
+++ b/src/Matrix/SparseMatrix.hpp
@@ -10,9 +10,9 @@
#ifdef HAVE_EIGEN
#include <Eigen/Sparse>
-#define DEFAULT_MATRIX = Eigen::SparseMatrix<T,0,id_t>
+#define DEFAULT_MATRIX = EIGEN_BASE
#else
-#define DEFAULT_MATRIX = void
+#define DEFAULT_MATRIX = 0
#endif
/*! \brief It store the non zero elements of the matrix
@@ -69,7 +69,7 @@ template<typename T, int impl> struct triplet
* \tparam Mi implementation
*
*/
-template<typename T,typename id_t ,typename Mi DEFAULT_MATRIX>
+template<typename T,typename id_t ,unsigned int Mi DEFAULT_MATRIX>
class SparseMatrix
{
public:
@@ -103,4 +103,8 @@ public:
#include "SparseMatrix_Eigen.hpp"
#endif
+#ifdef HAVE_PETSC
+#include "SparseMatrix_petsc.hpp"
+#endif
+
#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_HPP_ */
diff --git a/src/Matrix/SparseMatrix_Eigen.hpp b/src/Matrix/SparseMatrix_Eigen.hpp
index dde9edd02fbee6e663182f4ab27f1b179b59320a..c919690d4dcb419cae79330000c68d01151dada7 100644
--- a/src/Matrix/SparseMatrix_Eigen.hpp
+++ b/src/Matrix/SparseMatrix_Eigen.hpp
@@ -74,15 +74,15 @@ public:
*
*/
template<typename T, typename id_t>
-class SparseMatrix<T,id_t,Eigen::SparseMatrix<T,0,id_t>>
+class SparseMatrix<T,id_t,EIGEN_BASE>
{
public:
//! Triplet implementation id
- typedef boost::mpl::int_<EIGEN_TRIPLET> triplet_impl;
+ typedef boost::mpl::int_<EIGEN_BASE> triplet_impl;
//! Triplet type
- typedef triplet<T,EIGEN_TRIPLET> triplet_type;
+ typedef triplet<T,EIGEN_BASE> triplet_type;
private:
@@ -188,7 +188,7 @@ public:
* \param col number of colums
*
*/
- void resize(size_t row, size_t col)
+ void resize(size_t row, size_t col,size_t l_row, size_t l_col)
{
mat.resize(row,col);
}
diff --git a/src/Matrix/SparseMatrix_unit_tests.hpp b/src/Matrix/SparseMatrix_unit_tests.hpp
index 4c9154636e67397ba6c46ed6d0d7b27ad8d17944..6b93bf6b9969624df668af4fbc582d88b6553e14 100644
--- a/src/Matrix/SparseMatrix_unit_tests.hpp
+++ b/src/Matrix/SparseMatrix_unit_tests.hpp
@@ -11,6 +11,11 @@
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
+#include "Solvers/petsc_solver.hpp"
+
+#ifdef HAVE_PETSC
+#include <petscmat.h>
+#endif
BOOST_AUTO_TEST_SUITE( sparse_matrix_test_suite )
@@ -163,7 +168,246 @@ BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_parallel)
}
}
+BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc)
+{
+ Vcluster & vcl = create_vcluster();
+
+ if (vcl.getProcessingUnits() != 3)
+ return;
+
+ // 3 Processors 9x9 Matrix to invert
+
+ SparseMatrix<double,int,PETSC_BASE> sm(9,9,3);
+ Vector<double,PETSC_BASE> v(9,3);
+
+ typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
+
+ auto & triplets = sm.getMatrixTriplets();
+
+ if (vcl.getProcessUnitID() == 0)
+ {
+ // row 1
+ triplets.add(triplet(0,0,-2));
+ triplets.add(triplet(0,1,1));
+
+ // row 2
+ triplets.add(triplet(1,0,1));
+ triplets.add(triplet(1,1,-2));
+ triplets.add(triplet(1,2,1));
+
+ // row 3
+ triplets.add(triplet(2,1,1));
+ triplets.add(triplet(2,2,-2));
+ triplets.add(triplet(2,3,1));
+
+ // v row1
+ v.insert(0,1);
+ v.insert(1,1);
+ v.insert(2,1);
+ }
+ else if (vcl.getProcessUnitID() == 1)
+ {
+ // row 4
+ triplets.add(triplet(3,2,1));
+ triplets.add(triplet(3,3,-2));
+ triplets.add(triplet(3,4,1));
+
+ // row 5
+ triplets.add(triplet(4,3,1));
+ triplets.add(triplet(4,4,-2));
+ triplets.add(triplet(4,5,1));
+
+ // row 6
+ triplets.add(triplet(5,4,1));
+ triplets.add(triplet(5,5,-2));
+ triplets.add(triplet(5,6,1));
+
+ v.insert(3,1);
+ v.insert(4,1);
+ v.insert(5,1);
+ }
+ else if (vcl.getProcessUnitID() == 2)
+ {
+ // row 7
+ triplets.add(triplet(6,5,1));
+ triplets.add(triplet(6,6,-2));
+ triplets.add(triplet(6,7,1));
+
+ // row 8
+ triplets.add(triplet(7,6,1));
+ triplets.add(triplet(7,7,-2));
+ triplets.add(triplet(7,8,1));
+
+ // row 9
+ triplets.add(triplet(8,7,1));
+ triplets.add(triplet(8,8,-2));
+
+ v.insert(6,1);
+ v.insert(7,1);
+ v.insert(8,1);
+ }
+
+ // Get the petsc Matrix
+ Mat & matp = sm.getMat();
+
+ if (vcl.getProcessUnitID() == 0)
+ {
+ PetscInt r0[1] = {0};
+ PetscInt r1[1] = {1};
+ PetscInt r2[1] = {2};
+ PetscInt r0cols[3] = {0,1};
+ PetscInt r1cols[3] = {0,1,2};
+ PetscInt r2cols[3] = {1,2,3};
+ PetscScalar y[3];
+
+ MatGetValues(matp,1,r0,2,r0cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],-2);
+ BOOST_REQUIRE_EQUAL(y[1],1);
+
+ MatGetValues(matp,1,r1,3,r1cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ MatGetValues(matp,1,r2,3,r2cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ }
+ else if (vcl.getProcessUnitID() == 1)
+ {
+ PetscInt r0[1] = {3};
+ PetscInt r1[1] = {4};
+ PetscInt r2[1] = {5};
+ PetscInt r0cols[3] = {2,3,4};
+ PetscInt r1cols[3] = {3,4,5};
+ PetscInt r2cols[3] = {4,5,6};
+ PetscScalar y[3];
+
+ MatGetValues(matp,1,r0,3,r0cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ MatGetValues(matp,1,r1,3,r1cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ MatGetValues(matp,1,r2,3,r2cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+ }
+ else if (vcl.getProcessUnitID() == 2)
+ {
+ PetscInt r0[1] = {6};
+ PetscInt r1[1] = {7};
+ PetscInt r2[1] = {8};
+ PetscInt r0cols[3] = {5,6,7};
+ PetscInt r1cols[3] = {6,7,8};
+ PetscInt r2cols[3] = {7,8};
+ PetscScalar y[3];
+
+ MatGetValues(matp,1,r0,3,r0cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ MatGetValues(matp,1,r1,3,r1cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+ BOOST_REQUIRE_EQUAL(y[2],1);
+
+ MatGetValues(matp,1,r2,2,r2cols,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],1);
+ BOOST_REQUIRE_EQUAL(y[1],-2);
+
+ }
+
+ petsc_solver<double> solver;
+ solver.solve(sm,v);
+}
+
+BOOST_AUTO_TEST_CASE(sparse_matrix_eigen_petsc_solve)
+{
+ Vcluster & vcl = create_vcluster();
+
+ const int loc = 2000000;
+
+ // 3 Processors 9x9 Matrix to invert
+
+ SparseMatrix<double,int,PETSC_BASE> sm(3*loc,3*loc,loc);
+ Vector<double,PETSC_BASE> v(3*loc,loc);
+
+ typedef SparseMatrix<double,int,PETSC_BASE>::triplet_type triplet;
+
+ auto & triplets = sm.getMatrixTriplets();
+
+ if (vcl.getProcessUnitID() == 0)
+ {
+ // row 1
+ triplets.add(triplet(0,0,-2));
+ triplets.add(triplet(0,1,1));
+ v.insert(0,1);
+ // row 2
+ for (size_t i = 1 ; i < loc ; i++)
+ {
+ triplets.add(triplet(i,i-1,1));
+ triplets.add(triplet(i,i,-2));
+ triplets.add(triplet(i,i+1,1));
+
+ v.insert(i,1);
+ }
+ }
+ else if (vcl.getProcessUnitID() == 1)
+ {
+ // row 2
+ for (size_t i = loc ; i < 2*loc ; i++)
+ {
+ triplets.add(triplet(i,i-1,1));
+ triplets.add(triplet(i,i,-2));
+ triplets.add(triplet(i,i+1,1));
+
+ v.insert(i,1);
+ }
+ }
+ else if (vcl.getProcessUnitID() == 2)
+ {
+ // row 2
+ for (size_t i = 2*loc ; i < 3*loc-1 ; i++)
+ {
+ triplets.add(triplet(i,i-1,1));
+ triplets.add(triplet(i,i,-2));
+ triplets.add(triplet(i,i+1,1));
+
+ v.insert(i,1);
+ }
+
+ // row 9
+ triplets.add(triplet(3*loc-1,3*loc-2,1));
+ triplets.add(triplet(3*loc-1,3*loc-1,-2));
+ v.insert(3*loc-1,1);
+ }
+
+ // Get the petsc Matrix
+ Mat & matp = sm.getMat();
+
+
+ petsc_solver<double> solver;
+ solver.solve(sm,v);
+}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Solvers/umfpack_solver.hpp b/src/Solvers/umfpack_solver.hpp
index 9964728cdbb900b0685184b4b4a5254f31f55b8b..988412c5f2fca1e36974aa6cdb1a2f8455e3bedc 100644
--- a/src/Solvers/umfpack_solver.hpp
+++ b/src/Solvers/umfpack_solver.hpp
@@ -50,7 +50,7 @@ public:
* \tparam impl Implementation of the SparseMatrix
*
*/
- template<typename impl> static Vector<double> solve(SparseMatrix<double,int,impl> & A, const Vector<double> & b, size_t opt = UMFPACK_NONE)
+ static Vector<double,EIGEN_BASE> solve(SparseMatrix<double,int,EIGEN_BASE> & A, const Vector<double,EIGEN_BASE> & b, size_t opt = UMFPACK_NONE)
{
Vcluster & vcl = create_vcluster();
diff --git a/src/Vector/Vector.hpp b/src/Vector/Vector.hpp
index 61973c8dc715f28c806d43055003e8d83ca3817f..453ae1e78084c83a82975908a98b59d5ee443441 100644
--- a/src/Vector/Vector.hpp
+++ b/src/Vector/Vector.hpp
@@ -8,6 +8,8 @@
#ifndef OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
#define OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_
+#include "util/linalgebra_lib.hpp"
+
/*! \brief It store one row value of a vector
*
* Given a row, store a value
@@ -20,17 +22,19 @@ class rval
};
+// Eigen::Matrix<T, Eigen::Dynamic, 1>
+
#ifdef HAVE_EIGEN
#include <Eigen/Sparse>
-#define DEFAULT_VECTOR = Eigen::Matrix<T, Eigen::Dynamic, 1>
+#define DEFAULT_VECTOR = EIGEN_BASE
#else
-#define DEFAULT_VECTOR = void
+#define DEFAULT_VECTOR = 0
#endif
/* ! \brief Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support
*
*/
-template<typename T,typename impl DEFAULT_VECTOR >
+template<typename T,unsigned int impl DEFAULT_VECTOR >
class Vector
{
T stub;
@@ -60,10 +64,7 @@ public:
#endif
#ifdef HAVE_PETSC
-#define OFPM_PETSC_VEC Vec
#include "Vector_petsc.hpp"
-#else
-#define OFPM_PETSC_VEC void
#endif
#endif /* OPENFPM_NUMERICS_SRC_VECTOR_VECTOR_HPP_ */
diff --git a/src/Vector/Vector_eigen.hpp b/src/Vector/Vector_eigen.hpp
index 66c9333eea373d27e656dea2d8c3d1212c1bd6b3..3febf779f86c94cfd7dc05606e738c10d764678c 100644
--- a/src/Vector/Vector_eigen.hpp
+++ b/src/Vector/Vector_eigen.hpp
@@ -68,7 +68,7 @@ public:
};
template<typename T>
-class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
+class Vector<T, EIGEN_BASE>
{
mutable Eigen::Matrix<T, Eigen::Dynamic, 1> v;
@@ -162,7 +162,7 @@ public:
*/
Vector(size_t n)
{
- resize(n);
+ resize(n,0);
}
/*! \brief Create a vector with 0 elements
@@ -175,9 +175,10 @@ public:
/*! \brief Resize the Vector
*
* \param row numbers of row
+ * \param l_row unused
*
*/
- void resize(size_t row)
+ void resize(size_t row, size_t l_row)
{
v.resize(row);
}
diff --git a/src/Vector/Vector_petsc.hpp b/src/Vector/Vector_petsc.hpp
index 4bfff2fb024eedabc5adf189a5cbb00598c01ffd..4fa2d1bd317add2f2dac8684f9b7e4607a2ccc13 100644
--- a/src/Vector/Vector_petsc.hpp
+++ b/src/Vector/Vector_petsc.hpp
@@ -71,19 +71,19 @@ constexpr unsigned int row_id = 0;
constexpr unsigned int val_id = 1;
template<typename T>
-class Vector<T,Vec>
+class Vector<T,PETSC_BASE>
{
// n_row
size_t n_row;
- // n_row_local
+ // Number of local rows
size_t n_row_local;
// Mutable vector
mutable Vec v;
// Mutable row value vector
- openfpm::vector<rval<T,PETSC_RVAL>,HeapMemory,typename memory_traits_inte<rval<T,PETSC_RVAL>>::type > row_val;
+ mutable openfpm::vector<rval<T,PETSC_RVAL>,HeapMemory,typename memory_traits_inte<rval<T,PETSC_RVAL>>::type > row_val;
// Global to local map
mutable std::unordered_map<size_t,size_t> map;
@@ -100,13 +100,15 @@ class Vector<T,Vec>
// Create the vector
PETSC_SAFE_CALL(VecCreate(PETSC_COMM_WORLD,&v));
PETSC_SAFE_CALL(VecSetSizes(v,n_row_local,n_row));
+ PETSC_SAFE_CALL(VecSetFromOptions(v));
// set the vector
- PETSC_SAFE_CALL(VecSetValues(v,n_row_local,&row_val.template get<row_id>(0),&row_val.template get<val_id>(0),INSERT_VALUES))
+ if (row_val.size() != 0)
+ PETSC_SAFE_CALL(VecSetValues(v,row_val.size(),&row_val.template get<row_id>(0),&row_val.template get<val_id>(0),INSERT_VALUES))
-// for (size_t i = 0 ; i < row_val.size() ; i++)
-// v[row_val.get(i).row()] = row_val.get(i).value();
+ VecAssemblyBegin(v);
+ VecAssemblyEnd(v);
}
public:
@@ -116,7 +118,7 @@ public:
* \param v vector to copy
*
*/
- Vector(const Vector<T> & v)
+ Vector(const Vector<T,PETSC_BASE> & v)
{
this->operator=(v);
}
@@ -126,7 +128,7 @@ public:
* \param v vector to copy
*
*/
- Vector(const Vector<T> && v)
+ Vector(const Vector<T,PETSC_BASE> && v)
{
this->operator=(v);
}
@@ -134,28 +136,33 @@ public:
/*! \brief Create a vector with n elements
*
* \param n number of elements in the vector
+ * \param n_row_loc number
*
*/
- Vector(size_t n)
+ Vector(size_t n, size_t n_row_local)
+ :n_row_local(n_row_local)
{
- resize(n);
+ resize(n,n_row_local);
}
/*! \brief Create a vector with 0 elements
*
*/
Vector()
+ :n_row(0),n_row_local(0)
{
}
/*! \brief Resize the Vector
*
* \param row numbers of row
+ * \param l_row number of local row
*
*/
- void resize(size_t row)
+ void resize(size_t row, size_t l_row)
{
n_row = row;
+ n_row_local = l_row;
}
/*! \brief Return a reference to the vector element
@@ -207,8 +214,8 @@ public:
// Map
map[i] = row_val.size()-1;
- row_val.last().row() = i;
- return row_val.last().value();
+ row_val.last().template get<row_id>() = i;
+ return row_val.last().template get<val_id>();
}
/*! \brief Return a reference to the vector element
@@ -253,9 +260,9 @@ public:
return insert(i);
}
- /*! \brief Get the Eigen Vector object
+ /*! \brief Get the PETSC Vector object
*
- * \return the Eigen Vector
+ * \return the PETSC Vector
*
*/
const Vec & getVec() const
@@ -265,9 +272,9 @@ public:
return v;
}
- /*! \brief Get the Eigen Vector object
+ /*! \brief Get the PETSC Vector object
*
- * \return the Eigen Vector
+ * \return the PETSC Vector
*
*/
Vec & getVec()
@@ -277,18 +284,41 @@ public:
return v;
}
-
- /*! \brief Copy the vector
- *
- * \param v vector to copy
+ /*! \brief Update the Vector with the PETSC object
*
*/
- Vector<T,Vec> & operator=(const Vector<T,Vec> & v)
+ void update()
{
- map = v.map;
- row_val = v.row_val;
+ PetscInt n_row;
+ PetscInt n_row_local;
- return *this;
+ // Get the size of the vector from PETSC
+ VecGetSize(v,&n_row);
+ VecGetLocalSize(v,&n_row_local);
+
+ this->n_row = n_row;
+ this->n_row_local = n_row_local;
+
+ row_val.resize(n_row_local);
+
+ //
+
+ PetscInt low;
+ PetscInt high;
+
+ VecGetOwnershipRange(v,&low,&high);
+
+ // Fill the index and construct the map
+
+ size_t k = 0;
+ for (size_t i = low ; i < high ; i++)
+ {
+ row_val.template get<row_id>(k) = i;
+ map[i] = k;
+ k++;
+ }
+
+ PETSC_SAFE_CALL(VecGetValues(v,row_val.size(),&row_val.template get<row_id>(0),&row_val.template get<val_id>(0)))
}
/*! \brief Copy the vector
@@ -296,7 +326,7 @@ public:
* \param v vector to copy
*
*/
- Vector<T,Vec> & operator=(const Vector<T,Vec> && v)
+ Vector<T,PETSC_BASE> & operator=(const Vector<T,PETSC_BASE> & v)
{
map = v.map;
row_val = v.row_val;
@@ -304,17 +334,15 @@ public:
return *this;
}
- /*! \brief Copy the vector (it is used for special purpose)
- *
- * \warning v MUST contain at least all the elements of the vector
+ /*! \brief Copy the vector
*
- * \param v base eigen vector to copy
+ * \param v vector to copy
*
*/
- Vector<T> & operator=(Eigen::Matrix<T, Eigen::Dynamic, 1> & v)
+ Vector<T,PETSC_BASE> & operator=(const Vector<T,PETSC_BASE> && v)
{
- for (size_t i = 0 ; i < row_val.size() ; i++)
- row_val.get(i).value() = v(row_val.get(i).row());
+ map.swap(v.map);
+ row_val.swap(v.row_val);
return *this;
}
diff --git a/src/Vector/Vector_unit_tests.hpp b/src/Vector/Vector_unit_tests.hpp
index 70d8cab740a189536e424a043889676365af3d66..26483a6f814f82ef205474148a516b1fdea0a9f3 100644
--- a/src/Vector/Vector_unit_tests.hpp
+++ b/src/Vector/Vector_unit_tests.hpp
@@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
// 3 Processors 9x9 Matrix to invert
- Vector<double,OFPM_PETSC_VEC> v(9);
+ Vector<double,PETSC_BASE> v(9,3);
if (vcl.getProcessUnitID() == 0)
{
@@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
v.insert(8,8);
}
- Vector<double,OFPM_PETSC_VEC> v3;
+ Vector<double,PETSC_BASE> v3;
v3 = v;
if (vcl.getProcessUnitID() == 0)
@@ -202,9 +202,9 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
BOOST_REQUIRE_EQUAL(v(4),4);
BOOST_REQUIRE_EQUAL(v(5),3);
- BOOST_REQUIRE_EQUAL(v3(0),3);
- BOOST_REQUIRE_EQUAL(v3(1),4);
- BOOST_REQUIRE_EQUAL(v3(2),5);
+ BOOST_REQUIRE_EQUAL(v3(3),3);
+ BOOST_REQUIRE_EQUAL(v3(4),4);
+ BOOST_REQUIRE_EQUAL(v3(5),5);
}
else if (vcl.getProcessUnitID() == 2)
{
@@ -212,12 +212,50 @@ BOOST_AUTO_TEST_CASE(vector_petsc_parallel)
BOOST_REQUIRE_EQUAL(v(7),1);
BOOST_REQUIRE_EQUAL(v(8),0);
- BOOST_REQUIRE_EQUAL(v3(0),6);
- BOOST_REQUIRE_EQUAL(v3(1),7);
- BOOST_REQUIRE_EQUAL(v3(2),9);
+ BOOST_REQUIRE_EQUAL(v3(6),6);
+ BOOST_REQUIRE_EQUAL(v3(7),7);
+ BOOST_REQUIRE_EQUAL(v3(8),8);
+ }
+
+ // Here we get the petsc vector
+ auto & vp = v.getVec();
+
+ // We check the correctness of the PETSC vector
+
+ if (vcl.getProcessUnitID() == 0)
+ {
+ PetscInt ix[] = {0,1,2};
+ PetscScalar y[3];
+
+ VecGetValues(vp,3,ix,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],8);
+ BOOST_REQUIRE_EQUAL(y[1],7);
+ BOOST_REQUIRE_EQUAL(y[2],6);
+ }
+ else if (vcl.getProcessUnitID() == 1)
+ {
+ PetscInt ix[] = {3,4,5};
+ PetscScalar y[3];
+
+ VecGetValues(vp,3,ix,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],5);
+ BOOST_REQUIRE_EQUAL(y[1],4);
+ BOOST_REQUIRE_EQUAL(y[2],3);
+ }
+ else if (vcl.getProcessUnitID() == 2)
+ {
+ PetscInt ix[] = {6,7,8};
+ PetscScalar y[3];
+
+ VecGetValues(vp,3,ix,y);
+
+ BOOST_REQUIRE_EQUAL(y[0],2);
+ BOOST_REQUIRE_EQUAL(y[1],1);
+ BOOST_REQUIRE_EQUAL(y[2],0);
}
- auto & v2 = v.getVec();
}
BOOST_AUTO_TEST_SUITE_END()