From: bangerth Date: Thu, 21 Aug 2008 15:29:50 +0000 (+0000) Subject: Move tool functions all in the same place. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=34990e9eba19679d273daacbfd5f1f46333c3662;p=dealii-svn.git Move tool functions all in the same place. git-svn-id: https://svn.dealii.org/trunk@16626 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 9582f59b82..12c759d12f 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -437,6 +437,249 @@ namespace EquationData } +namespace LinearSolvers +{ + + + + // @sect3{Linear solvers and preconditioners} + + // This section introduces some + // objects that are used for the + // solution of the linear equations of + // Stokes system that we need to + // solve in each time step. The basic + // structure is still the same as + // in step-20, where Schur complement + // based preconditioners and solvers + // have been introduced, with the + // actual interface taken from step-22. + + // @sect4{The InverseMatrix class template} + + // This class is an interface to + // calculate the action of an + // "inverted" matrix on a vector + // (using the vmult + // operation) + // in the same way as the corresponding + // function in step-22: when the + // product of an object of this class + // is requested, we solve a linear + // equation system with that matrix + // using the CG method, accelerated + // by a preconditioner of (templated) class + // Preconditioner. + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); + + void vmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + const Preconditioner &preconditioner; + }; + + + template + InverseMatrix::InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (preconditioner) + {} + + + + template + void InverseMatrix::vmult (Vector &dst, + const Vector &src) const + { + SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); + SolverCG<> cg (solver_control); + + dst = 0; + + try + { + cg.solve (*matrix, dst, src, preconditioner); + } + catch (std::exception &e) + { + Assert (false, ExcMessage(e.what())); + } + } + + // @sect4{Schur complement preconditioner} + + // This is the implementation + // of the Schur complement + // preconditioner as described + // in the section on improved + // solvers in step-22. + // + // The basic + // concept of the preconditioner is + // different to the solution + // strategy used in step-20 and + // step-22. There, the Schur + // complement was used for a + // two-stage solution of the linear + // system. Recall that the process + // in the Schur complement solver is + // a Gaussian elimination of + // a 2x2 block matrix, where each + // block is solved iteratively. + // Here, the idea is to let + // an iterative solver act on the + // whole system, and to use + // a Schur complement for + // preconditioning. As usual when + // dealing with preconditioners, we + // don't intend to exacly set up a + // Schur complement, but rather use + // a good approximation to the + // Schur complement for the purpose of + // preconditioning. + // + // So the question is how we can + // obtain a good preconditioner. + // Let's have a look at the + // preconditioner matrix P + // acting on the block system, built + // as + // @f{eqnarray*} + // P^{-1} + // = + // \left(\begin{array}{cc} + // A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1} + // \end{array}\right) + // @f} + // using the Schur complement + // $S = B A^{-1} B^T$. If we apply + // this matrix in the solution of + // a linear system, convergence of + // an iterative Krylov-based solver + // will be governed by the matrix + // @f{eqnarray*} + // P^{-1}\left(\begin{array}{cc} + // A & B^T \\ B & 0 + // \end{array}\right) + // = + // \left(\begin{array}{cc} + // I & A^{-1} B^T \\ 0 & 0 + // \end{array}\right), + // @f} + // which turns out to be very simple. + // A GMRES solver based on exact + // matrices would converge in two + // iterations, since there are + // only two distinct eigenvalues. + // Such a preconditioner for the + // blocked Stokes system has been + // proposed by Silvester and Wathen, + // Fast iterative solution of + // stabilised Stokes systems part II. + // Using general block preconditioners. + // (SIAM J. Numer. Anal., 31 (1994), + // pp. 1352-1367). + // + // The deal.II users who have already + // gone through the step-20 and step-22 + // tutorials can certainly imagine + // how we're going to implement this. + // We replace the inverse matrices + // in $P^{-1}$ using the InverseMatrix + // class, and the inverse Schur + // complement will be approximated + // by the pressure mass matrix $M_p$. + // Having this in mind, we define a + // preconditioner class with a + // vmult functionality, + // which is all we need for the + // interaction with the usual solver + // functions further below in the + // program code. + // + // First the declarations. These + // are similar to the definition of + // the Schur complement in step-20, + // with the difference that we need + // some more preconditioners in + // the constructor. + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner (const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner); + + void vmult (BlockVector &dst, + const BlockVector &src) const; + + private: + const SmartPointer > stokes_matrix; + const SmartPointer, + PreconditionerMp > > m_inverse; +const PreconditionerA &a_preconditioner; + +mutable Vector tmp; + +}; + + + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner(const BlockSparseMatrix &S, + const InverseMatrix,PreconditionerMp> &Mpinv, + const PreconditionerA &Apreconditioner) + : + stokes_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (S.block(1,1).m()) + {} + + + // This is the vmult + // function. We implement + // the action of $P^{-1}$ as described + // above in three successive steps. + // The first step multiplies + // the velocity vector by a + // preconditioner of the matrix A. + // The resuling velocity vector + // is then multiplied by $B$ and + // subtracted from the pressure. + // This second step only acts on + // the pressure vector and is + // accomplished by the command + // SparseMatrix::residual. Next, + // we change the sign in the + // temporary pressure vector and + // finally multiply by the pressure + // mass matrix to get the final + // pressure vector. + template + void BlockSchurPreconditioner::vmult ( + BlockVector &dst, + const BlockVector &src) const + { + a_preconditioner.vmult (dst.block(0), src.block(0)); + stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); + tmp *= -1; + m_inverse->vmult (dst.block(1), tmp); + } +} + + // @sect3{The BoussinesqFlowProblem class template} @@ -473,6 +716,20 @@ class BoussinesqFlowProblem void output_results () const; void refine_mesh (const unsigned int max_grid_level); + static void compute_viscosity(const std::vector &old_temperature, + const std::vector &old_old_temperature, + const std::vector > &old_temperature_grads, + const std::vector > &old_old_temperature_grads, + const std::vector > &old_temperature_hessians, + const std::vector > &old_old_temperature_hessians, + const std::vector > &present_stokes_values, + const std::vector &gamma_values, + const double global_u_infty, + const double global_T_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step); + Triangulation triangulation; const unsigned int stokes_degree; @@ -509,8 +766,8 @@ class BoussinesqFlowProblem double old_time_step; unsigned int timestep_number; - boost::shared_ptr Amg_preconditioner; - boost::shared_ptr > Mp_preconditioner; + boost::shared_ptr Amg_preconditioner; + boost::shared_ptr > Mp_preconditioner; bool rebuild_stokes_matrix; bool rebuild_temperature_matrices; @@ -518,247 +775,6 @@ class BoussinesqFlowProblem }; - - - - // @sect3{Linear solvers and preconditioners} - - // This section introduces some - // objects that are used for the - // solution of the linear equations of - // Stokes system that we need to - // solve in each time step. The basic - // structure is still the same as - // in step-20, where Schur complement - // based preconditioners and solvers - // have been introduced, with the - // actual interface taken from step-22. - - // @sect4{The InverseMatrix class template} - - // This class is an interface to - // calculate the action of an - // "inverted" matrix on a vector - // (using the vmult - // operation) - // in the same way as the corresponding - // function in step-22: when the - // product of an object of this class - // is requested, we solve a linear - // equation system with that matrix - // using the CG method, accelerated - // by a preconditioner of (templated) class - // Preconditioner. -template -class InverseMatrix : public Subscriptor -{ - public: - InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer matrix; - const Preconditioner &preconditioner; -}; - - -template -InverseMatrix::InverseMatrix (const Matrix &m, - const Preconditioner &preconditioner) - : - matrix (&m), - preconditioner (preconditioner) -{} - - - -template -void InverseMatrix::vmult (Vector &dst, - const Vector &src) const -{ - SolverControl solver_control (src.size(), 1e-6*src.l2_norm()); - SolverCG<> cg (solver_control); - - dst = 0; - - try - { - cg.solve (*matrix, dst, src, preconditioner); - } - catch (std::exception &e) - { - Assert (false, ExcMessage(e.what())); - } -} - - // @sect4{Schur complement preconditioner} - - // This is the implementation - // of the Schur complement - // preconditioner as described - // in the section on improved - // solvers in step-22. - // - // The basic - // concept of the preconditioner is - // different to the solution - // strategy used in step-20 and - // step-22. There, the Schur - // complement was used for a - // two-stage solution of the linear - // system. Recall that the process - // in the Schur complement solver is - // a Gaussian elimination of - // a 2x2 block matrix, where each - // block is solved iteratively. - // Here, the idea is to let - // an iterative solver act on the - // whole system, and to use - // a Schur complement for - // preconditioning. As usual when - // dealing with preconditioners, we - // don't intend to exacly set up a - // Schur complement, but rather use - // a good approximation to the - // Schur complement for the purpose of - // preconditioning. - // - // So the question is how we can - // obtain a good preconditioner. - // Let's have a look at the - // preconditioner matrix P - // acting on the block system, built - // as - // @f{eqnarray*} - // P^{-1} - // = - // \left(\begin{array}{cc} - // A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1} - // \end{array}\right) - // @f} - // using the Schur complement - // $S = B A^{-1} B^T$. If we apply - // this matrix in the solution of - // a linear system, convergence of - // an iterative Krylov-based solver - // will be governed by the matrix - // @f{eqnarray*} - // P^{-1}\left(\begin{array}{cc} - // A & B^T \\ B & 0 - // \end{array}\right) - // = - // \left(\begin{array}{cc} - // I & A^{-1} B^T \\ 0 & 0 - // \end{array}\right), - // @f} - // which turns out to be very simple. - // A GMRES solver based on exact - // matrices would converge in two - // iterations, since there are - // only two distinct eigenvalues. - // Such a preconditioner for the - // blocked Stokes system has been - // proposed by Silvester and Wathen, - // Fast iterative solution of - // stabilised Stokes systems part II. - // Using general block preconditioners. - // (SIAM J. Numer. Anal., 31 (1994), - // pp. 1352-1367). - // - // The deal.II users who have already - // gone through the step-20 and step-22 - // tutorials can certainly imagine - // how we're going to implement this. - // We replace the inverse matrices - // in $P^{-1}$ using the InverseMatrix - // class, and the inverse Schur - // complement will be approximated - // by the pressure mass matrix $M_p$. - // Having this in mind, we define a - // preconditioner class with a - // vmult functionality, - // which is all we need for the - // interaction with the usual solver - // functions further below in the - // program code. - // - // First the declarations. These - // are similar to the definition of - // the Schur complement in step-20, - // with the difference that we need - // some more preconditioners in - // the constructor. -template -class BlockSchurPreconditioner : public Subscriptor -{ - public: - BlockSchurPreconditioner (const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner); - - void vmult (BlockVector &dst, - const BlockVector &src) const; - - private: - const SmartPointer > stokes_matrix; - const SmartPointer, - PreconditionerMp > > m_inverse; - const PreconditionerA &a_preconditioner; - - mutable Vector tmp; - -}; - -template -BlockSchurPreconditioner::BlockSchurPreconditioner( - const BlockSparseMatrix &S, - const InverseMatrix,PreconditionerMp> &Mpinv, - const PreconditionerA &Apreconditioner - ) - : - stokes_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (S.block(1,1).m()) -{ -} - - - // This is the vmult - // function. We implement - // the action of $P^{-1}$ as described - // above in three successive steps. - // The first step multiplies - // the velocity vector by a - // preconditioner of the matrix A. - // The resuling velocity vector - // is then multiplied by $B$ and - // subtracted from the pressure. - // This second step only acts on - // the pressure vector and is - // accomplished by the command - // SparseMatrix::residual. Next, - // we change the sign in the - // temporary pressure vector and - // finally multiply by the pressure - // mass matrix to get the final - // pressure vector. -template -void BlockSchurPreconditioner::vmult ( - BlockVector &dst, - const BlockVector &src) const -{ - a_preconditioner.vmult (dst.block(0), src.block(0)); - stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); - tmp *= -1; - m_inverse->vmult (dst.block(1), tmp); -} - - - // @sect3{BoussinesqFlowProblem class implementation} // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem} @@ -799,6 +815,151 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) + // @sect4{BoussinesqFlowProblem::get_maximal_velocity} +template +double BoussinesqFlowProblem::get_maximal_velocity () const +{ + const QGauss quadrature_formula(stokes_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (stokes_fe, quadrature_formula, update_values); + std::vector > stokes_values(n_q_points, + Vector(dim+1)); + double max_velocity = 0; + + typename DoFHandler::active_cell_iterator + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (stokes_solution, stokes_values); + + for (unsigned int q=0; q velocity; + for (unsigned int i=0; i +std::pair +BoussinesqFlowProblem::get_extrapolated_temperature_range () const +{ + QGauss quadrature_formula(temperature_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (temperature_fe, quadrature_formula, + update_values); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + + double min_temperature = (1. + time_step/old_time_step) * + old_temperature_solution.linfty_norm() + + + time_step/old_time_step * + old_old_temperature_solution.linfty_norm(), + max_temperature = -min_temperature; + + typename DoFHandler::active_cell_iterator + cell = temperature_dof_handler.begin_active(), + endc = temperature_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_temperature_solution, old_temperature_values); + fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); + + for (unsigned int q=0; q +double +BoussinesqFlowProblem:: +compute_viscosity(const std::vector &old_temperature, + const std::vector &old_old_temperature, + const std::vector > &old_temperature_grads, + const std::vector > &old_old_temperature_grads, + const std::vector > &old_temperature_hessians, + const std::vector > &old_old_temperature_hessians, + const std::vector > &present_stokes_values, + const std::vector &gamma_values, + const double global_u_infty, + const double global_T_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step) +{ + const double beta = 0.03; + const double alpha = 1; + + if (global_u_infty == 0) + return 5e-3 * cell_diameter; + + const unsigned int n_q_points = old_temperature.size(); + + // Stage 1: calculate residual + double max_residual = 0; + double max_velocity = 0; + + for (unsigned int q=0; q < n_q_points; ++q) + { + Tensor<1,dim> u; + for (unsigned int d=0; d::build_stokes_preconditioner () // copied over to Trilinos. we need to // keep the (1,1) block, though - Mp_preconditioner - = boost::shared_ptr > - (new SparseILU); + Mp_preconditioner = boost::shared_ptr > + (new SparseILU); Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1), SparseILU::AdditionalData()); @@ -1592,68 +1752,6 @@ void BoussinesqFlowProblem::assemble_stokes_system () -template -double compute_viscosity( - const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector > &old_temperature_hessians, - const std::vector > &old_old_temperature_hessians, - const std::vector > &present_stokes_values, - const std::vector &gamma_values, - const double global_u_infty, - const double global_T_variation, - const double global_Omega_diameter, - const double cell_diameter, - const double old_time_step -) -{ - const double beta = 0.03; - const double alpha = 1; - - if (global_u_infty == 0) - return 5e-3 * cell_diameter; - - const unsigned int n_q_points = old_temperature.size(); - - // Stage 1: calculate residual - double max_residual = 0; - double max_velocity = 0; - - for (unsigned int q=0; q < n_q_points; ++q) - { - Tensor<1,dim> u; - for (unsigned int d=0; d::assemble_temperature_system () temperature_fe_values.reinit (cell); stokes_fe_values.reinit (stokes_cell); - temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values); - temperature_fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); + temperature_fe_values.get_function_values (old_temperature_solution, + old_temperature_values); + temperature_fe_values.get_function_values (old_old_temperature_solution, + old_old_temperature_values); - temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads); - temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads); + temperature_fe_values.get_function_gradients (old_temperature_solution, + old_temperature_grads); + temperature_fe_values.get_function_gradients (old_old_temperature_solution, + old_old_temperature_grads); - temperature_fe_values.get_function_hessians (old_temperature_solution, old_temperature_hessians); - temperature_fe_values.get_function_hessians (old_old_temperature_solution, old_old_temperature_hessians); + temperature_fe_values.get_function_hessians (old_temperature_solution, + old_temperature_hessians); + temperature_fe_values.get_function_hessians (old_old_temperature_solution, + old_old_temperature_hessians); temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(), gamma_values); - stokes_fe_values.get_function_values (stokes_solution, present_stokes_values); + stokes_fe_values.get_function_values (stokes_solution, + present_stokes_values); - // build matrix contributions - - // define diffusion. take the - // maximum of what we really - // want and the minimal amount - // of diffusion (determined - // impirically) to keep the - // scheme stable const double nu = compute_viscosity (old_temperature_values, old_old_temperature_values, @@ -2193,87 +2290,6 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) - // @sect4{BoussinesqFlowProblem::get_maximal_velocity} -template -double BoussinesqFlowProblem::get_maximal_velocity () const -{ - const QGauss quadrature_formula(stokes_degree+2); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (stokes_fe, quadrature_formula, update_values); - std::vector > stokes_values(n_q_points, - Vector(dim+1)); - double max_velocity = 0; - - typename DoFHandler::active_cell_iterator - cell = stokes_dof_handler.begin_active(), - endc = stokes_dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - fe_values.get_function_values (stokes_solution, stokes_values); - - for (unsigned int q=0; q velocity; - for (unsigned int i=0; i -std::pair -BoussinesqFlowProblem::get_extrapolated_temperature_range () const -{ - QGauss quadrature_formula(temperature_degree+2); - const unsigned int n_q_points = quadrature_formula.size(); - - FEValues fe_values (temperature_fe, quadrature_formula, - update_values); - std::vector old_temperature_values(n_q_points); - std::vector old_old_temperature_values(n_q_points); - - double min_temperature = (1. + time_step/old_time_step) * - old_temperature_solution.linfty_norm() - + - time_step/old_time_step * - old_old_temperature_solution.linfty_norm(), - max_temperature = -min_temperature; - - typename DoFHandler::active_cell_iterator - cell = temperature_dof_handler.begin_active(), - endc = temperature_dof_handler.end(); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - fe_values.get_function_values (old_temperature_solution, old_temperature_values); - fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values); - - for (unsigned int q=0; q void BoussinesqFlowProblem::run ()