/* * petsc_solver.hpp * * Created on: Apr 26, 2016 * Author: i-bird */ #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ #define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_HPP_ #ifdef HAVE_PETSC #include "Vector/Vector.hpp" #include "Eigen/UmfPackSupport" #include #include #include #include "Plot/GoogleChart.hpp" #include "Matrix/SparseMatrix.hpp" #include "Vector/Vector.hpp" #define UMFPACK_NONE 0 template class petsc_solver { public: template static Vector solve(const SparseMatrix & A, const Vector & b) { std::cerr << "Error Petsc only suppor double precision" << "/n"; } }; #define SOLVER_NOOPTION 0 #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 { // It contain statistic of the error of the calculated solution struct solError { // infinity norm of the error PetscReal err_inf; // L1 norm of the error PetscReal err_norm; // Number of iterations 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; } // Default constructor itError() {} }; // 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 res; }; // KSP Maximum number of iterations PetscInt maxits; // Main parallel solver KSP ksp; // Temporal variable used for calculation of static members size_t tmp; // The full set of solvers openfpm::vector solvs; bool try_solve = false; // It contain the solver benchmark results openfpm::vector bench; /*! \brief Calculate the residual error at time t for one method * * \param t time * \param slv solver * * \return the residual * */ static double calculate_it(double t, solv_bench_info & slv) { double s_int = slv.time / slv.res.size(); // Calculate the discrete point in time size_t pt = std::floor(t / s_int); if (pt < slv.res.size()) return slv.res.get(pt).err_inf; return slv.res.last().err_inf; } /*! \brief Here we write the benchmark report * * */ void write_bench_report() { if (create_vcluster().getProcessUnitID() != 0) return; openfpm::vector x; openfpm::vector> y; openfpm::vector yn; openfpm::vector xd; for (size_t i = 0 ; i < bench.size() ; i++) x.add(bench.get(i).smethod); // 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"); // Each colums can have multiple data-set 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}); // Google charts options GCoptions options; 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'}"); GoogleChart cg; std::stringstream ss; cg.addHTML("

## Speed

"); y.clear(); yn.clear(); yn.add("Time in millseconds"); options.title = std::string("Time"); for (size_t i = 0 ; i < bench.size() ; i++) y.add({bench.get(i).time}); cg.AddHistGraph(x,y,yn,options); // Convergence in time x.clear(); y.clear(); yn.clear(); xd.clear(); // Get the maximum in time across all the solvers double max_time = 0.0; for (size_t i = 0 ; i < bench.size() ; i++) { if (bench.get(i).time > max_time) {max_time = bench.get(i).time;} yn.add(bench.get(i).smethod); } size_t n_int = maxits; // calculate dt double dt = max_time / n_int; // Add // For each solver we have a convergence plot 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))); } std::stringstream ss2; ss2 << "

## Convergence with time

" << std::endl; 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'}"); cg.addHTML(ss2.str()); cg.AddLinesGraph(xd,y,yn,options); ss << "

## Solvers Tested

" << std::endl; for (size_t i = 0 ; i < bench.size() ; i++) ss << bench.get(i).smethod << " = " << bench.get(i).method << "