diff --git a/src/FiniteDifference/eq_unit_test_3d.hpp b/src/FiniteDifference/eq_unit_test_3d.hpp
index 0ca6d6cb5e0d6457775d179676c685fd214ff5c0..a0939906aa0da3c43bab788f11dbf4eba7dae016 100644
--- a/src/FiniteDifference/eq_unit_test_3d.hpp
+++ b/src/FiniteDifference/eq_unit_test_3d.hpp
@@ -98,8 +98,8 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
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});
@@ -107,7 +107,7 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
// Ghost
Ghost<3,float> g(0.01);
- long int sz[] = {96,64,64};
+ long int sz[] = {36,12,12};
size_t szu[3];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
@@ -194,46 +194,22 @@ template<typename solver_type,typename lid_nn_3d> void lid_driven_cavity_3d()
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});
- 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.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.getProcessingUnits() == 1)
- {
- // Check that match
- bool test = compare("lid_driven_cavity_3d_p1_grid_0_test.vtk","lid_driven_cavity_3d_p1_grid_0.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- }
- else if (v_cl.getProcessingUnits() == 2)
- {
- // Check that match
- bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- }
- else if (v_cl.getProcessingUnits() == 3)
- {
- // Check that match
- bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
- BOOST_REQUIRE_EQUAL(test,true);
- }
- }*/
+ std::string s = std::string(demangle(typeid(solver_type).name()));
+ s += "_";
+
+ g_dist.write(s + "lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()));
+
+ // Check that match
+ bool test = compare(std::string(std::string("test/") + s + "lid_driven_cavity_3d_p3_grid_" + std::to_string(v_cl.getProcessUnitID()) + "_test.vtk"),std::string(s + "lid_driven_cavity_3d_p3_grid_" + std::to_string(v_cl.getProcessUnitID()) + ".vtk"));
+ BOOST_REQUIRE_EQUAL(test,true);
+
}
// Lid driven cavity, uncompressible fluid
diff --git a/src/Solvers/petsc_solver.hpp b/src/Solvers/petsc_solver.hpp
index b499de66fd857d1ac32ea8c97be87080bda655b7..0566014f1ea23b039cfb7da0b2bcc180c9543f1b 100644
--- a/src/Solvers/petsc_solver.hpp
+++ b/src/Solvers/petsc_solver.hpp
@@ -13,9 +13,11 @@
#include <Eigen/SparseLU>
#include <petscksp.h>
#include <petsctime.h>
+#include "Plot/GoogleChart.hpp"
#define UMFPACK_NONE 0
+
template<typename T>
class petsc_solver
{
@@ -31,6 +33,14 @@ public:
#define SOLVER_PRINT_RESIDUAL_NORM_INFINITY 1
#define SOLVER_PRINT_DETERMINANT 2
+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
+ *
+ */
template<>
class petsc_solver<double>
{
@@ -47,6 +57,40 @@ class petsc_solver<double>
PetscInt its;
};
+ //It contain statistic of the error at each iteration
+ struct itError
+ {
+ 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;
+ }
+ };
+
+ // It contain the benchmark
+ struct solv_bench_info
+ {
+ // Method name
+ std::string method;
+
+ //Method name (short form)
+ std::string smethod;
+
+ // time to converge in milliseconds
+ double time;
+
+ // Solution error
+ solError err;
+
+ // Convergence per iteration
+ openfpm::vector<itError> res;
+ };
+
// KSP Relative tolerance
PetscReal rtol;
@@ -59,9 +103,12 @@ class petsc_solver<double>
// KSP Maximum number of iterations
PetscInt maxits;
- // Parallel solver
+ // Main parallel solver
KSP ksp;
+ // Temporal variable used for calculation of static members
+ size_t tmp;
+
// Solver used
char s_type[32];
@@ -78,380 +125,447 @@ class petsc_solver<double>
bool try_solve = false;
- // Residual behavior per iteration
- openfpm::vector<PetscReal> vres;
+ // It contain the solver benchmark results
+ openfpm::vector<solv_bench_info> bench;
- /*! \brief procedure to collect the preconditioned residue at each iteration
+ /*! \brief Calculate the residual error at time t for one method
*
+ * \param t time
+ * \param slv solver
+ *
+ * \return the residual
*
*/
- static PetscErrorCode monitor(KSP ksp,PetscInt it,PetscReal res,void* data)
+ static double calculate_it(double t, solv_bench_info & slv)
{
- openfpm::vector<PetscReal> * vres = static_cast<openfpm::vector<PetscReal> *>(data);
+ double s_int = slv.time / slv.res.size();
- vres->add(res);
+ // Calculate the discrete point in time
+ size_t pt = std::floor(t / s_int);
- return 0;
+ if (pt < slv.res.size())
+ return slv.res.get(pt).err_inf;
+
+ return slv.res.last().err_inf;
}
- /*! \brief try to solve the system with block jacobi pre-conditioner
- *
+ /*! \brief Here we write the benchmark report
*
*
*/
- void try_solve_complex_bj(Mat & A_, const Vec & b_, Vec & x_)
+ void write_bench_report()
{
- PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,5));
- PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
+ if (create_vcluster().getProcessUnitID() != 0)
+ return;
- solve_complex(A_,b_,x_);
- }
+ openfpm::vector<std::string> x;
+ openfpm::vector<openfpm::vector<double>> y;
+ openfpm::vector<std::string> yn;
+ openfpm::vector<double> xd;
- /*! \brief Try to solve the system using all the Krylov solver with simple preconditioners
- *
- *
- *
- */
- void try_solve_simple(Mat & A_, const Vec & b_, Vec & x_)
- {
- PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,3000));
- for (size_t i = 0 ; i < solvs.size() ; i++)
- {
- setSolver(solvs.get(i).c_str());
+ for (size_t i = 0 ; i < bench.size() ; i++)
+ x.add(bench.get(i).smethod);
- // Disable convergence check
- PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
+ // Each colum can have multiple data set (in this case 4 dataset)
+ // Each dataset can have a name
+ yn.add("Norm infinity");
+ yn.add("Norm average");
+ yn.add("Number of iterations");
- solve_simple(A_,b_,x_);
- }
- }
+ // Each colums can have multiple data-set
- /*! \brief solve complex use Krylov solver in combination with Block-Jacobi + Nested
- *
- *
- * Solve the system using block-jacobi preconditioner + nested solver for each block
- *
- */
- void solve_complex(Mat & A_, const Vec & b_, Vec & x_)
- {
- // We get the size of the Matrix A
- PetscInt row;
- PetscInt col;
- PetscInt row_loc;
- PetscInt col_loc;
- PetscInt * blks;
- PetscInt nlocal;
- PetscInt first;
- KSP *subksp;
- PC subpc;
+ for (size_t i = 0 ; i < bench.size() ; i++)
+ y.add({bench.get(i).err.err_inf,bench.get(i).err.err_norm,(double)bench.get(i).err.its});
- PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
- PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+ // Google charts options
+ GCoptions options;
- // Set the preconditioner to Block-jacobi
- PC pc;
- KSPGetPC(ksp,&pc);
- PCSetType(pc,PCBJACOBI);
- PETSC_SAFE_CALL(KSPSetType(ksp,KSPGMRES));
+ options.title = std::string("Residual error");
+ options.yAxis = std::string("Error");
+ options.xAxis = std::string("Method");
+ options.stype = std::string("bars");
+ options.more = std::string("vAxis: {scaleType: 'log'}");
- PetscInt m = 1;
+ GoogleChart cg;
- PetscMalloc1(m,&blks);
- for (size_t i = 0 ; i < m ; i++) blks[i] = row_loc;
+ std::stringstream ss;
- PCBJacobiSetLocalBlocks(pc,m,blks);
- PetscFree(blks);
+ cg.addHTML("<h2>Robustness</h2>");
+ cg.AddHistGraph(x,y,yn,options);
+ cg.addHTML("<h2>Speed</h2>");
- KSPSetUp(ksp);
- PCSetUp(pc);
+ y.clear();
+ yn.clear();
- PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);
+ yn.add("Time in millseconds");
+ options.title = std::string("Time");
- for (size_t i = 0; i < nlocal; i++)
- {
- KSPGetPC(subksp[i],&subpc);
+ for (size_t i = 0 ; i < bench.size() ; i++)
+ y.add({bench.get(i).time});
- PCSetType(subpc,PCLU);
+ cg.AddHistGraph(x,y,yn,options);
-// PCSetType(subpc,PCJACOBI);
- KSPSetType(subksp[i],KSPPREONLY);
-// PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedDefault,NULL,NULL));
-// KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
+ // Convergence in time
-/* if (!rank)
- {
- if (i%2)
- {
- PCSetType(subpc,PCILU);
- }
- else
- {
- PCSetType(subpc,PCNONE);
- KSPSetType(subksp[i],KSPBCGS);
- KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
- }
- }
- else
- {
- PCSetType(subpc,PCJACOBI);
- KSPSetType(subksp[i],KSPGMRES);
- KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
- }*/
- }
+ x.clear();
+ y.clear();
+ yn.clear();
+ xd.clear();
+ // Get the maximum in time across all the solvers
- KSPSolve(ksp,b_,x_);
+ double max_time = 0.0;
- auto & v_cl = create_vcluster();
- if (try_solve == true)
+ for (size_t i = 0 ; i < bench.size() ; i++)
{
- // calculate error statistic about the solution
- solError err = statSolutionError(A_,b_,x_);
-
- if (v_cl.getProcessUnitID() == 0)
- {
- std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
- std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
- }
+ if (bench.get(i).time > max_time) {max_time = bench.get(i).time;}
+ yn.add(bench.get(i).smethod);
}
- }
-
- void solve_ASM(Mat & A_, const Vec & b_, Vec & x_)
- {
-
-#if 0
- PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
+ size_t n_int = 300;
- // We set the size of x according to the Matrix A
- PetscInt row;
- PetscInt col;
- PetscInt row_loc;
- PetscInt col_loc;
+ // calculate dt
- PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
- PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+ double dt = max_time / n_int;
- PC pc;
+ // Add
- // We set the Matrix operators
- PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
+ // For each solver we have a convergence plot
- // We set the pre-conditioner
- PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
- PETSC_SAFE_CALL(PCSetType(pc,PCASM));
+ for (double t = dt ; t <= max_time + 0.05 * max_time ; t += dt)
+ {
+ y.add();
+ xd.add(t);
+ for (size_t j = 0 ; j < bench.size() ; j++)
+ y.last().add(calculate_it(t,bench.get(j)));
+ }
- PCASMSetOverlap(pc,1);
+ std::stringstream ss2;
+ ss2 << "<h2>Convergence with time</h2><br>" << std::endl;
- KSP *subksp; /* array of KSP contexts for local subblocks */
- PetscInt nlocal,first; /* number of local subblocks, first local subblock */
- PC subpc; /* PC context for subblock */
- PetscBool isasm;
+ options.title = std::string("Residual norm infinity");
+ options.yAxis = std::string("residual");
+ options.xAxis = std::string("Time in milliseconds");
+ options.lineWidth = 1.0;
+ options.more = std::string("vAxis: {scaleType: 'log'},hAxis: {scaleType: 'log'}");
- /*
- Set runtime options
- */
- KSPSetFromOptions(ksp);
+ cg.addHTML(ss2.str());
+ cg.AddLinesGraph(xd,y,yn,options);
- /*
- Flag an error if PCTYPE is changed from the runtime options
- */
- PetscObjectTypeCompare((PetscObject)pc,PCASM,&isasm);
- if (!isasm) SETERRQ(PETSC_COMM_WORLD,1,"Cannot Change the PCTYPE when manually changing the subdomain solver settings");
+ ss << "<h2>Solvers Tested</h2><br><br>" << std::endl;
+ for (size_t i = 0 ; i < bench.size() ; i++)
+ ss << bench.get(i).smethod << " = " << bench.get(i).method << "<br>" << std::endl;
- /*
- Call KSPSetUp() to set the block Jacobi data structures (including
- creation of an internal KSP context for each block).
+ cg.addHTML(ss.str());
- Note: KSPSetUp() MUST be called before PCASMGetSubKSP().
- */
- KSPSetUp(ksp);
+ cg.write("gc_solver.html");
+ }
- /*
- Extract the array of KSP contexts for the local blocks
- */
- PCASMGetSubKSP(pc,&nlocal,&first,&subksp);
+ /*! \brief Print a progress bar on standard out
+ *
+ *
+ */
+ void print_progress_bar()
+ {
+ auto & v_cl = create_vcluster();
- /*
- Loop over the local blocks, setting various KSP options
- for each block.
- */
- for (i=0; i<nlocal; i++)
- {
- KSPGetPC(subksp[i],&subpc);
- PCSetType(subpc,PCILU);
- KSPSetType(subksp[i],KSPGMRES);
- KSPSetTolerances(subksp[i],1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
- }
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout << "-----------------------25-----------------------50-----------------------75----------------------100" << std::endl;
+ }
- // if we are on on best solve set-up a monitor function
+ /*! \brief procedure to collect the absolute residue at each iteration
+ *
+ * \param ksp Solver
+ * \param it Iteration number
+ * \param res resudual
+ * \param custom pointer to data
+ *
+ */
+ static PetscErrorCode monitor(KSP ksp,PetscInt it,PetscReal res,void* data)
+ {
+ petsc_solver<double> * pts = static_cast<petsc_solver<double> *>(data);
- if (try_solve == true)
- {
- // for bench-mark we are interested in non-preconditioned norm
- PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
+ Mat A;
+ Mat P;
+ Vec x;
+ Vec b;
- // Disable convergence check
- PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
- }
-
- KSPGMRESSetRestart(ksp,200);
+ PETSC_SAFE_CALL(KSPGetOperators(ksp,&A,&P));
+ PETSC_SAFE_CALL(KSPGetSolution(ksp,&x));
+ PETSC_SAFE_CALL(KSPGetRhs(ksp,&b));
- // Solve the system
- PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
+ solError err = statSolutionError(A,b,x,ksp);
+ itError erri(err);
+ erri.it = it;
- auto & v_cl = create_vcluster();
-// if (try_solve == true)
-// {
- // calculate error statistic about the solution
- solError err = statSolutionError(A_,b_,x_);
+ // Add the error per iteration
+ pts->bench.last().res.add(erri);
- if (v_cl.getProcessUnitID() == 0)
- {
- std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
- std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
- }
-// }
+ pts->progress(it);
-#endif
+ return 0;
}
- void solve_AMG(Mat & A_, const Vec & b_, Vec & x_)
+ /*! \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_SAFE_CALL(KSPSetType(ksp,s_type));
- PETSC_SAFE_CALL(KSPSetType(ksp,KSPRICHARDSON));
-
- // -pc_hypre_boomeramg_tol <tol>
-
- // We set the size of x according to the Matrix A
- PetscInt row;
- PetscInt col;
- PetscInt row_loc;
- PetscInt col_loc;
+ petsc_solver<double> * pts = (petsc_solver *)data;
- PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
- PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+ pts->progress(it);
- PC pc;
+ return 0;
+ }
- // We set the Matrix operators
- PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
+ /*! \brief This function print an "*" showing the progress of the solvers
+ *
+ *
+ */
+ void progress(PetscInt it)
+ {
+ PetscInt n_max_it;
- // We set the pre-conditioner
- PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
+ PETSC_SAFE_CALL(KSPGetTolerances(ksp,NULL,NULL,NULL,&n_max_it));
- // PETSC_SAFE_CALL(PCSetType(pc,PCJACOBI));
+ auto & v_cl = create_vcluster();
- PETSC_SAFE_CALL(PCSetType(pc,PCHYPRE));
- // PCGAMGSetNSmooths(pc,0);
- PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
- PCFactorSetShiftAmount(pc, PETSC_DECIDE);
- PCHYPRESetType(pc, "boomeramg");
- MatSetBlockSize(A_,4);
- PetscOptionsSetValue("-pc_hypre_boomeramg_print_statistics","2");
- PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter","1000");
- PetscOptionsSetValue("-pc_hypre_boomeramg_nodal_coarsen","true");
- PetscOptionsSetValue("-pc_hypre_boomeramg_relax_type_all","SOR/Jacobi");
- PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type","Falgout");
- PetscOptionsSetValue("-pc_hypre_boomeramg_cycle_type","W");
- PetscOptionsSetValue("-pc_hypre_boomeramg_max_levels","20");
- PetscOptionsSetValue("-pc_hypre_boomeramg_vec_interp_variant","0");
- KSPSetFromOptions(ksp);
+ if (std::floor(it * 100.0 / n_max_it) > tmp)
+ {
+ tmp = it * 100.0 / n_max_it;
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout << "*" << std::flush;
+ }
+ }
- // if we are on on best solve set-up a monitor function
+ /*! \brief Print statistic about the solution error and method used
+ *
+ * \param err structure that contain the solution errors
+ *
+ */
+ void print_stat(solError & err)
+ {
+ auto & v_cl = create_vcluster();
- if (try_solve == true)
- {
- // for bench-mark we are interested in non-preconditioned norm
- PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
+ if (v_cl.getProcessUnitID() == 0)
+ {
+ std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
+ std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
+ }
+ }
- // Disable convergence check
- PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
- }
+ /*! \brief It convert the KSP type into a human read-able string
+ *
+ *
+ */
+ std::string to_string_method(const std::string & solv)
+ {
+ if (solv == std::string(KSPBCGS))
+ {
+ return std::string("BiCGStab (Stabilized version of BiConjugate Gradient Squared)");
+ }
+ else if (solv == std::string(KSPIBCGS))
+ {
+ return std::string("IBiCGStab (Improved Stabilized version of BiConjugate Gradient Squared)");
+ }
+ else if (solv == std::string(KSPFBCGS))
+ {
+ return std::string("FBiCGStab (Flexible Stabilized version of BiConjugate Gradient Squared)");
+ }
+ else if (solv == std::string(KSPFBCGSR))
+ {
+ return std::string("FBiCGStab-R (Modified Flexible Stabilized version of BiConjugate Gradient Squared)");
+ }
+ else if (solv == std::string(KSPBCGSL))
+ {
+ return std::string("BiCGStab(L) (Stabilized version of BiConjugate Gradient Squared with L search directions)");
+ }
+ else if (solv == std::string(KSPGMRES))
+ {
+ return std::string("GMRES(M) (Generalized Minimal Residual method with restart)");
+ }
+ else if (solv == std::string(KSPFGMRES))
+ {
+ return std::string("FGMRES(M) (Flexible Generalized Minimal Residual with restart)");
+ }
+ else if (solv == std::string(KSPLGMRES))
+ {
+ return std::string("LGMRES(M) (Augment Generalized Minimal Residual method with restart)");
+ }
+ else if (solv == std::string(KSPPGMRES))
+ {
+ return std::string("PGMRES(M) (Pipelined Generalized Minimal Residual method)");
+ }
+ else if (solv == std::string(KSPGCR))
+ {
+ return std::string("GCR (Generalized Conjugate Residual)");
+ }
- // Solve the system
- PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
+ return std::string("UNKNOWN");
+ }
- auto & v_cl = create_vcluster();
- // if (try_solve == true)
- // {
- // calculate error statistic about the solution
- solError err = statSolutionError(A_,b_,x_);
+ /*! \brief Allocate a new benchmark slot for a method
+ *
+ * \param str Method name
+ *
+ */
+ void new_bench(const std::string & str)
+ {
+ // Add a new benchmark result
+ bench.add();
+ bench.last().method = to_string_method(str);
+ bench.last().smethod = std::string(str);
+ }
- if (v_cl.getProcessUnitID() == 0)
- {
- std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
- std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
- }
- // }
+ /*! \brief Copy the solution if better
+ *
+ * \param res Residual of the solution
+ * \param sol solution
+ * \param best_res residual of the best solution
+ * \param best_sol best solution
+ *
+ */
+ void copy_if_better(double res, Vec & sol, double & best_res, Vec & best_sol)
+ {
+ if (res < best_res)
+ {
+ VecCopy(sol,best_sol);
+ best_res = res;
+ }
}
- /*! \brief solve simple use a Krylov solver + Simple selected Parallel Pre-conditioner
+ /*! \brief Try to solve the system x=inv(A)*b using all the Krylov solvers with simple Jacobi pre-conditioner
+ *
+ * It try to solve the system using JACOBI pre-conditioner and all the Krylov solvers available at the end it write
+ * a report
+ *
+ * \param A_ Matrix
+ * \param b_ vector of coefficents
+ * \param x solution
*
*/
- void solve_Krylov_simple(Mat & A_, const Vec & b_, Vec & x_)
+ void try_solve_simple(Mat & A_, const Vec & b_, Vec & x_)
{
- PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
+ Vec best_sol;
+ PETSC_SAFE_CALL(VecDuplicate(x_,&best_sol));
- // We set the size of x according to the Matrix A
- PetscInt row;
- PetscInt col;
- PetscInt row_loc;
- PetscInt col_loc;
+ // Best residual
+ double best_res = std::numeric_limits<double>::max();
- PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
- PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
+ // Create a new VCluster
+ auto & v_cl = create_vcluster();
- PC pc;
+ destroyKSP();
- // We set the Matrix operators
- PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
+ // for each solver
+ for (size_t i = 0 ; i < solvs.size() ; i++)
+ {
+ initKSPForTest();
+ setSolver(solvs.get(i).c_str());
- // We set the pre-conditioner
- PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
+ // Setup for BCGSL, GMRES
+ if (solvs.get(i) == std::string(KSPBCGSL))
+ {
+ // we try from 2 to 6 as search direction
+ for (size_t j = 2 ; j < 6 ; j++)
+ {
+ new_bench(solvs.get(i));
- // PETSC_SAFE_CALL(PCSetType(pc,PCJACOBI));
-
- PETSC_SAFE_CALL(PCSetType(pc,PCHYPRE));
-// PCGAMGSetNSmooths(pc,0);
- PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
- PCFactorSetShiftAmount(pc, PETSC_DECIDE);
- PCHYPRESetType(pc, "boomeramg");
- MatSetBlockSize(A_,4);
- PetscOptionsSetValue("-pc_hypre_boomeramg_print_statistics","2");
- PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter","1000");
- PetscOptionsSetValue("-pc_hypre_boomeramg_nodal_coarsen","true");
- PetscOptionsSetValue("-pc_hypre_boomeramg_relax_type_all","SOR/Jacobi");
- PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type","Falgout");
- PetscOptionsSetValue("-pc_hypre_boomeramg_cycle_type","W");
- PetscOptionsSetValue("-pc_hypre_boomeramg_max_levels","10");
- KSPSetFromOptions(ksp);
- // if we are on on best solve set-up a monitor function
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout << "L = " << j << std::endl;
+ bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
+ searchDirections(j);
+ bench_solve_simple(A_,b_,x_,bench.last());
- if (try_solve == true)
- {
- // for bench-mark we are interested in non-preconditioned norm
- PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
+ copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
+ }
+ }
+ else if (solvs.get(i) == std::string(KSPGMRES) ||
+ solvs.get(i) == std::string(std::string(KSPFGMRES)) ||
+ solvs.get(i) == std::string(std::string(KSPLGMRES)) )
+ {
+ // we try from 2 to 6 as search direction
+ for (size_t j = 50 ; j < 300 ; j += 50)
+ {
+ // Add a new benchmark result
+ new_bench(solvs.get(i));
- // Disable convergence check
- PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
- }
+ if (v_cl.getProcessUnitID() == 0)
+ std::cout << "Restart = " << j << std::endl;
+ bench.last().smethod += std::string("(") + std::to_string(j) + std::string(")");
+ setRestart(j);
+ bench_solve_simple(A_,b_,x_,bench.last());
- // Solve the system
- PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
+ copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
+ }
+ }
+ else
+ {
+ // Add a new benchmark result
+ new_bench(solvs.get(i));
- auto & v_cl = create_vcluster();
-// if (try_solve == true)
-// {
- // calculate error statistic about the solution
- solError err = statSolutionError(A_,b_,x_);
+ bench_solve_simple(A_,b_,x_,bench.last());
- if (v_cl.getProcessUnitID() == 0)
- {
- std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
- std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
+ copy_if_better(bench.last().err.err_inf,x_,best_res,best_sol);
}
-// }
+
+ destroyKSP();
+ }
+
+ write_bench_report();
+
+ // Copy the best solution to x
+ PETSC_SAFE_CALL(VecCopy(best_sol,x_));
+ }
+
+ /*! \brief Benchmark solve simple solving x=inv(A)*b
+ *
+ * \param A_ Matrix A
+ * \param b_ vector b
+ * \param x_ solution x
+ * \param bench structure that store the benchmark information
+ *
+ */
+ void bench_solve_simple(Mat & A_, const Vec & b_, Vec & x_, solv_bench_info & bench)
+ {
+ // timer for benchmark
+ timer t;
+ t.start();
+ // Enable progress monitor
+ tmp = 0;
+ PETSC_SAFE_CALL(KSPMonitorCancel(ksp));
+ PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor_progress,this,NULL));
+ print_progress_bar();
+ solve_simple(A_,b_,x_);
+
+ // New line
+ if (create_vcluster().getProcessUnitID() == 0)
+ std::cout << std::endl;
+
+ t.stop();
+ bench.time = t.getwct() * 1000.0;
+
+ // calculate error statistic about the solution and print
+ solError err = statSolutionError(A_,b_,x_,ksp);
+ print_stat(err);
+
+ 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
@@ -483,30 +597,14 @@ class petsc_solver<double>
if (try_solve == true)
{
- // for bench-mark we are interested in non-preconditioned norm
- PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
-
// Disable convergence check
PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
}
- KSPGMRESSetRestart(ksp,300);
+ PETSC_SAFE_CALL(KSPSetFromOptions(ksp));
// Solve the system
PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
-
- auto & v_cl = create_vcluster();
-// if (try_solve == true)
-// {
- // calculate error statistic about the solution
- solError err = statSolutionError(A_,b_,x_);
-
- if (v_cl.getProcessUnitID() == 0)
- {
- std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
- std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
- }
-// }
}
/*! \brief Calculate statistic on the error solution
@@ -516,7 +614,7 @@ class petsc_solver<double>
* \brief x Solution
*
*/
- solError statSolutionError(Mat & A_, const Vec & b_, Vec & x_)
+ static solError statSolutionError(Mat & A_, const Vec & b_, Vec & x_, KSP ksp)
{
PetscScalar neg_one = -1.0;
@@ -548,36 +646,94 @@ class petsc_solver<double>
PETSC_SAFE_CALL(KSPGetIterationNumber(ksp,&its));
solError err;
- err.err_norm = norm;
+ err.err_norm = norm / col;
err.err_inf = norm_inf;
err.its = its;
return err;
}
+ /*! \brief initialize the KSP object
+ *
+ *
+ */
+ void initKSP()
+ {
+ PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
+ PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
+ }
+
+ /*! \brief initialize the KSP object for solver testing
+ *
+ *
+ */
+ void initKSPForTest()
+ {
+ PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
+ PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
+ PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,300));
+
+ // Disable convergence check
+ PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
+ }
+
+ /*! \brief Destroy the KSP object
+ *
+ *
+ */
+ void destroyKSP()
+ {
+ KSPDestroy(&ksp);
+ }
public:
petsc_solver()
{
- strcpy(s_type,KSPBCGSL);
- PETSC_SAFE_CALL(KSPCreate(PETSC_COMM_WORLD,&ksp));
- PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
+ initKSP();
// Add the solvers
-// solvs.add(std::string(KSPBCGS));
-// solvs.add(std::string(KSPIBCGS));
-// solvs.add(std::string(KSPFBCGS));
-// solvs.add(std::string(KSPFBCGSR));
-// solvs.add(std::string(KSPBCGSL));
+ solvs.add(std::string(KSPBCGS));
+// solvs.add(std::string(KSPIBCGS)); <--- Produce invalid argument
+ solvs.add(std::string(KSPFBCGS));
+// solvs.add(std::string(KSPFBCGSR)); <--- Nan production problems
+ solvs.add(std::string(KSPBCGSL));
solvs.add(std::string(KSPGMRES));
-// solvs.add(std::string(KSPFGMRES));
-// solvs.add(std::string(KSPLGMRES));
-// solvs.add(std::string(KSPPGMRES));
+ solvs.add(std::string(KSPFGMRES));
+ solvs.add(std::string(KSPLGMRES));
+// solvs.add(std::string(KSPPGMRES)); <--- Take forever
// solvs.add(std::string(KSPGCR));
}
+ /*! \brief Add a test solver
+ *
+ * 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))
+ *
+ */
+ void addTestSolver(std::string & solver)
+ {
+ solvs.add(solver);
+ }
+
+ /*! \brief Remove a test solver
+ *
+ * The try solve function use the most robust solvers in PETSC, if you want to remove
+ * a solver like use removeTestSolver(std::string(KSPIBCGS))
+ *
+ */
+ void removeTestSolver(const std::string & solver)
+ {
+ // search for the test solver and remove it
+
+ for (size_t i = 0 ; i < solvs.size() ; i++)
+ {
+ if (solvs.get(i) == solver)
+ return;
+ }
+ }
+
/*! \brief Set the solver to test all the solvers and generate a report
*
* In this mode the system will try different Solvers, Preconditioner and
@@ -624,7 +780,7 @@ public:
void setAbsTol(PetscReal abstol_)
{
abstol = abstol_;
- KSPSetTolerances(ksp,0.0,abstol,30000.0,1000);
+ KSPSetTolerances(ksp,rtol,abstol,dtol,maxits);
}
/*! \brief Set the divergence tolerance
@@ -646,12 +802,22 @@ public:
* to avoid this problem. Such method has a parameter L (2 by default) Bigger is L
* more the system will try to avoid the breakdown
*
- * \size_t l Increasing L should reduce the probability of failure of the solver because of break-down of the base
+ * \param l Increasing L should reduce the probability of failure of the solver because of break-down of the base
*
*/
void searchDirections(PetscInt l)
{
- KSPBCGSLSetEll(ksp, l);
+ PetscOptionsSetValue("-ksp_bcgsl_ell",std::to_string(l).c_str());
+ }
+
+ /*! \brief For GMRES based method, the number of Krylov directions to orthogonalize against
+ *
+ * \param n number of directions
+ *
+ */
+ void setRestart(PetscInt n)
+ {
+ PetscOptionsSetValue("-ksp_gmres_restart",std::to_string(n).c_str());
}
/*! \brief Here we invert the matrix and solve the system
@@ -681,14 +847,9 @@ public:
Vec & x_ = x.getVec();
if (try_solve == true)
- {
try_solve_simple(A_,b_,x_);
-// try_solve_complex_bj(A_,b_,x_);
- }
else
- {
solve_simple(A_,b_,x_);
- }
x.update();
return x;
diff --git a/src/unit_test_init_cleanup.hpp b/src/unit_test_init_cleanup.hpp
index ca73aa358938f05a668dac4c3b41ac919fe82378..5128804068f705666fba9e5a79b0cb258bd161cb 100644
--- a/src/unit_test_init_cleanup.hpp
+++ b/src/unit_test_init_cleanup.hpp
@@ -11,8 +11,15 @@
#include "VCluster.hpp"
struct ut_start {
- ut_start() { BOOST_TEST_MESSAGE("Initialize global VCluster"); openfpm_init(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv); }
- ~ut_start() { BOOST_TEST_MESSAGE("Delete global VClster"); openfpm_finalize(); }
+ ut_start() {
+ BOOST_TEST_MESSAGE("Initialize global VCluster");
+ openfpm_init(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
+ }
+
+ ~ut_start() {
+ BOOST_TEST_MESSAGE("Delete global VClster");
+ openfpm_finalize();
+ }
};
//____________________________________________________________________________//