From: kronbichler Date: Thu, 14 Aug 2008 20:27:39 +0000 (+0000) Subject: Changed interface to Trilinos to make it generic. Now the initialize call provides... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3d251afefbe16beffb7f1cfc1cd8fcd7ffb79327;p=dealii-svn.git Changed interface to Trilinos to make it generic. Now the initialize call provides some options to be able to adapt the preconditioner to the problem. git-svn-id: https://svn.dealii.org/trunk@16555 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 9e7a0e84ee..d726a58865 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -138,8 +138,11 @@ class TrilinosAmgPreconditioner : public Subscriptor TrilinosAmgPreconditioner (); void initialize (const SparseMatrix &preconditioner_matrix, - const DoFHandler &dof_handler, - const std::vector &precondition_components); + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool higher_order_elements, + const bool elliptic, + const bool output_details); void vmult (Vector &dst, const Vector &src) const; @@ -160,9 +163,12 @@ TrilinosAmgPreconditioner::TrilinosAmgPreconditioner () template void TrilinosAmgPreconditioner::initialize ( - const SparseMatrix &preconditioner_matrix, - const DoFHandler &dof_handler, - const std::vector &precondition_components + const SparseMatrix &preconditioner_matrix, + const std::vector &null_space, + const unsigned int null_space_dimension, + const bool elliptic, + const bool higher_order_elements, + const bool output_details ) { unsigned int n_u = preconditioner_matrix.m(); @@ -229,106 +235,52 @@ void TrilinosAmgPreconditioner::initialize ( Matrix->FillComplete(); } - // And now build the AMG - // preconditioner. - const bool output_amg_info = false; + // Build the AMG preconditioner. Teuchos::ParameterList MLList; - // set default values for - // smoothed aggregation in MLList, - // the standard choice for elliptic - // (Laplace-type) problems - ML_Epetra::SetDefaults("SA",MLList); + // The implementation is able + // to distinguish between + // matrices from elliptic problems + // and convection dominated + // problems. We use the standard + // options for elliptic problems, + // except that we use a + // Chebyshev smoother instead + // of a symmetric Gauss-Seidel + // smoother. For most elliptic + // problems, Chebyshev is better + // than Gauss-Seidel (SSOR). + if (elliptic) + { + ML_Epetra::SetDefaults("SA",MLList); + MLList.set("smoother: type", "Chebyshev"); + MLList.set("smoother: sweeps", 4); + } + else + ML_Epetra::SetDefaults("NSSA",MLList); - if (output_amg_info) + if (output_details) MLList.set("ML output", 10); else MLList.set("ML output", 0); - // modify some AMG parameters from default to - // get better performance on FE induced Laplace - // type matrices - MLList.set("max levels",10); - MLList.set("increasing or decreasing", "increasing"); - MLList.set("aggregation: type", "MIS"); - MLList.set("smoother: type", "Chebyshev"); - MLList.set("smoother: sweeps", 4); - MLList.set("smoother: damping factor", 4./3.); - MLList.set("smoother: pre or post", "both"); - MLList.set("coarse: type","Amesos-KLU"); + if (higher_order_elements) + MLList.set("aggregation: type", "MIS"); - unsigned int n_total_components = precondition_components.size(); - Assert (n_total_components == dof_handler.get_fe().n_components(), - ExcDimensionMismatch (n_total_components, - dof_handler.get_fe().n_components() ) ); + Assert (n_u * null_space_dimension == null_space.size(), + ExcDimensionMismatch(n_u * null_space_dimension, + null_space.size())); - unsigned int n_precondition_components = 0; - for (unsigned int counter = 0; counter < n_total_components; ++counter) - n_precondition_components += precondition_components[counter]; - - // Build null space, i.e. build - // n_preconditioner_components vectors - // of ones in each velocity component. - std::vector null_vectors (n_precondition_components * n_u, 0.); - - if (n_precondition_components > 1) + if (null_space_dimension > 1) { - // TODO: is this n_ud really necessary? - // Or does it just bring new difficulties - // and confusion into play and limits - // the usability? - unsigned int n_ud = n_u/n_precondition_components; - Assert (n_ud * n_precondition_components == n_u, - ExcMessage("Cannot detect single problem components! " - "No. dofs not the same in all components!")); - - std::vector precondition_dof_list (dof_handler.n_dofs(), false); - std::vector precondition_mask (n_total_components, false); - unsigned int d = 0; - - for (unsigned int component=0; component < n_total_components; - component++) - { - if (precondition_components[component]) - { - precondition_mask[component] = true; - DoFTools::extract_dofs (dof_handler, precondition_mask, - precondition_dof_list); - precondition_mask[component] = false; - - // TODO: The current implementation - // assumes that we are working on - // the first components of a system when - // writing into the null vector. - // Change this to the general case, - // probably use something similar as - // for block vectors. - unsigned int counter = 0; - for (unsigned int i=0; iPrintUnused(0); } @@ -1706,13 +1658,44 @@ void BoussinesqFlowProblem::assemble_system () boost::shared_ptr > (new TrilinosAmgPreconditioner()); - std::vector prec_comp(dof_handler.get_fe().n_components(), false); - for (unsigned int component=0; component < dim; ++component) - prec_comp[component] = true; + const unsigned int n_u = preconditioner_matrix.block(0,0).m(); + std::vector null_space (dim * n_u, 0.); + + std::vector precondition_dof_list (dof_handler.n_dofs(), false); + std::vector precondition_mask (dim + 2, false); + + for (unsigned int component=0; component < dim; component++) + { + precondition_mask[component] = true; + DoFTools::extract_dofs (dof_handler, precondition_mask, + precondition_dof_list); + precondition_mask[component] = false; + + // TODO: The current implementation + // assumes that we are working on + // the first components of a system when + // writing into the null vector. + // Change this to the general case, + // probably use something similar as + // for block vectors. + unsigned int counter = 0; + for (unsigned int i=0; iinitialize(preconditioner_matrix.block(0,0), - dof_handler, - prec_comp); + null_space, dim, true, true, false); // TODO: we could throw away the (0,0) // block here since things have been