Commit 55add63f authored by Pietro Incardona's avatar Pietro Incardona

Testing all PETSC_SOLVERS

parent 71fc012b
......@@ -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;
......
......@@ -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()
......
......@@ -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_ */
......@@ -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);
}
......
......@@ -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()
......
......@@ -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();
......
......@@ -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_ */
......@@ -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);
}
......
......@@ -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);