From e0860bf8f1f6be7c64276257ddb9b8373336701f Mon Sep 17 00:00:00 2001 From: steigemann Date: Mon, 24 Sep 2012 21:44:07 +0000 Subject: [PATCH] Interface to MatShell matrix-type in PETSc git-svn-id: https://svn.dealii.org/trunk@26689 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/base/quadrature_lib.h | 74 +---- deal.II/include/deal.II/base/table_handler.h | 19 -- .../include/deal.II/lac/petsc_precondition.h | 259 +----------------- deal.II/include/deal.II/lac/petsc_solver.h | 12 - deal.II/source/base/quadrature_lib.cc | 87 ------ deal.II/source/base/table_handler.cc | 130 --------- deal.II/source/distributed/tria.cc | 74 +---- deal.II/source/lac/petsc_precondition.cc | 217 --------------- deal.II/source/lac/petsc_solver.cc | 138 ---------- 9 files changed, 6 insertions(+), 1004 deletions(-) diff --git a/deal.II/include/deal.II/base/quadrature_lib.h b/deal.II/include/deal.II/base/quadrature_lib.h index 678bb73933..a73a6d2f26 100644 --- a/deal.II/include/deal.II/base/quadrature_lib.h +++ b/deal.II/include/deal.II/base/quadrature_lib.h @@ -172,7 +172,7 @@ class QTrapez : public Quadrature /** * Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5. - * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102 + * See Stoer: Einführung in die Numerische Mathematik I, p. 102 */ template class QMilne : public Quadrature @@ -184,7 +184,7 @@ class QMilne : public Quadrature /** * Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7. - * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102 + * See Stoer: Einführung in die Numerische Mathematik I, p. 102 */ template class QWeddle : public Quadrature @@ -430,75 +430,6 @@ class QGaussOneOverR : public Quadrature -/** - * Gauss Quadrature Formula with $1/R^{3/2}$ weighting function. This formula - * can be used to to integrate $1/R^{3/2} \ f(x)$ on the reference - * element $[0,1]^2$, where $f$ is a smooth function without - * singularities, and $R$ is the distance from the point $x$ to the vertex - * $\xi$, given at construction time by specifying its index. Notice that - * this distance is evaluated in the reference element. - * - * This quadrature formula is a specialization of QGaussOneOverR. - * We apply a second transformation $R = t^2$ to cancel a singularity - * of order $R^{1/2}$: - * \f[ - * \int_0^1 \int_0^1 \frac{1}{R^{1/2}} f(x,y) \biggr|_{x < y} dxdy = - * \int_0^1 \int_0^{r(\pi/4 v)} \frac{1}{R^{1/2}} f(R \cos(\pi/4 v), R \sin(\pi/4 v)) \frac{\pi}{4} R dr dv = - * 2 \int_0^1 \int_0^{\sqrt{r(\pi/4 v)}} f(t^2 \cos(\pi/4 v), t^2 \sin(\pi/4 v)) \frac{\pi}{4} t^2 dt dv = - * 2 \int_0^1 \int_0^1 f(t^2, t^2 \tan(\pi/4 v)) \frac{\pi}{4} \frac{t^2}{(\cos(\pi/4 v))^{3/2}} dt dv - * \f] - * - * Upon construction it is possible to specify wether we want the - * singularity removed, or not. In other words, this quadrature can be - * used to integrate $g(x) = 1/R^{3/2}\ f(x)$, or simply $f(x)$, with the $1/R^{3/2}$ - * factor already included in the quadrature weights. - */ -template -class QGaussOneOverRThreeHalfs : public Quadrature -{ - public: - /** - * The constructor takes three arguments: the order of the Gauss - * formula, the index of the vertex where the singularity is - * located, and whether we include the weighting singular function - * inside the quadrature, or we leave it in the user function to - * be integrated. Notice that this constructor only works for the - * vertices of the quadrilateral. - * - * Traditionally, quadrature formulas include their weighting - * function, and the last argument is set to false by - * default. There are cases, however, where this is undesirable - * (for example when you only know that your singularity has the - * same order of 1/R, but cannot be written exactly in this - * way). - * - * In other words, you can use this function in either of - * the following way, obtaining the same result: - * - * @code - * QGaussOneOverRThreeHalfs singular_quad(order, vertex_id, false); - * // This will produce the integral of f(x)/R^{3/2} - * for(unsigned int i=0; ivertex(vertex_id)).norm(); - * integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R^{3/2}; - * } - * @endcode - */ - QGaussOneOverRThreeHalfs(const unsigned int n, - const unsigned int vertex_index, - const bool factor_out_singular_weight=false); -}; - - - /*@}*/ /* -------------- declaration of explicit specializations ------------- */ @@ -529,7 +460,6 @@ template <> QWeddle<1>::QWeddle (); template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert); template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag); template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag); -template <> QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs (const unsigned int n, const unsigned int index, const bool flag); diff --git a/deal.II/include/deal.II/base/table_handler.h b/deal.II/include/deal.II/base/table_handler.h index 8bcceffebf..0c6366717e 100644 --- a/deal.II/include/deal.II/base/table_handler.h +++ b/deal.II/include/deal.II/base/table_handler.h @@ -505,23 +505,6 @@ class TableHandler */ void write_tex (std::ostream &file, const bool with_header=true) const; - /** - * Write table as a tex file - * and use the booktabs package. - * Values are written in - * math mode: $value$ - * If with_header is set to false - * (it is true by default), then - * no "\documentclass{...}", - * "\begin{document}" and - * "\end{document}" are used. In - * this way the file can be - * included into an existing tex - * file using a command like - * "\input{table_file}". - */ - void write_results (std::ostream &file, const bool with_header=true) const; - /** * Read or write the data of this * object to or from a stream for @@ -888,8 +871,6 @@ namespace internal } } - - template void TableHandler::add_value (const std::string &key, const T value) diff --git a/deal.II/include/deal.II/lac/petsc_precondition.h b/deal.II/include/deal.II/lac/petsc_precondition.h index 43730abcf1..af3e05acab 100644 --- a/deal.II/include/deal.II/lac/petsc_precondition.h +++ b/deal.II/include/deal.II/lac/petsc_precondition.h @@ -804,107 +804,6 @@ namespace PETScWrappers bool output_details; }; - /** - * Empty Constructor. You need to call - * initialize() before using this - * object. - */ - PreconditionBoomerAMG (); - - /** - * Constructor. Take the matrix which - * is used to form the preconditioner, - * and additional flags if there are - * any. - */ - PreconditionBoomerAMG (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - /** - * Initializes the preconditioner - * object and calculate all data that - * is necessary for applying it in a - * solver. This function is - * automatically called when calling - * the constructor with the same - * arguments and is only used if you - * create the preconditioner without - * arguments. - */ - void initialize (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - protected: - /** - * Store a copy of the flags for this - * particular preconditioner. - */ - AdditionalData additional_data; - }; - - - -/** - * A class that implements the interface to use the ParaSails sparse - * approximate inverse preconditioner from the HYPRE suite. Note that - * PETSc has to be configured with HYPRE (e.g. with --download-hypre=1). - * - * @ingroup PETScWrappers - * @author Martin Steigemann, 2011 - */ - class PreconditionParaSails : public PreconditionerBase - { - public: - /** - * Standardized data struct to - * pipe additional flags to the - * preconditioner. - */ - struct AdditionalData - { - /** - * Constructor. - */ - AdditionalData ( - const unsigned int symmetric = 0, - const unsigned int n_levels = 1, - const double threshold = 0.1, - const double filter = 0.05, - const double load_bal = 0., - const bool output_details = false - ); - - /** - * This parameter has the following meanings, - * to indicate the symmetry and definiteness - * of the problem, and to specify the type - * of the preconditioner to construct: - *
    - *
  • @p 0: nonsymmetric and/or indefinite problem, and nonsymmetric preconditioner - *
  • @p 1: SPD problem, and SPD (factored) preconditioner - *
  • @p 2: nonsymmetric, definite problem, and SPD (factored) preconditioner - *
