It looks like it was hidden bug inside my code. I've wrote the coarse grid
direct solver for matrix-free operators. It ignores sparsity pattern and
assembles the matrix by doing N vmults. The interface may be extended for
optimalization. I think it might be an useful tool. Usage, for step-37:
constructor needs vector partitioning:
CoarseDirectSolver<double> coarse_solver(owned_dofs_u,relevant_dofs_u );
This
coarse_solver.initialize(mg_matrices[0]/*, dof_handler*/);
or
coarse_solver.initialize(mg_matrices[0], dof_handler);
assembles the marix. In second case the matrix is initilized with sparsity
pattern, it could be optimized to reduce the number of vmults.
Michał
W dniu wtorek, 10 marca 2020 18:28:15 UTC+1 użytkownik Michał Wichrowski
napisał:
>
>
>
> W dniu wtorek, 10 marca 2020 17:33:43 UTC+1 użytkownik Wolfgang Bangerth
> napisał:
>>
>> On 3/9/20 8:26 AM, Michał Wichrowski wrote:
>> > I've got matrix-free multigrid solver for Stokes problem. The main
>> > bottleneck is solution of coarse problem, so I tried to assemble the
>> > regular sparse matrix and use direct solver. Since the coarse problem
>> is
>> > (relatively) small, I used vmults by unit vector to obtain columns of
>> > the matrix.
>>
>> This is unrelated to your original question, but worth mentioning
>> anyway: Conceptually, you are multiplying a unit vector by a sparse
>> matrix, and what you should get is a vector with nonzero entries in only
>> those rows for which the matrix has a nonzero entry in the column
>> corresponding to the unit vector.
>>
>> This is inefficient, because you will have to do N such matrix-vector
>> products. But if you know which entries of the matrix are non-zero
>> (which you do, because you know which DoFs couple with each other), then
>> you can multiply with a sum of unit vectors, get an output vector, and
>> you will know which entries of the output vector correspond to the
>> product with each of these unit vectors. You should be able to get the
>> number of matrix-vector products down from N to N/100 or N/1000 for
>> sufficiently large problems -- i.e., *much much faster*.
>>
>> (In fact, I think that we should at one point add a function that takes
>> an operator and that computes a matrix representation. If we would also
>> give that function a SparsityPattern, it could do the algorithm outlined
>> above automatically. As always, any help would be welcome!)
>>
>
> Yes, I thought about that, but I my case is not worth an effort. I'm doing
> it to assemble coarse matrix, so the problem is relatively small (N_dof
> around 500). The assemble time is unnoticeable comparing to solution time,
> even in case of small problems. Moreover, I need to apply inverse of
> some-kind Schur complement matrix (on coarse level only), and for that it
> would be extremely hard to guess sparsity pattern (it is still sparse).
> Fortunately, this matrix is even smaller.
>
> Your strategy (if implemented correctly) would require constant number of
> matrix-vector multiplications so the complexity will be linear. The number
> of vmults will depend on the mesh connectivity and type of FE. So the
> overall complexity will be linear with the problem size (instead (O(N^2),
> as I understand from Your email). But indeed, it will be much faster.
> Best,
> Michał
>
>>
>> Best
>> W.
>>
>> --
>> ------------------------------------------------------------------------
>> Wolfgang Bangerth email: [email protected]
>> www: http://www.math.colostate.edu/~bangerth/
>>
>
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/f8204afb-1750-4459-a7d1-67c33e4a256e%40googlegroups.com.
/*
* assemble_coarse_matrix.h
*
* Created on: Mar 10, 2020
* Author: mwichro
*/
#ifndef ASSEMBLE_COARSE_MATRIX_H_
#define ASSEMBLE_COARSE_MATRIX_H_
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <fstream>
#include <iostream>
#include <limits>
#include <locale>
#include <string>
using namespace dealii;
template <typename LevelNumber>
class CoarseDirectSolver
:
public MGCoarseGridBase<LinearAlgebra::distributed::Vector<LevelNumber>>{
public:
typedef LevelNumber number;
typedef LinearAlgebra::distributed::Vector<LevelNumber> VectorType;
CoarseDirectSolver(const IndexSet& owned_dofs_,
const IndexSet& relevant_dofs_);
template<class Operator>
void initialize(const Operator & vmult_operator,
const unsigned int& lvl=0);
template<class Operator, int dim>
void initialize(const Operator & vmult_operator,
const DoFHandler<dim> &dof_handler,
const unsigned int& lvl=0);
void vmult(VectorType &dst,
const VectorType &src) const{
inverse.solve(dst, src);
}
virtual void
operator()(const unsigned int lvl,
VectorType & dst,
const VectorType &src) const;
private:
int level;
IndexSet owned_dofs;
IndexSet relevant_dofs;
std::unique_ptr<TrilinosWrappers::SparseMatrix> matrix;
SolverControl solver_control;
mutable TrilinosWrappers::SolverDirect inverse;
template<int dim>
void setup_matrix(const DoFHandler<dim>& dof_handler_u);
template<class Operator>
void do_assembly(const Operator & vmult_operator);
};
template < typename LevelNumber>
CoarseDirectSolver<LevelNumber>::CoarseDirectSolver(const IndexSet& owned_dofs_,
const IndexSet& relevant_dofs_)
: level(-1)
, owned_dofs(owned_dofs_)
, relevant_dofs(relevant_dofs_)
, solver_control(10, 1e-14)
, inverse(solver_control)
{}
template < typename LevelNumber>
template<class Operator>
void
CoarseDirectSolver<LevelNumber>::initialize(const Operator & vmult_operator
, const unsigned int& lvl){
matrix=std::make_unique<TrilinosWrappers::SparseMatrix> (owned_dofs,
owned_dofs,
MPI_COMM_WORLD,
100);
level=lvl;
do_assembly( vmult_operator);
}
template < typename LevelNumber>
template<class Operator, int dim>
void
CoarseDirectSolver<LevelNumber>::initialize(const Operator & vmult_operator,
const DoFHandler<dim> &dof_handler,
const unsigned int& lvl){
level=lvl;
matrix=std::make_unique<TrilinosWrappers::SparseMatrix> ();
TrilinosWrappers::SparsityPattern A_sparsity( owned_dofs,
owned_dofs,
relevant_dofs,
MPI_COMM_WORLD,
100 );
typename DoFHandler<dim>::cell_iterator cell =
dof_handler.begin(level),
endc = dof_handler.end(level);
std::vector<types::global_dof_index> local_dof_indices_u(dof_handler.get_fe().dofs_per_cell);
for (; cell != endc; ++cell)
{
if(!cell->is_locally_owned_on_level())
continue;
cell->get_mg_dof_indices(local_dof_indices_u);
for(types::global_dof_index & i : local_dof_indices_u ){
A_sparsity.add_entries(i,
local_dof_indices_u.begin(),
local_dof_indices_u.end() );
}
}
A_sparsity.compress();
matrix->reinit(A_sparsity);
// std::cout<<"Sparisty pattern done"<<std::endl;
do_assembly( vmult_operator);
}
template < typename LevelNumber>
template<class Operator>
void
CoarseDirectSolver< LevelNumber>::do_assembly(const Operator & vmult_operator){
VectorType dst;
VectorType src;
// vmult_operator.get_matrix_free()->initialize_dof_vector(
// src, component);
src.reinit(owned_dofs, relevant_dofs, MPI_COMM_WORLD);
dst.reinit(owned_dofs, relevant_dofs, MPI_COMM_WORLD);
for(types::global_dof_index i =0; i< owned_dofs.size(); ++i){
src=0;
dst=0;
if(owned_dofs.is_element(i) )
src(i)=1;
src.compress(VectorOperation::insert);
vmult_operator.vmult(dst, src);
dst.update_ghost_values();
for(IndexSet::ElementIterator index_iter =owned_dofs.begin();
index_iter != owned_dofs.end();
++index_iter){
if(dst(*index_iter)!=0 )
matrix->set(*index_iter,i, dst(*index_iter) );
}
}
matrix->compress(VectorOperation::unknown);
inverse.initialize(*matrix);
// std::cout<<"Assembly&inverse done"<<std::endl;
}
template <typename LevelNumber>
void
CoarseDirectSolver< LevelNumber>::operator()(const unsigned int lvl,
VectorType & dst,
const VectorType &src) const {
Assert(lvl==level,ExcMessage("You are trying to use coarse solver on the wrong level"));
this->vmult(dst,src);
}
#endif /* ASSEMBLE_COARSE_MATRIX_H_ */
