>
>
> I've started working on new version MGTransferBlockMatrixFree and I now
> see the problem. MGTransferBlockMatrixFree is easy to rewrite using vector
> of dof handres etc, but PreconditionMG requires MGTransfer that have
> copy_to_mg and copy_from_mg that work on single dof handre. I am thinking
> of 3 way to resolve it:
> 1) Rewrite MGTransferBlockMatrixFree using vector of dof handler,
> implement PreconditionBlockMG and probably Multigrid will also need block
> version
> 2) Use standard MGTransfer for individual blocks and implement
> PreconditionBlockMG and BlockMultigrid to work with several MGTransfer s.
> 3) Since MGTransferBlockMatrixFree requires dof handler in build (const
> DoFHandler< dim, dim > &mg_dof) method, the pointers to dof_handlers may be
> stored and used in copy_from_mg and copy_to_mg.
> I think 3) is simplest but also risky solution. 2) seems quite rational to
> me.
>
>>
>> I've decided to implement approach 3). You can find it in attachment.
Second file contains test on block Laplace problem. Main difference between
this file and one of deal.II test are this lines:
std::vector< const MGConstrainedDoFs * > mg_constrained_dofs_vector(nb,
&mg_constrained_dofs );
std::vector< const DoFHandler<dim> * > mg_dofs_vector(nb, &dof );
MWichrowski::MGTransferBlockMatrixFree<dim,dim,number>
mg_transfer(mg_constrained_dofs_vector);
mg_transfer.build(mg_dofs_vector);
> Best,
> Michał
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/*
* MGTransferBlockMatrixFree.h
*
* Created on: Sep 26, 2017
* Author: mwichro
*/
#ifndef MGTRANSFERBLOCKMATRIXFREE_H_
#define MGTRANSFERBLOCKMATRIXFREE_H_
#include <deal.II/base/config.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
using namespace dealii;
namespace MWichrowski {
template <int dim, int spacedim, typename Number>
class MGTransferBlockMatrixFree : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
{
public:
/**
* Constructor without constraint matrices. Use this constructor only with
* discontinuous finite elements or with no local refinement.
*/
MGTransferBlockMatrixFree (const unsigned int &n_b) //= default;
:
matrix_free_transfer(n_b),
n_blocks(n_b)
{}
/**
* Constructor with constraints. Equivalent to the default constructor
* followed by initialize_constraints().
*/
MGTransferBlockMatrixFree (const std::vector< const MGConstrainedDoFs * > &mg_constrained_dofs);
/**
* Destructor.
*/
virtual ~MGTransferBlockMatrixFree () = default;
/**
* Initialize the constraints to be used in build().
*/
void initialize_constraints (const std::vector< const MGConstrainedDoFs* > &mg_constrained_dofs);
/**
* Reset the object to the state it had right after the default constructor.
*/
void clear ();
/**
* Actually build the information for the prolongation for each level.
*/
void build (const std::vector<const DoFHandler<dim,spacedim> * > &mg_dofs);
/**
* Prolongate a vector from level <tt>to_level-1</tt> to level
* <tt>to_level</tt> using the embedding matrices of the underlying finite
* element. The previous content of <tt>dst</tt> is overwritten.
*
* @param to_level The index of the level to prolongate to, which is the
* level of @p dst.
*
* @param src is a vector with as many elements as there are degrees of
* freedom on the coarser level involved.
*
* @param dst has as many elements as there are degrees of freedom on the
* finer level.
*/
virtual void prolongate (const unsigned int to_level,
LinearAlgebra::distributed::BlockVector<Number> &dst,
const LinearAlgebra::distributed::BlockVector<Number> &src) const;
/**
* Restrict a vector from level <tt>from_level</tt> to level
* <tt>from_level-1</tt> using the transpose operation of the prolongate()
* method. If the region covered by cells on level <tt>from_level</tt> is
* smaller than that of level <tt>from_level-1</tt> (local refinement), then
* some degrees of freedom in <tt>dst</tt> are active and will not be
* altered. For the other degrees of freedom, the result of the restriction
* is added.
*
* @param from_level The index of the level to restrict from, which is the
* level of @p src.
*
* @param src is a vector with as many elements as there are degrees of
* freedom on the finer level involved.
*
* @param dst has as many elements as there are degrees of freedom on the
* coarser level.
*/
virtual void restrict_and_add (const unsigned int from_level,
LinearAlgebra::distributed::BlockVector<Number> &dst,
const LinearAlgebra::distributed::BlockVector<Number> &src) const;
/**
* Transfer from a block-vector on the global grid to block-vectors defined
* on each of the levels separately for active degrees of freedom.
* In particular, for a globally refined mesh only the finest level in @p dst
* is filled as a plain copy of @p src. All the other level objects are left
* untouched.
*
* This function will initialize @p dst accordingly if needed as required by
* the Multigrid class.
*/
template <typename Number2>
void
copy_to_mg (const DoFHandler<dim,spacedim> &first_mg_dof,
MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
const LinearAlgebra::distributed::BlockVector<Number2> &src) const;
/**
* Transfer from multi-level block-vector to normal vector.
*/
template <typename Number2>
void
copy_from_mg (const DoFHandler<dim,spacedim> &first_mg_dof,
LinearAlgebra::distributed::BlockVector<Number2> &dst,
const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const;
/**
* Memory used by this object.
*/
std::size_t memory_consumption () const;
private:
/**
* A storage of non-block matrix-free version of transfer operation.
*/
std::vector<MGTransferMatrixFree<dim,Number> > matrix_free_transfer;
std::vector<const DoFHandler<dim,spacedim> * > dof_handlers;
unsigned int n_blocks;
bool has_dof_handlers;
};
template <int dim, int spacedim, typename Number>
MGTransferBlockMatrixFree<dim,spacedim,Number>::MGTransferBlockMatrixFree (const std::vector< const MGConstrainedDoFs * > &mg_constrained_dofs)
:
matrix_free_transfer(mg_constrained_dofs.size() ),
dof_handlers(mg_constrained_dofs.size()),
n_blocks(mg_constrained_dofs.size()),
has_dof_handlers(false)
{
initialize_constraints(mg_constrained_dofs);
}
template <int dim, int spacedim, typename Number>
void MGTransferBlockMatrixFree<dim,spacedim,Number>::initialize_constraints
(const std::vector< const MGConstrainedDoFs* > &mg_constrained_dofs)
{
AssertDimension(mg_constrained_dofs.size(),n_blocks );
for(unsigned int i=0; i<n_blocks; ++i)
matrix_free_transfer[i].initialize_constraints(*mg_constrained_dofs[i]);
}
template <int dim, int spacedim, typename Number>
void MGTransferBlockMatrixFree<dim,spacedim,Number>::clear ()
{
for(unsigned int i=0; i<n_blocks; ++i)
matrix_free_transfer[i].clear();
dof_handlers=std::vector<const DoFHandler<dim,spacedim> * >(n_blocks);
has_dof_handlers=false;
}
template <int dim, int spacedim, typename Number>
void MGTransferBlockMatrixFree<dim,spacedim,Number>::build
(const std::vector<const DoFHandler<dim,spacedim> * > &mg_dofs)
{
// matrix_free_transfer.build(mg_dof);
AssertDimension(mg_dofs.size(),n_blocks );
for(unsigned int i=0; i<n_blocks; ++i)
matrix_free_transfer[i].build(*mg_dofs[i]);
dof_handlers=mg_dofs;
has_dof_handlers=true;
}
template <int dim, int spacedim, typename Number>
void MGTransferBlockMatrixFree<dim,spacedim,Number>
::prolongate (const unsigned int to_level,
LinearAlgebra::distributed::BlockVector<Number> &dst,
const LinearAlgebra::distributed::BlockVector<Number> &src) const
{
AssertDimension(dst.n_blocks(), src.n_blocks());
AssertDimension(dst.n_blocks(), n_blocks);
for (unsigned int b = 0; b < n_blocks; ++b)
matrix_free_transfer[b].prolongate(to_level, dst.block(b), src.block(b));
}
template <int dim, int spacedim, typename Number>
void MGTransferBlockMatrixFree<dim,spacedim,Number>
::restrict_and_add (const unsigned int from_level,
LinearAlgebra::distributed::BlockVector<Number> &dst,
const LinearAlgebra::distributed::BlockVector<Number> &src) const
{
AssertDimension(dst.n_blocks(), src.n_blocks());
AssertDimension(dst.n_blocks(), n_blocks);
for (unsigned int b = 0; b < n_blocks; ++b)
matrix_free_transfer[b].restrict_and_add(from_level, dst.block(b), src.block(b));
}
template <int dim, int spacedim, typename Number>
std::size_t
MGTransferBlockMatrixFree<dim,spacedim,Number>::memory_consumption() const
{
std::size_t size=0;
for (unsigned int b = 0; b < n_blocks; ++b)
size+= matrix_free_transfer[b].memory_consumption();
return size;
}
template <int dim, int spacedim, typename Number>
template <typename Number2>
void
MGTransferBlockMatrixFree<dim,spacedim,Number>::
copy_to_mg (const DoFHandler<dim,spacedim> &first_mg_dof,
MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
const LinearAlgebra::distributed::BlockVector<Number2> &src) const
{
const unsigned int min_level = dst.min_level();
const unsigned int max_level = dst.max_level();
Assert(has_dof_handlers,ExcNotInitialized () );
Assert(first_mg_dof == *dof_handlers[0],
ExcMessage ("DoF handler does not match with DoF handler used in build"));
// this function is normally called within the Multigrid class with
// dst == defect level block vector. At first run this vector is not
// initialized. Do this below:
{
std::vector<const parallel::Triangulation<dim,spacedim> * > tria(n_blocks);
for (unsigned int b = 0; b < n_blocks; ++b)
tria[b]= (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
(& dof_handlers[b]->get_triangulation()));
for (unsigned int level = min_level; level <= max_level; ++level)
{
dst[level].reinit(n_blocks);
bool collect_size = false;
for (unsigned int b = 0; b < n_blocks; ++b)
{
LinearAlgebra::distributed::Vector<Number> &v = dst[level].block(b);
if (v.size() != dof_handlers[b]->locally_owned_mg_dofs(level).size() ||
v.local_size() != dof_handlers[b]->locally_owned_mg_dofs(level).n_elements())
{
v.reinit(dof_handlers[b]->locally_owned_mg_dofs(level), tria[b] != nullptr ?
tria[b]->get_communicator() : MPI_COMM_SELF);
collect_size = true;
}
else
v = 0.;
}
if (collect_size)
dst[level].collect_sizes ();
}
}
// FIXME: this a quite ugly as we need a temporary object:
MGLevelObject<LinearAlgebra::distributed::Vector<Number>> dst_non_block(min_level, max_level);
for (unsigned int b = 0; b < n_blocks; ++b)
{
for (unsigned int l = min_level; l <= max_level; ++l)
dst_non_block[l].reinit(dst[l].block(b));
matrix_free_transfer[b].copy_to_mg(*dof_handlers[b], dst_non_block, src.block(b));
for (unsigned int l = min_level; l <= max_level; ++l)
dst[l].block(b) = dst_non_block[l];
}
}
template <int dim, int spacedim, typename Number>
template <typename Number2>
void
MGTransferBlockMatrixFree<dim,spacedim,Number>::
copy_from_mg (const DoFHandler<dim,spacedim> &first_mg_dof,
LinearAlgebra::distributed::BlockVector<Number2> &dst,
const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src) const
{
Assert(has_dof_handlers,ExcNotInitialized () );
Assert(first_mg_dof == *dof_handlers[0],
ExcMessage ("DoF handler does not match with DoF handler used in build"));
const unsigned int min_level = src.min_level();
const unsigned int max_level = src.max_level();
for (unsigned int l = min_level; l <= max_level; ++l){
AssertDimension(src[l].n_blocks(), dst.n_blocks());
AssertDimension(src[l].n_blocks(), n_blocks);
}
// FIXME: this a quite ugly as we need a temporary object:
MGLevelObject<LinearAlgebra::distributed::Vector<Number>> src_non_block(min_level, max_level);
for (unsigned int b = 0; b < n_blocks; ++b)
{
for (unsigned int l = min_level; l <= max_level; ++l)
{
src_non_block[l].reinit(src[l].block(b));
src_non_block[l] = src[l].block(b);
}
matrix_free_transfer[b].copy_from_mg(*dof_handlers[b], dst.block(b), src_non_block);
}
}
} /* namespace MWichrowski */
#endif /* MGTRANSFERBLOCKMATRIXFREE_H_ */
// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 2017 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.
//
// ---------------------------------------------------------------------
// similar to parallel_multigrid_adaptive_06ref but using MGTransferBlockMatrixFree
// to solve block-diagonal matrix with Laplace operator on diagonals.
// As expected, when we use a block vector with a single block, we get
// the same results as the reference, non-block solution, i.e.
//
// DEAL:2d:cg::Starting value 21.93
// DEAL:2d:cg::Convergence step 7 value 1.961e-07
//
//
// What this test is really for is block operations. As expected we see
// exactly the same number of iterations and the ratio of \sqrt(2) in values, i.e.
//
// DEAL:2d:cg::Starting value 31.02
// DEAL:2d:cg::Convergence step 7 value 2.774e-07
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <iostream>
#include <complex>
#include <fstream>
#include"MGTransferBlockMatrixFree.h"
std::ofstream logfile("output");
using namespace dealii;
using namespace dealii::MatrixFreeOperators;
/**
* Block-Laplace operator with Block vector.
*/
template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1, typename BlockVectorType = LinearAlgebra::distributed::BlockVector<double> >
class BlockLaplace : public Subscriptor
{
public:
typedef typename BlockVectorType::value_type value_type;
typedef typename BlockVectorType::size_type size_type;
BlockLaplace ()
: Subscriptor()
{
}
void initialize (std::shared_ptr<const MatrixFree<dim,value_type>> data)
{
laplace.initialize(data);
}
void initialize (std::shared_ptr<const MatrixFree<dim,value_type>> data, const MGConstrainedDoFs &mg_constrained_dofs,
const unsigned int level)
{
laplace.initialize(data, mg_constrained_dofs, level);
}
void vmult_interface_down (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.vmult_interface_down(dst.block(b), src.block(b));
}
void vmult_interface_up (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.vmult_interface_up(dst.block(b), src.block(b));
}
void vmult (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.vmult(dst.block(b), src.block(b));
}
void Tvmult (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.Tvmult(dst.block(b), src.block(b));
}
void vmult_add (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.vmult_add(dst.block(b), src.block(b));
}
void Tvmult_add (BlockVectorType &dst, const BlockVectorType &src) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.Tvmult_add(dst.block(b), src.block(b));
}
void precondition_Jacobi (BlockVectorType &dst, const BlockVectorType &src, const value_type omega) const
{
for (unsigned int b = 0; b < src.n_blocks(); ++b)
laplace.precondition_Jacobi(dst.block(b), src.block(b), omega);
}
void compute_diagonal ()
{
laplace.compute_diagonal ();
}
virtual void clear()
{
laplace.clear();
}
private:
MatrixFreeOperators::LaplaceOperator<dim, fe_degree, n_q_points_1d, 1, typename BlockVectorType::BlockType> laplace;
};
template <typename MatrixType, typename Number>
class MGCoarseIterative : public MGCoarseGridBase<LinearAlgebra::distributed::BlockVector<Number> >
{
public:
MGCoarseIterative() {}
void initialize(const MatrixType &matrix)
{
coarse_matrix = &matrix;
}
virtual void operator() (const unsigned int level,
LinearAlgebra::distributed::BlockVector<Number> &dst,
const LinearAlgebra::distributed::BlockVector<Number> &src) const
{
ReductionControl solver_control (1e4, 1e-50, 1e-10);
SolverCG<LinearAlgebra::distributed::BlockVector<Number> > solver_coarse (solver_control);
solver_coarse.solve (*coarse_matrix, dst, src, PreconditionIdentity());
}
const MatrixType *coarse_matrix;
};
template <int dim, int fe_degree, int n_q_points_1d, typename number>
void do_test (const DoFHandler<dim> &dof, const unsigned int nb)
{
if (types_are_equal<number,float>::value == true)
{
deallog.push("float");
}
else
{
}
deallog << "Testing " << dof.get_fe().get_name();
deallog << std::endl;
deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs(dof, locally_relevant_dofs);
// Dirichlet BC
ZeroFunction<dim> zero_function;
typename FunctionMap<dim>::type dirichlet_boundary;
dirichlet_boundary[0] = &zero_function;
// fine-level constraints
ConstraintMatrix constraints;
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof, constraints);
VectorTools::interpolate_boundary_values(dof, dirichlet_boundary,
constraints);
constraints.close();
// level constraints:
MGConstrainedDoFs mg_constrained_dofs;
mg_constrained_dofs.initialize(dof, dirichlet_boundary);
MappingQ<dim> mapping(fe_degree+1);
BlockLaplace<dim,fe_degree,n_q_points_1d,LinearAlgebra::distributed::BlockVector<number>> fine_matrix;
std::shared_ptr<MatrixFree<dim,number> > fine_level_data(new MatrixFree<dim,number> ());
typename MatrixFree<dim,number>::AdditionalData fine_level_additional_data;
fine_level_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
fine_level_additional_data.tasks_block_size = 3;
fine_level_data->reinit (mapping, dof, constraints, QGauss<1>(n_q_points_1d),
fine_level_additional_data);
fine_matrix.initialize(fine_level_data);
fine_matrix.compute_diagonal();
LinearAlgebra::distributed::BlockVector<number> in(nb), sol(nb);
for (unsigned int b=0; b < nb; ++b)
{
fine_level_data->initialize_dof_vector(in.block(b));
fine_level_data->initialize_dof_vector(sol.block(b));
}
in.collect_sizes();
sol.collect_sizes();
// set constant rhs vector
{
// this is to make it consistent with parallel_multigrid_adaptive.cc
ConstraintMatrix hanging_node_constraints;
hanging_node_constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof, hanging_node_constraints);
hanging_node_constraints.close();
for (unsigned int i=0; i<in.block(0).local_size(); ++i)
if (!hanging_node_constraints.is_constrained(in.block(0).get_partitioner()->local_to_global(i)))
in.block(0).local_element(i) = 1.;
for (unsigned int b = 1; b < nb; ++b)
in.block(b) = in.block(0);
}
// set up multigrid in analogy to step-37
typedef BlockLaplace<dim,fe_degree,n_q_points_1d,LinearAlgebra::distributed::BlockVector<number>> LevelMatrixType;
MGLevelObject<LevelMatrixType> mg_matrices;
MGLevelObject<MatrixFree<dim,number> > mg_level_data;
mg_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
mg_level_data.resize(0, dof.get_triangulation().n_global_levels()-1);
for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
{
typename MatrixFree<dim,number>::AdditionalData mg_additional_data;
mg_additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
mg_additional_data.tasks_block_size = 3;
mg_additional_data.level_mg_handler = level;
ConstraintMatrix level_constraints;
IndexSet relevant_dofs;
DoFTools::extract_locally_relevant_level_dofs(dof, level,
relevant_dofs);
level_constraints.reinit(relevant_dofs);
level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level));
level_constraints.close();
mg_level_data[level].reinit (mapping, dof, level_constraints, QGauss<1>(n_q_points_1d),
mg_additional_data);
mg_matrices[level].initialize(std::make_shared<MatrixFree<dim,number> >(
mg_level_data[level]),
mg_constrained_dofs,
level);
mg_matrices[level].compute_diagonal();
}
MGLevelObject<MGInterfaceOperator<LevelMatrixType> > mg_interface_matrices;
mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1);
for (unsigned int level=0; level<dof.get_triangulation().n_global_levels(); ++level)
mg_interface_matrices[level].initialize(mg_matrices[level]);
// MGTransferBlockMatrixFree<dim,number> mg_transfer(mg_constrained_dofs);
std::vector< const MGConstrainedDoFs * > mg_constrained_dofs_vector(nb, &mg_constrained_dofs );
std::vector< const DoFHandler<dim> * > mg_dofs_vector(nb, &dof );
MWichrowski::MGTransferBlockMatrixFree<dim,dim,number> mg_transfer(mg_constrained_dofs_vector);
mg_transfer.build(mg_dofs_vector);
MGCoarseIterative<LevelMatrixType,number> mg_coarse;
mg_coarse.initialize(mg_matrices[0]);
typedef PreconditionJacobi<LevelMatrixType> SMOOTHER;
MGSmootherPrecondition<LevelMatrixType, SMOOTHER, LinearAlgebra::distributed::BlockVector<number> >
mg_smoother;
mg_smoother.initialize(mg_matrices, typename SMOOTHER::AdditionalData(0.8));
mg::Matrix<LinearAlgebra::distributed::BlockVector<number> >
mg_matrix(mg_matrices);
mg::Matrix<LinearAlgebra::distributed::BlockVector<number> >
mg_interface(mg_interface_matrices);
Multigrid<LinearAlgebra::distributed::BlockVector<number> > mg(
mg_matrix,
mg_coarse,
mg_transfer,
mg_smoother,
mg_smoother);
mg.set_edge_matrices(mg_interface, mg_interface);
PreconditionMG<dim, LinearAlgebra::distributed::BlockVector<number>,
MWichrowski::MGTransferBlockMatrixFree<dim,dim,number> >
preconditioner(dof, mg, mg_transfer);
{
// avoid output from inner (coarse-level) solver
deallog.depth_file(3);
ReductionControl control(30, 1e-20, 1e-7);
SolverCG<LinearAlgebra::distributed::BlockVector<number> > solver(control);
solver.solve(fine_matrix, sol, in, preconditioner);
}
if (types_are_equal<number,float>::value == true)
deallog.pop();
fine_matrix.clear();
for (unsigned int level = 0; level<dof.get_triangulation().n_global_levels(); ++level)
mg_matrices[level].clear();
}
template <int dim, int fe_degree>
void test (const unsigned int nbands = 1)
{
parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD,
Triangulation<dim>::limit_level_difference_at_vertices,
parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
GridGenerator::hyper_cube (tria);
tria.refine_global(6-dim);
const unsigned int n_runs = fe_degree == 1 ? 6-dim : 5-dim;
for (unsigned int i=0; i<n_runs; ++i)
{
for (typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
cell != tria.end(); ++cell)
if (cell->is_locally_owned() &&
(cell->center().norm() < 0.5 && (cell->level() < 5 ||
cell->center().norm() > 0.45)
||
(dim == 2 && cell->center().norm() > 1.2)))
cell->set_refine_flag();
tria.execute_coarsening_and_refinement();
FE_Q<dim> fe (fe_degree);
DoFHandler<dim> dof (tria);
dof.distribute_dofs(fe);
dof.distribute_mg_dofs(fe);
do_test<dim, fe_degree, fe_degree+1, double> (dof, nbands);
}
}
int main (int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
{
deallog.attach(logfile);
deallog << std::setprecision (4);
}
{
deallog.push("2d");
test<2,1>();
test<2,3>();
deallog.pop();
deallog.push("3d");
test<3,1>();
test<3,2>();
deallog.pop();
}
// 2 blocks
{
deallog.push("2d");
test<2,1>(2);
test<2,3>(2);
deallog.pop();
deallog.push("3d");
test<3,1>(2);
test<3,2>(2);
deallog.pop();
}
}
