Commit 58c0aea7 authored by incardon's avatar incardon

Fixes on AMG + documentation + AMG report

parent 43a521cb
......@@ -8,11 +8,6 @@
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
#define CENTRAL 0
#define CENTRAL_B_ONE_SIDE 1
#define FORWARD 2
#define BACKWARD 3
#define CENTRAL_SYM 4
#include "util/mathutil.hpp"
#include "Vector/map_vector.hpp"
......@@ -20,6 +15,7 @@
#include "FiniteDifference/util/common.hpp"
#include "util/util_num.hpp"
#include <unordered_map>
#include "FD_util_include.hpp"
/*! \brief Derivative second order on h (spacing)
*
......
......@@ -23,21 +23,23 @@ constexpr unsigned int V = 0;
struct sys_nn
{
// dimensionaly of the equation (2D problem 3D problem ...)
//! dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system
//! number of degree of freedoms
static const unsigned int nvar = 1;
static const unsigned int ord = EQS_FIELD;
// boundary at X and Y
//! boundary at X and Y
static const bool boundary[];
// type of space float, double, ...
//! type of space float, double, ...
typedef float stype;
// Base grid
//! Base grid
typedef grid_dist_id<dims,stype,scalar<float>,CartDecomposition<2,stype> > b_grid;
//! specify that we are on testing
typedef void testing;
};
......@@ -45,21 +47,22 @@ const bool sys_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
struct sys_pp
{
// dimensionaly of the equation (2D problem 3D problem ...)
//! dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system
//! number of degree of freedom in the system
static const unsigned int nvar = 1;
static const unsigned int ord = EQS_FIELD;
// boundary at X and Y
//! boundary at X and Y
static const bool boundary[];
// type of space float, double, ...
//! type of space float, double, ...
typedef float stype;
// Base grid
//! Base grid
typedef grid_dist_id<dims,stype,scalar<float>,CartDecomposition<2,stype> > b_grid;
//! Indicate we are on testing
typedef void testing;
};
......@@ -69,23 +72,26 @@ const bool sys_pp::boundary[] = {PERIODIC,PERIODIC};
struct syss_nn
{
// dimensionaly of the equation (2D problem 3D problem ...)
//! dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system
//! Degree of freedom in the system
static const unsigned int nvar = 1;
static const unsigned int ord = EQS_FIELD;
//! Indicate that the type of grid is staggered
static const unsigned int grid_type = STAGGERED_GRID;
// boundary at X and Y
//! boundary at X and Y
static const bool boundary[];
// type of space float, double, ...
//! type of space float, double, ...
typedef float stype;
// Base grid
//! Base grid
typedef grid_dist_id<dims,stype,scalar<float>,CartDecomposition<2,stype> > b_grid;
//! Indicate we are on testing
typedef void testing;
};
......@@ -93,23 +99,25 @@ const bool syss_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
struct syss_pp
{
// dimensionaly of the equation (2D problem 3D problem ...)
//! dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system
//! number of fields in the system
static const unsigned int nvar = 1;
static const unsigned int ord = EQS_FIELD;
//! Indicate the grid is staggered
static const unsigned int grid_type = STAGGERED_GRID;
// boundary at X and Y
//! boundary at X and Y
static const bool boundary[];
// type of space float, double, ...
//! type of space float, double, ...
typedef float stype;
// Base grid
//! Base grid
typedef grid_dist_id<dims,stype,scalar<float>,CartDecomposition<2,stype> > b_grid;
//! Indicate we are on testing
typedef void testing;
};
......
......@@ -8,6 +8,9 @@
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_LAPLACIAN_HPP_
#include "FD_util_include.hpp"
#include "util/util_num.hpp"
#include "FiniteDifference/eq.hpp"
/*! \brief Laplacian second order on h (spacing)
*
......
......@@ -193,8 +193,7 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
solver_type solver;
solver.best_solve();
auto x = solver.solve(fd.getA(),fd.getB());
auto x = solver.try_solve(fd.getA(),fd.getB());
// Bring the solution to grid
fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
......
......@@ -2,7 +2,7 @@
LINKLIBS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(OPENMP_LDFLAGS) $(LIBHILBERT_LIB) $(PETSC_LIB) $(SUITESPARSE_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(METIS_LIB) $(PARMETIS_LIB) $(DEFAULT_LIB) $(PTHREAD_LIBS) $(OPT_LIBS) $(BOOST_LDFLAGS) $(BOOST_PROGRAM_OPTIONS_LIB) $(BOOST_IOSTREAMS_LIB) $(LIBQUADMATH) $(OPENMP_LDFLAGS) $(LIBIFCORE)
noinst_PROGRAMS = numerics
numerics_SOURCES = main.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
numerics_SOURCES = main.cpp Solvers/petsc_solver_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp
numerics_CXXFLAGS = $(HDF5_CPPFLAGS) $(OPENMP_CFLAGS) $(LIBHILBERT_INCLUDE) $(AM_CXXFLAGS) $(INCLUDES_PATH) $(BOOST_CPPFLAGS) $(SUITESPARSE_INCLUDE) $(METIS_INCLUDE) $(PARMETIS_INCLUDE) $(EIGEN_INCLUDE) $(PETSC_INCLUDE) -Wno-deprecated-declarations -Wno-unused-local-typedefs
numerics_CFLAGS = $(CUDA_CFLAGS)
numerics_LDADD = $(LINKLIBS) -lparmetis -lmetis
......
......@@ -83,7 +83,7 @@ public:
triplet() {};
};
/* ! \brief Sparse Matrix implementation, that map over Eigen
/*! \brief Sparse Matrix implementation, that map over Eigen
*
* \tparam T Type of the sparse Matrix store on each row,colums
* \tparam id_t type of id
......@@ -194,6 +194,8 @@ private:
PETSC_SAFE_CALL(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY));
PETSC_SAFE_CALL(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY));
m_created = true;
}
......@@ -207,7 +209,7 @@ public:
*
*/
SparseMatrix(size_t N1, size_t N2, size_t n_row_local)
:l_row(n_row_local),l_col(n_row_local)
:g_row(N1),g_col(N2),l_row(n_row_local),l_col(n_row_local)
{
PETSC_SAFE_CALL(MatCreate(PETSC_COMM_WORLD,&mat));
PETSC_SAFE_CALL(MatSetType(mat,MATMPIAIJ));
......@@ -230,6 +232,7 @@ public:
*
*/
SparseMatrix()
:g_row(0),g_col(0),l_row(0l),l_col(0),start_row(0)
{
PETSC_SAFE_CALL(MatCreate(PETSC_COMM_WORLD,&mat));
PETSC_SAFE_CALL(MatSetType(mat,MATMPIAIJ));
......@@ -251,6 +254,8 @@ public:
*/
openfpm::vector<triplet_type> & getMatrixTriplets()
{
m_created = false;
return this->trpl;
}
......@@ -261,7 +266,8 @@ public:
*/
const Mat & getMat() const
{
fill_petsc();
if (m_created == false)
{fill_petsc();}
return mat;
}
......@@ -273,7 +279,8 @@ public:
*/
Mat & getMat()
{
fill_petsc();
if (m_created == false)
{fill_petsc();}
return mat;
}
......
......@@ -8,11 +8,11 @@
#ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
#define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_
#include "config.h"
#ifdef HAVE_PETSC
#include "Vector/Vector.hpp"
#include "Eigen/UmfPackSupport"
#include <Eigen/SparseLU>
#include <petscksp.h>
#include <petsctime.h>
#include "Plot/GoogleChart.hpp"
......@@ -20,12 +20,41 @@
#include "Vector/Vector.hpp"
#define UMFPACK_NONE 0
#define REL_DOWN 1
#define REL_UP 2
#define REL_ALL 3
#define PCHYPRE_BOOMERAMG "petsc_boomeramg"
enum AMG_type
{
HYPRE_AMG,
PETSC_AMG,
TRILINOS_ML
};
/*! \brief In case T does not match the PETSC precision compilation create a
* stub structure
*
* \param T precision
*
*/
template<typename T>
class petsc_solver
{
public:
/*! \brief Solve the linear system.
*
* In this case just return an error
*
* \param A sparse matrix
* \param b right-hand-side
*
* \return the solution
*
*/
template<typename impl> static Vector<T> solve(const SparseMatrix<T,impl> & A, const Vector<T> & b)
{
std::cerr << "Error Petsc only suppor double precision" << "/n";
......@@ -36,84 +65,87 @@ public:
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2
// It contain statistic of the error of the calculated solution
//! It contain statistic of the error of the calculated solution
struct solError
{
// infinity norm of the error
//! infinity norm of the error
PetscReal err_inf;
// L1 norm of the error
//! L1 norm of the error
PetscReal err_norm;
// Number of iterations
//! Number of iterations
PetscInt its;
};
size_t cpu_rank;
/*! \brief This class is able to do Matrix inversion in parallel with PETSC solvers
*
* #Example of how to solve a lid driven cavity 3D with a PETSC solver in paralle
* \snippet
* # Example of use of this class #
*
* \snippet eq_unit_test.hpp lid-driven cavity 2D
*
*/
template<>
class petsc_solver<double>
{
//It contain statistic of the error at each iteration
//! contain the infinity norm of the residual at each iteration
struct itError
{
//! Iteration
PetscInt it;
PetscReal err_inf;
PetscReal err_norm;
itError(const solError & e)
{
it = e.its;
err_inf = e.err_inf;
err_norm = e.err_norm;
}
// Default constructor
itError() {}
//! error norm at iteration it
PetscReal err_norm;
};
// It contain the benchmark
/*! \brief It contain the benchmark information for each solver
*
*
*/
struct solv_bench_info
{
// Method name
//! Solver Krylov name
std::string method;
//Method name (short form)
//! Method name (short form)
std::string smethod;
// time to converge in milliseconds
//! time to converge in milliseconds
double time;
// Solution error
//! Solution error
solError err;
// Convergence per iteration
//! Convergence per iteration
openfpm::vector<itError> res;
};
// KSP Maximum number of iterations
//! indicate if the preconditioner is set
bool is_preconditioner_set = false;
//! KSP Maximum number of iterations
PetscInt maxits;
// Main parallel solver
//! Main parallel solver
KSP ksp;
// Temporal variable used for calculation of static members
//! Temporal variable used for calculation of static members
size_t tmp;
// The full set of solvers
//! The full set of solvers
openfpm::vector<std::string> solvs;
bool try_solve = false;
// It contain the solver benchmark results
//! It contain the solver benchmark results
openfpm::vector<solv_bench_info> bench;
//! Type of the algebraic multi-grid preconditioner
AMG_type atype;
//! Block size
int block_sz = 0;
/*! \brief Calculate the residual error at time t for one method
*
* \param t time
......@@ -130,9 +162,9 @@ class petsc_solver<double>
size_t pt = std::floor(t / s_int);
if (pt < slv.res.size())
return slv.res.get(pt).err_inf;
return slv.res.get(pt).err_norm;
return slv.res.last().err_inf;
return slv.res.last().err_norm;
}
/*! \brief Here we write the benchmark report
......@@ -247,6 +279,29 @@ class petsc_solver<double>
cg.write("gc_solver.html");
}
/*! \brief It set up the solver based on the provided options
*
* \param A_ Matrix
* \param x_ solution
* \param b_ right-hand-side
*
*/
void pre_solve_impl(const Mat & A_, const Vec & b_, Vec & x_)
{
PETSC_SAFE_CALL(KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED));
if (atype == HYPRE_AMG)
{
PC pc;
// We set the pre-conditioner
PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
PCFactorSetShiftAmount(pc, PETSC_DECIDE);
}
}
/*! \brief Print a progress bar on standard out
*
*
......@@ -259,72 +314,56 @@ class petsc_solver<double>
std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
}
/*! \brief procedure to collect the absolute residue at each iteration
/*! \brief procedure print the progress of the solver in benchmark mode
*
* \param ksp Solver
* \param it Iteration number
* \param res resudual
* \param custom pointer to data
* \param data custom pointer to data
*
* \return always zero
*
*/
static PetscErrorCode monitor(KSP ksp,PetscInt it,PetscReal res,void* data)
static PetscErrorCode monitor_progress_residual(KSP ksp,PetscInt it,PetscReal res,void* data)
{
petsc_solver<double> * pts = static_cast<petsc_solver<double> *>(data);
Mat A;
Mat P;
Vec x;
Vec b;
petsc_solver<double> * pts = (petsc_solver *)data;
PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,&P));
PETSC_SAFE_CALL(KSPGetSolution(ksp,&x));
PETSC_SAFE_CALL(KSPGetRhs(ksp,&b));
pts->progress(it);
solError err = statSolutionError(A,b,x,ksp);
itError erri(err);
itError erri;
erri.it = it;
itError err_fill;
size_t old_size = pts->bench.last().res.size();
pts->bench.last().res.resize(it+1);
if (old_size > 0)
err_fill = pts->bench.last().res.get(old_size-1);
else
err_fill = erri;
for (long int i = old_size ; i < (long int)it ; i++)
pts->bench.last().res.get(i) = err_fill;
// Add the error per iteration
pts->bench.last().res.get(it) = erri;
Mat A;
Vec Br,v,w;
PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,NULL));
PETSC_SAFE_CALL(MatCreateVecs(A,&w,&v));
PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
PETSC_SAFE_CALL(KSPBuildResidual(ksp,v,w,&Br));
PETSC_SAFE_CALL(VecNorm(Br,NORM_INFINITY,&erri.err_norm));
itError err_fill;
pts->progress(it);
size_t old_size = pts->bench.last().res.size();
pts->bench.last().res.resize(it+1);
return 0;
}
if (old_size > 0)
err_fill = pts->bench.last().res.get(old_size-1);
else
err_fill = erri;
/*! \brief procedure print the progress of the solver in benchmark mode
*
* \param ksp Solver
* \param it Iteration number
* \param res resudual
* \param custom pointer to data
*
*/
static PetscErrorCode monitor_progress(KSP ksp,PetscInt it,PetscReal res,void* data)
{
petsc_solver<double> * pts = (petsc_solver *)data;
for (long int i = old_size ; i < (long int)it ; i++)
pts->bench.last().res.get(i) = err_fill;
pts->progress(it);
// Add the error per iteration
pts->bench.last().res.get(it) = erri;
return 0;
}
/*! \brief This function print an "*" showing the progress of the solvers
*
* \param it iteration number
*
*/
void progress(PetscInt it)
......@@ -368,6 +407,9 @@ class petsc_solver<double>
/*! \brief It convert the KSP type into a human read-able string
*
* \param solv solver (short form)
*
* \return the name of the solver in long form
*
*/
std::string to_string_method(const std::string & solv)
......@@ -453,10 +495,10 @@ class petsc_solver<double>
*
* \param A_ Matrix
* \param b_ vector of coefficents
* \param x solution
* \param x_ solution
*
*/
void try_solve_simple(Mat & A_, const Vec & b_, Vec & x_)
void try_solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
{
Vec best_sol;
PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
......@@ -473,6 +515,14 @@ class petsc_solver<double>
for (size_t i = 0 ; i < solvs.size() ; i++)
{
initKSPForTest();
// Here we solve without preconditioner
PC pc;
// We set the pre-conditioner to none
PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
PETSC_SAFE_CALL(PetscOptionsSetValue("-pc_type",PCNONE));
setSolver(solvs.get(i).c_str());
// Setup for BCGSL, GMRES
......@@ -531,6 +581,7 @@ class petsc_solver<double>
PETSC_SAFE_CALL(VecDestroy(&best_sol));
}
/*! \brief Benchmark solve simple solving x=inv(A)*b
*
* \param A_ Matrix A
......@@ -539,7 +590,7 @@ class petsc_solver<double>
* \param bench structure that store the benchmark information
*
*/
void bench_solve_simple(Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
void bench_solve_simple(const Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
{
// timer for benchmark
timer t;
......@@ -547,7 +598,7 @@ class petsc_solver<double>
// Enable progress monitor
tmp = 0;
PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress,this,NULL));
PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress_residual,this,NULL));
print_progress_bar();
solve_simple(A_,b_,x_);
......@@ -564,22 +615,19 @@ class petsc_solver<double>
bench.err = err;
// Solve getting the residual per iteration
tmp = 0;
PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,this,NULL));
print_progress_bar();
solve_simple(A_,b_,x_);
// New line
if (create_vcluster().getProcessUnitID() == 0)
std::cout << std::endl;
}
/*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
*
* \param A_ SparseMatrix
* \param b_ right-hand-side
* \param x_ solution
*
*/
void solve_simple(Mat & A_, const Vec & b_, Vec & x_)
void solve_simple(const Mat & A_, const Vec & b_, Vec & x_)
{
// PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
......@@ -592,37 +640,41 @@ class petsc_solver<double>
PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
PC pc;
// We set the Matrix operators
PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
// We set the pre-conditioner
PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
PETSC_SAFE_CALL(PCSetType(pc,PCNONE));
// if we are on on best solve set-up a monitor function
if (try_solve == true)
{
// Disable convergence check
PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
}
PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
PETSC_SAFE_CALL(KSPSetUp(ksp));
// Solve the system
PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
}
/*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
*
* \param b_ right hand side
* \param x_ solution
*
*/
void solve_simple(const Vec & b_, Vec & x_)
{
// Solve the system
PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
}
/*! \brief Calculate statistic on the error solution
*
* \brief A Matrix of the system
* \brief b Right hand side of the matrix
* \brief x Solution
* \param A_ Matrix of the system
* \param b_ Right hand side of the matrix
* \param x_ Solution
* \param ksp Krylov solver
*
* \return the solution error
*
*/
static solError statSolutionError(Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
static solError statSolutionError(const Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
{
solError err;
......@@ -654,9 +706,6 @@ class petsc_solver<double>
PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
setMaxIter(maxits);
// Disable convergence check
PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
}
/*! \brief Destroy the KSP object
......@@ -674,8 +723,10 @@ class petsc_solver<double>
* \param x_ the solution
* \param b_ the right-hand-side
*
* \return the solution error
*
*/
static solError getSolNormError(const Mat & A_, const Vec & x_, const Vec & b_)
static solError getSolNormError(const Mat & A_, const Vec & b_, const Vec & x_)
{
PetscScalar neg_one = -1.0;
......@@ -713,6 +764,7 @@ class petsc_solver<double>
public:
//! Type of the solution object
typedef Vector<double,PETSC_BASE> return_type;
~petsc_solver()
......@@ -746,6 +798,8 @@ public:
* The try solve function use the most robust solvers in PETSC, if you want to add
* additionally solver like KSPIBCGS,KSPFBCGSR,KSPPGMRES, use addTestSolver(std::string(KSPIBCGS))
*
* \param solver additional solver solver to test
*
*/
void addTestSolver(std::string & solver)
{
......@@ -757,6 +811,8 @@ public:
* The try solve function use the most robust solvers in PETSC, if you want to remove
* a solver like use removeTestSolver(std::string(KSPIBCGS))
*
* \param solver remove solver to test
*
*/
void removeTestSolver(const std::string & solver)
{
......@@ -769,18 +825,6 @@ public:
}
}
/*! \brief Set the solver to test all the solvers and generate a report
*
* In this mode the system will try different Solvers, Preconditioner and
* combination of solvers in order to find the best solver in speed, and
* precision. As output it will produce a performance report
*
*
*/