From 80cd59fc170a549ce9eb77d93963047a4f6203f6 Mon Sep 17 00:00:00 2001 From: steigemann Date: Mon, 24 Sep 2012 20:43:44 +0000 Subject: [PATCH] //--------------------------------------------------------------------------- // $Id: petsc_matrix_free.cc 26043 2012-09-24 19:25:57Z steigemann $ // Version: $Name$ // // Copyright (C) 2012 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer // to the file deal.II/doc/license.html for the text and // further information on this license. // //--------------------------------------------------------------------------- #include #ifdef DEAL_II_USE_PETSC DEAL_II_NAMESPACE_OPEN namespace PETScWrappers { MatrixFree::MatrixFree () : communicator (PETSC_COMM_SELF) { const int m=0; do_reinit (m, m, m, m); } MatrixFree::MatrixFree (const MPI_Comm &communicator, const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns) : communicator (communicator) { do_reinit (m, n, local_rows, local_columns); } MatrixFree::MatrixFree (const MPI_Comm &communicator, const unsigned int m, const unsigned int n, const std::vector &local_rows_per_process, const std::vector &local_columns_per_process, const unsigned int this_process) : communicator (communicator) { Assert (local_rows_per_process.size() == local_columns_per_process.size(), ExcDimensionMismatch (local_rows_per_process.size(), local_columns_per_process.size())); Assert (this_process < local_rows_per_process.size(), ExcInternalError()); do_reinit (m, n, local_rows_per_process[this_process], local_columns_per_process[this_process]); } MatrixFree::MatrixFree (const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns) : communicator (MPI_COMM_WORLD) { do_reinit (m, n, local_rows, local_columns); } MatrixFree::MatrixFree (const unsigned int m, const unsigned int n, const std::vector &local_rows_per_process, const std::vector &local_columns_per_process, const unsigned int this_process) : communicator (MPI_COMM_WORLD) { Assert (local_rows_per_process.size() == local_columns_per_process.size(), ExcDimensionMismatch (local_rows_per_process.size(), local_columns_per_process.size())); Assert (this_process < local_rows_per_process.size(), ExcInternalError()); do_reinit (m, n, local_rows_per_process[this_process], local_columns_per_process[this_process]); } void MatrixFree::reinit (const MPI_Comm &communicator, const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns) { this->communicator = communicator; // destroy the matrix and // generate a new one #if DEAL_II_PETSC_VERSION_LT(3,2,0) int ierr = MatDestroy (matrix); #else int ierr = MatDestroy (&matrix); #endif AssertThrow (ierr == 0, ExcPETScError(ierr)); do_reinit (m, n, local_rows, local_columns); } void MatrixFree::reinit (const MPI_Comm &communicator, const unsigned int m, const unsigned int n, const std::vector &local_rows_per_process, const std::vector &local_columns_per_process, const unsigned int this_process) { Assert (local_rows_per_process.size() == local_columns_per_process.size(), ExcDimensionMismatch (local_rows_per_process.size(), local_columns_per_process.size())); Assert (this_process < local_rows_per_process.size(), ExcInternalError()); this->communicator = communicator; #if DEAL_II_PETSC_VERSION_LT(3,2,0) int ierr = MatDestroy (matrix); #else int ierr = MatDestroy (&matrix); #endif AssertThrow (ierr == 0, ExcPETScError(ierr)); do_reinit (m, n, local_rows_per_process[this_process], local_columns_per_process[this_process]); } void MatrixFree::reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns) { reinit (MPI_COMM_WORLD, m, n, local_rows, local_columns); } void MatrixFree::reinit (const unsigned int m, const unsigned int n, const std::vector &local_rows_per_process, const std::vector &local_columns_per_process, const unsigned int this_process) { reinit (MPI_COMM_WORLD, m, n, local_rows_per_process, local_columns_per_process, this_process); } void MatrixFree::clear () { #if DEAL_II_PETSC_VERSION_LT(3,2,0) int ierr = MatDestroy (matrix); #else int ierr = MatDestroy (&matrix); #endif AssertThrow (ierr == 0, ExcPETScError(ierr)); const int m=0; do_reinit (m, m, m, m); } void MatrixFree::vmult (Vec &dst, const Vec &src) const { //TODO: Translate the given PETSc Vec* vector into a deal.II // vector so we can call the vmult function with the usual // interface; then convert back. This could be much more // efficient, if the PETScWrappers::*::Vector classes // had a way to simply generate such a vector object from // a given PETSc Vec* object without allocating new memory // and without taking ownership of the Vec* VectorBase *x = 0; VectorBase *y = 0; // because we do not know, // if dst and src are sequential // or distributed vectors, // we ask for the vector-type // and reinit x and y with // dealii::PETScWrappers::*::Vector: const char *vec_type; int ierr = VecGetType (src, &vec_type); PetscInt local_size; ierr = VecGetLocalSize (src, &local_size); AssertThrow (ierr == 0, ExcPETScError(ierr)); if (strcmp(vec_type,"mpi") == 0) { PetscInt size; ierr = VecGetSize (src, &size); AssertThrow (ierr == 0, ExcPETScError(ierr)); x = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size); y = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size); } else if (strcmp(vec_type,"seq") == 0) { x = new PETScWrappers::Vector (local_size); y = new PETScWrappers::Vector (local_size); } else AssertThrow (false, ExcMessage("PETScWrappers::MPI::MatrixFree::do_matrix_vector_action: " "This only works for Petsc Vec Type = VECMPI | VECSEQ")); // copy src to x x->equ(1., PETScWrappers::VectorBase(src)); // and call vmult(x,y) which must // be reimplemented in derived classes vmult (*y, *x); y->compress(); // copy the result back to dst ierr = VecCopy (&(*(*y)), dst); AssertThrow (ierr == 0, ExcPETScError(ierr)); delete (x); delete (y); } int MatrixFree::matrix_free_mult (Mat A, Vec src, Vec dst) { // create a pointer to this MatrixFree // object and link the given matrix A // to the matrix-vector multiplication // of this MatrixFree object, MatrixFree *this_object; int ierr = MatShellGetContext (A, &this_object); AssertThrow (ierr == 0, ExcPETScError(ierr)); // call vmult of this object: this_object->vmult (dst, src); return (0); } void MatrixFree::do_reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, const unsigned int local_columns) { Assert (local_rows <= m, ExcDimensionMismatch (local_rows, m)); Assert (local_columns <= n, ExcDimensionMismatch (local_columns, n)); int ierr; // create a PETSc MatShell matrix-type // object of dimension m x n and local size // local_rows x local_columns ierr = MatCreateShell(communicator, local_rows, local_columns, m, n, (void*)this, &matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); // register the MatrixFree::matrix_free_mult function // as the matrix multiplication used by this matrix ierr = MatShellSetOperation (matrix, MATOP_MULT, (void(*)(void))&dealii::PETScWrappers::MatrixFree::matrix_free_mult); AssertThrow (ierr == 0, ExcPETScError(ierr)); ierr = MatSetFromOptions (matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); } } DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_USE_PETSC git-svn-id: https://svn.dealii.org/trunk@26688 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, 1004 insertions(+), 6 deletions(-) diff --git a/deal.II/include/deal.II/base/quadrature_lib.h b/deal.II/include/deal.II/base/quadrature_lib.h index a73a6d2f26..678bb73933 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,6 +430,75 @@ 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 ------------- */ @@ -460,6 +529,7 @@ 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 0c6366717e..8bcceffebf 100644 --- a/deal.II/include/deal.II/base/table_handler.h +++ b/deal.II/include/deal.II/base/table_handler.h @@ -505,6 +505,23 @@ 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 @@ -871,6 +888,8 @@ 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 af3e05acab..43730abcf1 100644 --- a/deal.II/include/deal.II/lac/petsc_precondition.h +++ b/deal.II/include/deal.II/lac/petsc_precondition.h @@ -804,8 +804,6 @@ namespace PETScWrappers bool output_details; }; - - /** * Empty Constructor. You need to call * initialize() before using this @@ -843,6 +841,263 @@ namespace PETScWrappers */ 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; + }; + + + + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionParaSails (); + + /** + * Constructor. Take the matrix which + * is used to form the preconditioner, + * and additional flags if there are + * any. + */ + PreconditionParaSails (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 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 984ea7b6b5..0f3b3b5ab4 100644 --- a/deal.II/include/deal.II/lac/petsc_solver.h +++ b/deal.II/include/deal.II/lac/petsc_solver.h @@ -138,6 +138,14 @@ 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 @@ -169,6 +177,10 @@ 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 a38e566433..498f45f66b 100644 --- a/deal.II/source/base/quadrature_lib.cc +++ b/deal.II/source/base/quadrature_lib.cc @@ -925,6 +925,93 @@ 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 &) + copy_triangulation (const dealii::Triangulation &old_tria) { - Assert (false, ExcNotImplemented()); + 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 (); } diff --git a/deal.II/source/lac/petsc_precondition.cc b/deal.II/source/lac/petsc_precondition.cc index 2df82399d2..370d5c99b5 100644 --- a/deal.II/source/lac/petsc_precondition.cc +++ b/deal.II/source/lac/petsc_precondition.cc @@ -483,6 +483,223 @@ 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 d01dae47f0..0cbb09d2bc 100644 --- a/deal.II/source/lac/petsc_solver.cc +++ b/deal.II/source/lac/petsc_solver.cc @@ -129,6 +129,138 @@ 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) { @@ -302,6 +434,9 @@ 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); } @@ -404,6 +539,9 @@ 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