- * - */ - unsigned int symmetric; - - unsigned int n_levels; - - double threshold; - - double filter; - - double load_bal; - - /** - * Setting this flag to true - * produces output from HYPRE, - * when the preconditioner - * is constructed. - */ - bool output_details; - }; - /** @@ -912,7 +811,7 @@ namespace PETScWrappers * initialize() before using this * object. */ - PreconditionParaSails (); + PreconditionBoomerAMG (); /** * Constructor. Take the matrix which @@ -920,7 +819,7 @@ namespace PETScWrappers * and additional flags if there are * any. */ - PreconditionParaSails (const MatrixBase &matrix, + PreconditionBoomerAMG (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); /** @@ -944,160 +843,6 @@ namespace PETScWrappers */ AdditionalData additional_data; }; - - - -/** - * A class that implements the interface to use the scalable implementation - * of the Parallel ILU algorithm in the Euclid library from the HYPRE suite. - * Note that PETSc has to be configured with HYPRE (e.g. with --download-hypre=1). - * - * @ingroup PETScWrappers - * @author Martin Steigemann, 2011 - */ - class PreconditionEuclid : public PreconditionerBase - { - public: - /** - * Standardized data struct to - * pipe additional flags to the - * preconditioner. - */ - struct AdditionalData - { - /** - * Constructor. - */ - AdditionalData ( - const unsigned int level = 1, - const bool use_block_jacobi = false, - const bool output_details = false - ); - - /** - * Factorization level for ILU(k). Default is 1. - * For 2D convection-diffusion and similar problems, - * fastest solution time is typically obtained with - * levels 4 through 8. For 3D problems, fastest solution - * time is typically obtained with level 1. - */ - unsigned int level; - - /** - * Use Block Jacobi ILU preconditioning - * instead of PILU. Default is false. - * If subdomains contain relatively few nodes - * (less than 1000), or the problem is not - * well partitioned, Block Jacobi ILU - * may give faster solution time then PILU. - */ - unsigned int use_block_jacobi; - - /** - * Setting this flag to true - * produces debug output from - * HYPRE, when the preconditioner - * is constructed. - */ - bool output_details; - }; - - - - /** - * Empty Constructor. You need to call - * initialize() before using this - * object. - */ - PreconditionEuclid (); - - /** - * Constructor. Take the matrix which - * is used to form the preconditioner, - * and additional flags if there are - * any. - */ - PreconditionEuclid (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - /** - * Initializes the preconditioner - * object and calculate all data that - * is necessary for applying it in a - * solver. This function is - * automatically called when calling - * the constructor with the same - * arguments and is only used if you - * create the preconditioner without - * arguments. - */ - void initialize (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - protected: - /** - * Store a copy of the flags for this - * particular preconditioner. - */ - AdditionalData additional_data; - }; - - - -/** - * A class that implements a non-preconditioned Krylov method. - * - * @ingroup PETScWrappers - * @author Martin Steigemann, 2011 - */ - class PreconditionNone : public PreconditionerBase - { - public: - /** - * Standardized data struct to - * pipe additional flags to the - * preconditioner. - */ - struct AdditionalData - {}; - - /** - * Empty Constructor. You need to call - * initialize() before using this - * object. - */ - PreconditionNone (); - - /** - * Constructor. Take the matrix which - * is used to form the preconditioner, - * and additional flags if there are - * any. - */ - PreconditionNone (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - /** - * Initializes the preconditioner - * object and calculate all data that - * is necessary for applying it in a - * solver. This function is - * automatically called when calling - * the constructor with the same - * arguments and is only used if you - * create the preconditioner without - * arguments. - */ - void initialize (const MatrixBase &matrix, - const AdditionalData &additional_data = AdditionalData()); - - protected: - /** - * Store a copy of the flags for this - * particular preconditioner. - */ - AdditionalData additional_data; - }; } diff --git a/deal.II/include/deal.II/lac/petsc_solver.h b/deal.II/include/deal.II/lac/petsc_solver.h index 0f3b3b5ab4..984ea7b6b5 100644 --- a/deal.II/include/deal.II/lac/petsc_solver.h +++ b/deal.II/include/deal.II/lac/petsc_solver.h @@ -138,14 +138,6 @@ namespace PETScWrappers const PreconditionerBase &preconditioner); - - void - solve (const MatrixBase &A, - VectorBase &x, - const VectorBase &b, - const PreconditionerBase &preconditioner, - const std::vector &nullspace); - /** * Resets the contained preconditioner * and solver object. See class @@ -177,10 +169,6 @@ namespace PETScWrappers << "An error with error number " << arg1 << " occurred while calling a PETSc function"); - DeclException1 (ExcPETScSolverError, - char*, - << "PETSc solver failed: " << arg1); - protected: /** diff --git a/deal.II/source/base/quadrature_lib.cc b/deal.II/source/base/quadrature_lib.cc index 498f45f66b..a38e566433 100644 --- a/deal.II/source/base/quadrature_lib.cc +++ b/deal.II/source/base/quadrature_lib.cc @@ -925,93 +925,6 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, } - -template<> -QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs(const unsigned int n, - const unsigned int vertex_index, - const bool factor_out_singularity) : - Quadrature<2>(2*n*n) -{ - // This version of the constructor - // works only for the 4 - // vertices. If you need a more - // general one, you should use the - // one with the Point<2> in the - // constructor. - Assert(vertex_index <4, ExcIndexRange(vertex_index, 0, 4)); - - // Start with the gauss quadrature - // formula on the (u,v) reference - // element. - QGauss<2> gauss(n); - - Assert(gauss.size() == n*n, ExcInternalError()); - Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4)); - - // We create only the first one. All other pieces are rotation of - // this one. - // In this case the transformation is - // - // (x,y) = (u*u, u*u tan(pi/4 v)) - // - // with Jacobian - // - // J = 2 pi/4 R*R sqrt(cos(pi/4 v)) - // - // And we get rid of R to take into account the singularity, - // unless specified differently in the constructor. - std::vector > &ps = this->quadrature_points; - std::vector &ws = this->weights; - double pi4 = numbers::PI/4; - - for(unsigned int q=0; q &gp = gauss.point(q); - ps[q][0] = gp[0]*gp[0]; - ps[q][1] = gp[0]*gp[0]*std::tan(pi4 *gp[1]); - ws[q] = 2.*gauss.weight(q)*pi4/std::sqrt(std::cos(pi4 *gp[1])); - if(factor_out_singularity) { - const double abs_value = (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm(); - ws[q] *= abs_value; - ws[q] *= std::sqrt(abs_value); - } - // The other half of the quadrilateral is symmetric with - // respect to xy plane. - ws[gauss.size()+q] = ws[q]; - ps[gauss.size()+q][0] = ps[q][1]; - ps[gauss.size()+q][1] = ps[q][0]; - } - - // Now we distribute these vertices in the correct manner - double theta = 0; - switch(vertex_index) { - case 0: - theta = 0; - break; - case 1: - theta = numbers::PI/2; - break; - case 2: - theta = -numbers::PI/2; - break; - case 3: - theta = numbers::PI; - break; - } - - double R00 = std::cos(theta), R01 = -std::sin(theta); - double R10 = std::sin(theta), R11 = std::cos(theta); - - if(vertex_index != 0) - for(unsigned int q=0; q sel_columns; - get_selected_columns(sel_columns); - - // write the column formats - for (unsigned int j=0; j >::const_iterator - super_iter=supercolumns.find(key); - - if (super_iter!=supercolumns.end()) - { - const unsigned int n_subcolumns=super_iter->second.size(); - for (unsigned int k=0; k::const_iterator - col_iter=columns.find(super_iter->second[k]); - Assert(col_iter!=columns.end(), ExcInternalError()); - - out << col_iter->second.tex_format; - } - } - else - { - // avoid `columns[key]'; - const std::map::const_iterator - col_iter=columns.find(key); - Assert(col_iter!=columns.end(), ExcInternalError()); - out << col_iter->second.tex_format; - } - } - out << "@{}} \\toprule" << std::endl; - - // write the caption line of the table - for (unsigned int j=0; j >::const_iterator - super_iter=supercolumns.find(key); - - if (super_iter!=supercolumns.end()) - { - const unsigned int n_subcolumns=super_iter->second.size(); - // avoid use of `tex_supercaptions[key]' - std::map::const_iterator - tex_super_cap_iter=tex_supercaptions.find(key); - out << std::endl << "\\multicolumn{" << n_subcolumns << "}{c}{" - << tex_super_cap_iter->second << "}"; - } - else - { - // col_iter->second=columns[col]; - const std::map::const_iterator - col_iter=columns.find(key); - Assert(col_iter!=columns.end(), ExcInternalError()); - out << col_iter->second.tex_caption; - } - if (j::const_iterator - col_iter=columns.find(key); - Assert(col_iter!=columns.end(), ExcInternalError()); - - const Column &column=col_iter->second; - - out << std::setprecision(column.precision); - - if (col_iter->second.scientific) - out.setf(std::ios::scientific, std::ios::floatfield); - else - out.setf(std::ios::fixed, std::ios::floatfield); - - if (math_mode) - out << "$"; - - out << column.entries[i].value; - - if (math_mode) - out << "$"; - - if (j void Triangulation:: - copy_triangulation (const dealii::Triangulation &old_tria) + copy_triangulation (const dealii::Triangulation &) { - clear(); - - try - { - dealii::Triangulation:: - copy_triangulation (old_tria); - } - catch (const typename dealii::Triangulation::DistortedCellList &) - { - // the underlying - // triangulation should not - // be checking for - // distorted cells - AssertThrow (false, ExcInternalError()); - } - - // note that now we have some content in - // the p4est objects and call the - // functions that do the actual work - // (which are dimension dependent, so - // separate) - triangulation_has_content = true; - - Assert (old_tria.n_levels() == 1, - ExcMessage ("Parallel distributed triangulations can only be copied, " - "if they are not refined!")); - - if (dynamic_cast *>(&old_tria) != 0) - { - Assert (!(dynamic_cast&> - (old_tria).refinement_in_progress), - ExcMessage ("Parallel distributed triangulations can only " - "be copied, if no refinement is in progress!")); - - coarse_cell_to_p4est_tree_permutation = - dynamic_cast&> - (old_tria).coarse_cell_to_p4est_tree_permutation; - - p4est_tree_to_coarse_cell_permutation = - dynamic_cast&> - (old_tria).p4est_tree_to_coarse_cell_permutation; - - attached_data_size = - dynamic_cast&> - (old_tria).attached_data_size; - - n_attached_datas = - dynamic_cast&> - (old_tria).n_attached_datas; - } - else - { - setup_coarse_cell_to_p4est_tree_permutation (); - }; - - copy_new_triangulation_to_p4est (dealii::internal::int2type()); - - try - { - copy_local_forest_to_triangulation (); - } - catch (const typename Triangulation::DistortedCellList &) - { - // the underlying - // triangulation should not - // be checking for - // distorted cells - AssertThrow (false, ExcInternalError()); - } - - update_number_cache (); + Assert (false, ExcNotImplemented()); } diff --git a/deal.II/source/lac/petsc_precondition.cc b/deal.II/source/lac/petsc_precondition.cc index 370d5c99b5..2df82399d2 100644 --- a/deal.II/source/lac/petsc_precondition.cc +++ b/deal.II/source/lac/petsc_precondition.cc @@ -483,223 +483,6 @@ namespace PETScWrappers } -/* ----------------- PreconditionParaSails -------------------- */ - - PreconditionParaSails::AdditionalData:: - AdditionalData(const unsigned int symmetric, - const unsigned int n_levels, - const double threshold, - const double filter, - const double load_bal, - const bool output_details) - : - symmetric(symmetric), - n_levels(n_levels), - threshold(threshold), - filter(filter), - load_bal(load_bal), - output_details(output_details) - {} - - - PreconditionParaSails::PreconditionParaSails () - {} - - - PreconditionParaSails::PreconditionParaSails (const MatrixBase &matrix, - const AdditionalData &additional_data) - { - initialize(matrix, additional_data); - } - - - void - PreconditionParaSails::initialize (const MatrixBase &matrix_, - const AdditionalData &additional_data_) - { - matrix = static_cast(matrix_); - additional_data = additional_data_; - -#ifdef PETSC_HAVE_HYPRE - create_pc(); - - int ierr; - ierr = PCSetType (pc, const_cast(PCHYPRE)); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCHYPRESetType(pc, "parasails"); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - if (additional_data.output_details) - PetscOptionsSetValue("-pc_hypre_parasails_logging","1"); - - Assert ((additional_data.symmetric == 0 || - additional_data.symmetric == 1 || - additional_data.symmetric == 2), - ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!")); - - std::stringstream ssStream; - - switch (additional_data.symmetric) - { - case 0: - { - ssStream << "nonsymmetric"; - break; - } - - case 1: - { - ssStream << "SPD"; - break; - } - - case 2: - { - ssStream << "nonsymmetric,SPD"; - break; - } - - default: - Assert (false, - ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!")); - }; - - PetscOptionsSetValue("-pc_hypre_parasails_sym",ssStream.str().c_str()); - - PetscOptionsSetValue("-pc_hypre_parasails_nlevels", - Utilities::int_to_string( - additional_data.n_levels - ).c_str()); - - ssStream.str(""); // empty the stringstream - ssStream << additional_data.threshold; - PetscOptionsSetValue("-pc_hypre_parasails_thresh", ssStream.str().c_str()); - - ssStream.str(""); // empty the stringstream - ssStream << additional_data.filter; - PetscOptionsSetValue("-pc_hypre_parasails_filter", ssStream.str().c_str()); - - ssStream.str(""); // empty the stringstream - ssStream << additional_data.load_bal; - PetscOptionsSetValue("-pc_hypre_parasails_loadbal", ssStream.str().c_str()); - - ierr = PCSetFromOptions (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - -#else // PETSC_HAVE_HYPRE - (void)pc; - Assert (false, - ExcMessage ("Your PETSc installation does not include a copy of " - "the hypre package necessary for this preconditioner.")); -#endif - } - - -/* ----------------- PreconditionEuclid ----------------------- */ - - PreconditionEuclid::AdditionalData:: - AdditionalData(const unsigned int level, - const bool use_block_jacobi, - const bool output_details) - : - level(level), - use_block_jacobi(use_block_jacobi), - output_details(output_details) - {} - - - PreconditionEuclid::PreconditionEuclid () - {} - - - PreconditionEuclid::PreconditionEuclid (const MatrixBase &matrix, - const AdditionalData &additional_data) - { - initialize(matrix, additional_data); - } - - - void - PreconditionEuclid::initialize (const MatrixBase &matrix_, - const AdditionalData &additional_data_) - { - matrix = static_cast(matrix_); - additional_data = additional_data_; - -#ifdef PETSC_HAVE_HYPRE - create_pc(); - - int ierr; - ierr = PCSetType (pc, const_cast(PCHYPRE)); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCHYPRESetType(pc, "euclid"); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - PetscOptionsSetValue("-pc_hypre_euclid_levels", - Utilities::int_to_string( - additional_data.level - ).c_str()); - - if (additional_data.use_block_jacobi) - PetscOptionsSetValue("-pc_hypre_euclid_bj","1"); - - if (additional_data.output_details) - PetscOptionsSetValue("-pc_hypre_euclid_print_statistics","1"); - - ierr = PCSetFromOptions (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - -#else // PETSC_HAVE_HYPRE - (void)pc; - Assert (false, - ExcMessage ("Your PETSc installation does not include a copy of " - "the hypre package necessary for this preconditioner.")); -#endif - } - - -/* ----------------- PreconditionNone ------------------------- */ - - PreconditionNone::PreconditionNone () - {} - - - PreconditionNone::PreconditionNone (const MatrixBase &matrix, - const AdditionalData &additional_data) - { - initialize(matrix, additional_data); - } - - - void - PreconditionNone::initialize (const MatrixBase &matrix_, - const AdditionalData &additional_data_) - { - matrix = static_cast(matrix_); - additional_data = additional_data_; - - create_pc(); - - int ierr; - ierr = PCSetType (pc, const_cast(PCNONE)); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetFromOptions (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp (pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - } - - /* ----------------- PreconditionLU -------------------- */ PreconditionLU::AdditionalData:: diff --git a/deal.II/source/lac/petsc_solver.cc b/deal.II/source/lac/petsc_solver.cc index 0cbb09d2bc..d01dae47f0 100644 --- a/deal.II/source/lac/petsc_solver.cc +++ b/deal.II/source/lac/petsc_solver.cc @@ -129,138 +129,6 @@ namespace PETScWrappers } - void - SolverBase::solve (const MatrixBase &A, - VectorBase &x, - const VectorBase &b, - const PreconditionerBase &preconditioner, - const std::vector &nullspace) - { - int ierr; - // first create a solver object if this - // is necessary - if (solver_data.get() == 0) - { - solver_data.reset (new SolverData()); - - ierr = KSPCreate (mpi_communicator, &solver_data->ksp); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - // set the matrices involved. the - // last argument is irrelevant here, - // since we use the solver only once - // anyway - ierr = KSPSetOperators (solver_data->ksp, A, preconditioner, - SAME_PRECONDITIONER); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - // let derived classes set the solver - // type, and the preconditioning - // object set the type of - // preconditioner - set_solver_type (solver_data->ksp); - - ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc()); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - PC prec; - ierr = KSPGetPC (solver_data->ksp, &prec); - ierr = PCFactorSetShiftType (prec, MAT_SHIFT_POSITIVE_DEFINITE); - - // then a convergence monitor - // function. that function simply - // checks with the solver_control - // object we have in this object for - // convergence -#if DEAL_II_PETSC_VERSION_LT(3,0,0) - KSPSetConvergenceTest (solver_data->ksp, &convergence_test, - reinterpret_cast(&solver_control)); -#else - KSPSetConvergenceTest (solver_data->ksp, &convergence_test, - reinterpret_cast(&solver_control), - PETSC_NULL); -#endif - - } - - // then do the real work: set up solver - // internal data and solve the - // system. - ierr = KSPSetUp (solver_data->ksp); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - unsigned int dim_nullsp = nullspace.size(); - - if (dim_nullsp > 0) - { - Vec * nullsp_basis; - - PetscMalloc (dim_nullsp*sizeof(Vec), &nullsp_basis); - - for (unsigned int i=0; iksp, nullsp); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - }; - - ierr = KSPSolve (solver_data->ksp, b, x); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - KSPConvergedReason reason; - ierr = KSPGetConvergedReason (solver_data->ksp, &reason); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - if (reason < 0) - { - std::ostringstream error_code; - - if (reason == KSP_DIVERGED_NULL) - error_code << "KSP_DIVERGED_NULL"; - else if (reason == KSP_DIVERGED_ITS) - error_code << "KSP_DIVERGED_ITS"; - else if (reason == KSP_DIVERGED_ITS) - error_code << "KSP_DIVERGED_DTOL"; - else if (reason == KSP_DIVERGED_BREAKDOWN) - error_code << "KSP_DIVERGED_BREAKDOWN"; - else if (reason == KSP_DIVERGED_BREAKDOWN_BICG) - error_code << "KSP_DIVERGED_BREAKDOWN_BICG"; - else if (reason == KSP_DIVERGED_NONSYMMETRIC) - error_code << "KSP_DIVERGED_NONSYMMETRIC"; - else if (reason == KSP_DIVERGED_INDEFINITE_PC) - error_code << "KSP_DIVERGED_INDEFINITE_PC"; - else if (reason == KSP_DIVERGED_NAN) - error_code << "KSP_DIVERGED_NAN"; - else if (reason == KSP_DIVERGED_INDEFINITE_MAT) - error_code << "KSP_DIVERGED_INDEFINITE_MAT"; - else - error_code << "Unknown Error"; - - AssertThrow (false, ExcPETScSolverError(error_code.str().c_str())); - }; - - // do not destroy solver object -// solver_data.reset (); - - // in case of failure: throw - // exception - if (solver_control.last_check() != SolverControl::success) - throw SolverControl::NoConvergence (solver_control.last_step(), - solver_control.last_value()); - // otherwise exit as normal - } - - void SolverBase::set_prefix(const std::string &prefix) { @@ -434,9 +302,6 @@ namespace PETScWrappers // honor the initial guess in the // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); - - KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(), - PETSC_DEFAULT, this->solver_control.max_steps()+1); } @@ -539,9 +404,6 @@ namespace PETScWrappers // honor the initial guess in the // solution vector. do so here as well: KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); - - KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(), - PETSC_DEFAULT, this->solver_control.max_steps()+1); } -- 2.39.5