From: Guido Kanschat Date: Tue, 2 Aug 2016 12:55:01 +0000 (+0200) Subject: parallel RelaxationBlock X-Git-Tag: v8.5.0-rc1~792^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=264b4543468f577f0269e6a40c4dd9b23bb826fc;p=dealii.git parallel RelaxationBlock - template RelaxationBlock::step on vector type - instantiate for Trilinos - add test --- diff --git a/include/deal.II/lac/relaxation_block.h b/include/deal.II/lac/relaxation_block.h index ffeb83c3ec..ce5597255d 100644 --- a/include/deal.II/lac/relaxation_block.h +++ b/include/deal.II/lac/relaxation_block.h @@ -210,11 +210,11 @@ protected: * vectors). For the Jacobi step, the calling function must copy @p dst to * @p pref after this. */ - template + template void do_step ( - Vector &dst, - const Vector &prev, - const Vector &src, + VECTOR &dst, + const VECTOR &prev, + const VECTOR &src, const bool backward) const; /** * Pointer to the matrix. Make sure that the matrix exists as long as this @@ -300,14 +300,14 @@ public: /** * Perform one step of the Jacobi iteration. */ - template - void step (Vector &dst, const Vector &rhs) const; + template + void step (VECTOR &dst, const VECTOR &rhs) const; /** * Perform one step of the Jacobi iteration. */ - template - void Tstep (Vector &dst, const Vector &rhs) const; + template + void Tstep (VECTOR &dst, const VECTOR &rhs) const; /** * Return the memory allocated in this object. @@ -385,14 +385,14 @@ public: /** * Perform one step of the SOR iteration. */ - template - void step (Vector &dst, const Vector &rhs) const; + template + void step (VECTOR &dst, const VECTOR &rhs) const; /** * Perform one step of the transposed SOR iteration. */ - template - void Tstep (Vector &dst, const Vector &rhs) const; + template + void Tstep (VECTOR &dst, const VECTOR &rhs) const; }; @@ -459,16 +459,16 @@ public: using RelaxationBlock::inverse_svd; using PreconditionBlockBase::log_statistics; /** - * Perform one step of the SOR iteration. + * Perform one step of the SSOR iteration. */ - template - void step (Vector &dst, const Vector &rhs) const; + template + void step (VECTOR &dst, const VECTOR &rhs) const; /** - * Perform one step of the transposed SOR iteration. + * Perform one step of the transposed SSOR iteration. */ - template - void Tstep (Vector &dst, const Vector &rhs) const; + template + void Tstep (VECTOR &dst, const VECTOR &rhs) const; }; diff --git a/include/deal.II/lac/relaxation_block.templates.h b/include/deal.II/lac/relaxation_block.templates.h index 27c7254427..70a99df7a7 100644 --- a/include/deal.II/lac/relaxation_block.templates.h +++ b/include/deal.II/lac/relaxation_block.templates.h @@ -153,18 +153,18 @@ RelaxationBlock::invert_diagblocks () template -template +template inline void -RelaxationBlock::do_step (Vector &dst, - const Vector &prev, - const Vector &src, +RelaxationBlock::do_step (VECTOR &dst, + const VECTOR &prev, + const VECTOR &src, const bool backward) const { Assert (additional_data->invert_diagonal, ExcNotImplemented()); const MatrixType &M=*this->A; - Vector b_cell, x_cell; + Vector b_cell, x_cell; const bool permutation_empty = additional_data->order.size() == 0; const unsigned int n_permutations = (permutation_empty) @@ -213,19 +213,20 @@ RelaxationBlock::do_step (Vector &dst, dst(row->column()) += additional_data->relaxation * x_cell(row_cell); } } + dst.compress(dealii::VectorOperation::add); } //----------------------------------------------------------------------// template -template +template void RelaxationBlockJacobi::step -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { - GrowingVectorMemory > mem; - typename VectorMemory >::Pointer aux = mem; + GrowingVectorMemory mem; + typename VectorMemory::Pointer aux = mem; aux->reinit(dst, false); *aux = dst; this->do_step(dst, *aux, src, false); @@ -233,13 +234,13 @@ void RelaxationBlockJacobi::step template -template +template void RelaxationBlockJacobi::Tstep -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { - GrowingVectorMemory > mem; - typename VectorMemory >::Pointer aux = mem; + GrowingVectorMemory mem; + typename VectorMemory::Pointer aux = mem; aux->reinit(dst, false); *aux = dst; this->do_step(dst, *aux, src, true); @@ -249,20 +250,20 @@ void RelaxationBlockJacobi::Tstep //----------------------------------------------------------------------// template -template +template void RelaxationBlockSOR::step -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { this->do_step(dst, dst, src, false); } template -template +template void RelaxationBlockSOR::Tstep -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { this->do_step(dst, dst, src, true); } @@ -271,10 +272,10 @@ void RelaxationBlockSOR::Tstep //----------------------------------------------------------------------// template -template +template void RelaxationBlockSSOR::step -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { this->do_step(dst, dst, src, false); this->do_step(dst, dst, src, true); @@ -282,10 +283,10 @@ void RelaxationBlockSSOR::step template -template +template void RelaxationBlockSSOR::Tstep -(Vector &dst, - const Vector &src) const +(VECTOR &dst, + const VECTOR &src) const { this->do_step(dst, dst, src, true); this->do_step(dst, dst, src, false); diff --git a/source/lac/relaxation_block.cc b/source/lac/relaxation_block.cc index 3ab39515f8..f75a3684fd 100644 --- a/source/lac/relaxation_block.cc +++ b/source/lac/relaxation_block.cc @@ -15,6 +15,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN #include "relaxation_block.inst" diff --git a/source/lac/relaxation_block.inst.in b/source/lac/relaxation_block.inst.in index a013316e2c..47cf6b86fb 100644 --- a/source/lac/relaxation_block.inst.in +++ b/source/lac/relaxation_block.inst.in @@ -24,7 +24,7 @@ for (S1, S2 : REAL_SCALARS) } -for (S1, S2, S3 : REAL_SCALARS) +for (S1, S2, S3: REAL_SCALARS) { // ------------ RelaxationBlockJacobi ----------------- template @@ -51,3 +51,30 @@ for (S1, S2, S3 : REAL_SCALARS) (Vector &, const Vector &) const; } + +for (S2 : REAL_SCALARS) +{ +#ifdef DEAL_II_WITH_TRILINOS + template class RelaxationBlock; + template class RelaxationBlockJacobi; + template class RelaxationBlockSOR; + template class RelaxationBlockSSOR; + + template void RelaxationBlockJacobi::step + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; + template void RelaxationBlockJacobi::Tstep + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; + + template void RelaxationBlockSOR::step + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; + template void RelaxationBlockSOR::Tstep + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; + +// ------------ RelaxationBlockSSOR ----------------- + template void RelaxationBlockSSOR::step + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; + template void RelaxationBlockSSOR::Tstep + (TrilinosWrappers::MPI::Vector &, const TrilinosWrappers::MPI::Vector &) const; +#endif +} + diff --git a/tests/mpi/step-39-block.cc b/tests/mpi/step-39-block.cc new file mode 100644 index 0000000000..e24b17e9a3 --- /dev/null +++ b/tests/mpi/step-39-block.cc @@ -0,0 +1,746 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Author: Guido Kanschat, Texas A&M University, 2009 + */ + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace Step39 +{ + using namespace dealii; + + Functions::SlitSingularityFunction<2> exact_solution; + + + + + template + class MatrixIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void MatrixIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values()); + } + + + template + void MatrixIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::nitsche_matrix( + dinfo.matrix(0,false).matrix, info.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); + } + + template + void MatrixIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); + LocalIntegrators::Laplace::ip_matrix( + dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, + dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); + } + + template + class RHSIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void RHSIntegrator::cell(MeshWorker::DoFInfo &, typename MeshWorker::IntegrationInfo &) const + {} + + + template + void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + Vector &local_vector = dinfo.vector(0).block(0); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void RHSIntegrator::face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + typename MeshWorker::IntegrationInfo &, + typename MeshWorker::IntegrationInfo &) const + {} + + + template + class Estimator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void Estimator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + const std::vector > &DDuh = info.hessians[0][0]; + for (unsigned k=0; kdiameter() * trace(DDuh[k]); + dinfo.value(0) += t*t * fe.JxW(k); + } + dinfo.value(0) = std::sqrt(dinfo.value(0)); + } + + template + void Estimator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_values); + + const std::vector &uh = info.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void Estimator::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + const std::vector > &Duh1 = info1.gradients[0][0]; + const std::vector > &Duh2 = info2.gradients[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + const double h = dinfo1.face->measure(); + + for (unsigned k=0; kis_locally_owned()) + dinfo1.value(0) = 0.0; + if (!dinfo2.cell->is_locally_owned()) + dinfo2.value(0) = 0.0; + } + + + + template + class ErrorIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void ErrorIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + std::vector > exact_gradients(fe.n_quadrature_points); + std::vector exact_values(fe.n_quadrature_points); + + exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector > &Duh = info.gradients[0][0]; + const std::vector &uh = info.values[0][0]; + + for (unsigned k=0; k + void ErrorIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector exact_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector &uh = info.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); + + for (unsigned k=0; k + void ErrorIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + + const unsigned int deg = fe.get_fe().tensor_degree(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0; k + class InteriorPenaltyProblem + { + public: + typedef MeshWorker::IntegrationInfo CellInfo; + + InteriorPenaltyProblem(const FiniteElement &fe); + + void run(unsigned int n_steps); + + private: + void setup_system (); + void assemble_matrix (); + void assemble_mg_matrix (); + void assemble_right_hand_side (); + void error (); + double estimate (); + void solve (); + void output_results (const unsigned int cycle) const; + + parallel::distributed::Triangulation triangulation; + const MappingQGeneric mapping; + const FiniteElement &fe; + DoFHandler dof_handler; + + IndexSet locally_relevant_set; + + TrilinosWrappers::SparseMatrix matrix; + TrilinosWrappers::MPI::Vector solution; + TrilinosWrappers::MPI::Vector right_hand_side; + BlockVector estimates; + + MGLevelObject mg_matrix; + + MGLevelObject mg_matrix_dg_down; + MGLevelObject mg_matrix_dg_up; + }; + + + template + InteriorPenaltyProblem::InteriorPenaltyProblem(const FiniteElement &fe) + : + triangulation (MPI_COMM_WORLD,Triangulation:: + limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy), + mapping(1), + fe(fe), + dof_handler(triangulation), + estimates(1) + { + GridGenerator::hyper_L(triangulation, -1, 1); + } + + + template + void + InteriorPenaltyProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs (fe); + + DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_set); + solution.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + right_hand_side.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); + + DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); + matrix.reinit(dof_handler.locally_owned_dofs(), c_sparsity, MPI_COMM_WORLD, true); + + const unsigned int n_levels = triangulation.n_global_levels(); + mg_matrix.resize(0, n_levels-1); + mg_matrix.clear(); + mg_matrix_dg_up.resize(0, n_levels-1); + mg_matrix_dg_up.clear(); + mg_matrix_dg_down.resize(0, n_levels-1); + mg_matrix_dg_down.clear(); + + for (unsigned int level=mg_matrix.min_level(); + level<=mg_matrix.max_level(); ++level) + { + DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level); + mg_matrix[level].reinit(dof_handler.locally_owned_mg_dofs(level), + dof_handler.locally_owned_mg_dofs(level), + c_sparsity, + MPI_COMM_WORLD, true); + + if (level>0) + { + DynamicSparsityPattern ci_sparsity; + ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level); + + mg_matrix_dg_up[level].reinit(dof_handler.locally_owned_mg_dofs(level-1), + dof_handler.locally_owned_mg_dofs(level), + ci_sparsity, + MPI_COMM_WORLD, true); + mg_matrix_dg_down[level].reinit(dof_handler.locally_owned_mg_dofs(level-1), + dof_handler.locally_owned_mg_dofs(level), + ci_sparsity, + MPI_COMM_WORLD, true); + } + } + } + + + template + void + InteriorPenaltyProblem::assemble_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::MatrixSimple assembler; + assembler.initialize(matrix); + + FilteredIterator::active_cell_iterator> + begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); + FilteredIterator::active_cell_iterator> + end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + begin, end, + dof_info, info_box, + integrator, assembler); + + matrix.compress(VectorOperation::add); + } + + + template + void + InteriorPenaltyProblem::assemble_mg_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::MGMatrixSimple assembler; + assembler.initialize(mg_matrix); + assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + + FilteredIterator::level_cell_iterator> + begin(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.begin()); + FilteredIterator::level_cell_iterator> + end(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.end()); + + MatrixIntegrator integrator; + MeshWorker::integration_loop ( + begin, end, + dof_info, info_box, + integrator, assembler); + + for (unsigned int level=mg_matrix.min_level(); level <= mg_matrix.max_level(); ++level) + { + mg_matrix[level].compress(VectorOperation::add); + if (level > mg_matrix.min_level()) + { + mg_matrix_dg_up[level].compress(VectorOperation::add); + mg_matrix_dg_down[level].compress(VectorOperation::add); + } + } + } + + + template + void + InteriorPenaltyProblem::assemble_right_hand_side() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::ResidualSimple assembler; + AnyData data; + data.add(&right_hand_side, "RHS"); + assembler.initialize(data); + + FilteredIterator::active_cell_iterator> + begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); + FilteredIterator::active_cell_iterator> + end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); + + + RHSIntegrator integrator; + MeshWorker::integration_loop(begin, end, + dof_info, info_box, + integrator, assembler); + + right_hand_side.compress(VectorOperation::add); + right_hand_side *= -1.; + } + + + template + void + InteriorPenaltyProblem::solve() + { + SolverControl control(1000, 1.e-12); + SolverCG solver(control); + + MGTransferPrebuilt mg_transfer; + mg_transfer.build_matrices(dof_handler); + + SolverControl coarse_solver_control (1000, 1e-10, false, false); + SolverCG coarse_solver(coarse_solver_control); + PreconditionIdentity identity; + TrilinosWrappers::SparseMatrix &coarse_matrix = mg_matrix[0]; + MGCoarseGridLACIteration,TrilinosWrappers::MPI::Vector> + coarse_grid_solver(coarse_solver, coarse_matrix, identity); + + typedef RelaxationBlockJacobi + Smoother; + + MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_levels() - 1); + mg::SmootherRelaxation mg_smoother; + + for (unsigned int l = smoother_data.min_level() + 1; l <= smoother_data.max_level(); ++l) + { + DoFTools::make_cell_patches(smoother_data[l].block_list, this->dof_handler, l); + smoother_data[l].block_list.compress(); + smoother_data[l].relaxation = 0.7; + smoother_data[l].inversion = PreconditionBlockBase::svd; + } + + mg_smoother.initialize(mg_matrix, smoother_data); + mg_smoother.set_steps(1); + mg_smoother.set_variable(true); + + mg::Matrix mgmatrix(mg_matrix); + mg::Matrix mgdown(mg_matrix_dg_down); + mg::Matrix mgup(mg_matrix_dg_up); + + Multigrid mg(dof_handler, mgmatrix, + coarse_grid_solver, mg_transfer, + mg_smoother, mg_smoother); + mg.set_edge_flux_matrices(mgdown, mgup); + + PreconditionMG > + preconditioner(dof_handler, mg, mg_transfer); + solver.solve(matrix, solution, right_hand_side, preconditioner); + } + + + template + double + InteriorPenaltyProblem::estimate() + { + TrilinosWrappers::MPI::Vector ghost; + ghost.reinit(locally_relevant_set, MPI_COMM_WORLD); + ghost = solution; + + std::vector old_user_indices; + triangulation.save_user_indices(old_user_indices); + + estimates.block(0).reinit(triangulation.n_active_cells()); + unsigned int i=0; + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell,++i) + cell->set_user_index(i); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + AnyData solution_data; + solution_data.add(&ghost, "solution"); + + info_box.cell_selector.add("solution", false, false, true); + info_box.boundary_selector.add("solution", true, true, false); + info_box.face_selector.add("solution", true, true, false); + + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data, solution); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + AnyData out_data; + out_data.add*>(&estimates, "cells"); + assembler.initialize(out_data, false); + + Estimator integrator; + MeshWorker::LoopControl lctrl; + // assemble all faces adjacent to ghost cells to get the full + // information for all own cells without communication + lctrl.faces_to_ghost = MeshWorker::LoopControl::both; + + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler, lctrl); + + triangulation.load_user_indices(old_user_indices); + // estimates is a BlockVector (so serial) on each processor + // with one entry per active cell. Note that only the locally owned + // cells are !=0, so summing the contributions of l2_norm() over all + // processors is the right way to do this. + double local_norm = estimates.block(0).l2_norm(); + local_norm *= local_norm; + return std::sqrt(Utilities::MPI::sum(local_norm, MPI_COMM_WORLD)); + } + + template + void + InteriorPenaltyProblem::run(unsigned int n_steps) + { + deallog << "Element: " << fe.get_name() << std::endl; + for (unsigned int s=0; s fe1(2); + InteriorPenaltyProblem<2> test1(fe1); + test1.run(6); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/mpi/step-39-block.with_trilinos=true.mpirun=3.output b/tests/mpi/step-39-block.with_trilinos=true.mpirun=3.output new file mode 100644 index 0000000000..0519ecba6e --- /dev/null +++ b/tests/mpi/step-39-block.with_trilinos=true.mpirun=3.output @@ -0,0 +1 @@ + \ No newline at end of file