Loading...
Searching...
No Matches
novlpschwarz.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
#define DUNE_ISTL_NOVLPSCHWARZ_HH
#include <iostream> // for input/output to shell
#include <fstream> // for input/output to files
#include <vector> // STL vector class
#include <sstream>
#include <cmath> // Yes, we do some math here
#include <dune/common/timer.hh>
#include <dune/common/hybridutilities.hh>
#include "io.hh"
#include "bvector.hh"
#include "vbvector.hh"
#include "bcrsmatrix.hh"
#include "io.hh"
#include "gsetc.hh"
#include "ilu.hh"
#include "operators.hh"
#include "solvers.hh"
#include "preconditioners.hh"
#include "scalarproducts.hh"
#include "owneroverlapcopy.hh"
namespace Dune {
template<class M, class X, class Y, class C>
class NonoverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
{
public:
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
typedef C communication_type;
typedef typename C::PIS PIS;
typedef typename C::RI RI;
typedef typename RI::RemoteIndexList RIL;
typedef typename RI::const_iterator RIIterator;
typedef typename RIL::const_iterator RILIterator;
typedef typename M::ConstColIterator ColIterator;
typedef typename M::ConstRowIterator RowIterator;
typedef std::multimap<int,int> MM;
typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
typedef typename RIMap::iterator RIMapit;
NonoverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
: _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
{}
NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
: _A_(A), communication(com), buildcomm(true)
{}
void apply (const X& x, Y& y) const override
{
y = 0;
novlp_op_apply(x,y,1);
communication.addOwnerCopyToOwnerCopy(y,y);
}
void applyscaleadd (field_type alpha, const X& x, Y& y) const override
{
// only apply communication to alpha*A*x to make it consistent,
// y already has to be consistent.
Y y1(y);
y = 0;
novlp_op_apply(x,y,alpha);
communication.addOwnerCopyToOwnerCopy(y,y);
y += y1;
}
const M& getmat () const override
{
return *_A_;
}
void novlp_op_apply (const X& x, Y& y, field_type alpha) const
{
//get index sets
const PIS& pis=communication.indexSet();
const RI& ri = communication.remoteIndices();
// at the beginning make a multimap "bordercontribution".
// process has i and j as border dofs but is not the owner
// => only contribute to Ax if i,j is in bordercontribution
if (buildcomm) {
// set up mask vector
if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
mask.resize(x.size());
for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
mask[i] = 1;
for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
mask[i->local().local()] = 0;
else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
mask[i->local().local()] = 2;
}
for (MM::iterator iter = bordercontribution.begin();
iter != bordercontribution.end(); ++iter)
bordercontribution.erase(iter);
std::map<int,int> owner; //key: local index i, value: process, that owns i
RIMap rimap;
// for each local index make multimap rimap:
// key: local index i, data: pair of process that knows i and pointer to RI entry
for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
if (mask[i.index()] == 0)
for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
RIL& ril = *(remote->second.first);
for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
if (rindex->localIndexPair().local().local() == i.index()) {
rimap.insert
(std::make_pair(i.index(),
std::pair<int,RILIterator>(remote->first, rindex)));
if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
owner.insert(std::make_pair(i.index(),remote->first));
}
}
int iowner = 0;
for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
if (mask[i.index()] == 0) {
std::map<int,int>::iterator it = owner.find(i.index());
iowner = it->second;
std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
if (mask[j.index()] == 0) {
bool flag = true;
for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
if (foundj->second.first == foundi->second.first)
if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
|| foundj->second.first == iowner
|| foundj->second.first < communication.communicator().rank()) {
flag = false;
continue;
}
if (!flag)
continue;
}
// don´t contribute to Ax if
// 1. the owner of j has i as interior/border dof
// 2. iowner has j as interior/border dof
// 3. there is another process with smaller rank that has i and j
// as interor/border dofs
// if the owner of j does not have i as interior/border dof,
// it will not be taken into account
if (flag)
bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
}
}
}
}
buildcomm = false;
}
//compute alpha*A*x nonoverlapping case
for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
if (mask[i.index()] == 0) {
//dof doesn't belong to process but is border (not ghost)
for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
if (mask[j.index()] == 1) //j is owner => then sum entries
Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
else if (mask[j.index()] == 0) {
std::pair<MM::iterator, MM::iterator> itp =
bordercontribution.equal_range(i.index());
for (MM::iterator it = itp.first; it != itp.second; ++it)
if ((*it).second == (int)j.index())
Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
}
}
}
else if (mask[i.index()] == 1) {
for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j)
if (mask[j.index()] != 2)
Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
}
}
}
SolverCategory::Category category() const override
{
return SolverCategory::nonoverlapping;
}
const communication_type& getCommunication() const
{
return communication;
}
private:
std::shared_ptr<const matrix_type> _A_;
const communication_type& communication;
mutable bool buildcomm;
mutable std::vector<double> mask;
mutable std::multimap<int,int> bordercontribution;
};
namespace Amg
{
template<class T> struct ConstructionTraits;
}
template<class C, class P>
class NonoverlappingBlockPreconditioner
: public Preconditioner<typename P::domain_type,typename P::range_type> {
friend struct Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
using X = typename P::domain_type;
using Y = typename P::range_type;
public:
typedef typename P::domain_type domain_type;
typedef typename P::range_type range_type;
typedef C communication_type;
NonoverlappingBlockPreconditioner (P& p, const communication_type& c)
: _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
{ }
NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
: _preconditioner(p), _communication(c)
{ }
void pre (domain_type& x, range_type& b) override
{
_preconditioner->pre(x,b);
}
void apply (domain_type& v, const range_type& d) override
{
// block preconditioner equivalent to WrappedPreconditioner from
// pdelab/backend/ovlpistsolverbackend.hh,
// but not to BlockPreconditioner from schwarz.hh
_preconditioner->apply(v,d);
_communication.addOwnerCopyToOwnerCopy(v,v);
}
template<bool forward>
void apply (X& v, const Y& d)
{
_preconditioner->template apply<forward>(v,d);
_communication.addOwnerCopyToOwnerCopy(v,v);
}
void post (domain_type& x) override
{
_preconditioner->post(x);
}
SolverCategory::Category category() const override
{
return SolverCategory::nonoverlapping;
}
private:
std::shared_ptr<P> _preconditioner;
const communication_type& _communication;
};
} // end namespace
#endif
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
gsetc.hh
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
owneroverlapcopy.hh
Classes providing communication interfaces for overlapping Schwarz methods.
preconditioners.hh
Define general preconditioner interface.
solvers.hh
Implementations of the inverse operator interface.
io.hh
Some generic functions for pretty printing vectors and matrices.
operators.hh
Define general, extensible interface for operators. The available implementation wraps a matrix.
bvector.hh
This file implements a vector space as a tensor product of a given vector space. The number of compon...
scalarproducts.hh
Define base class for scalar product and norm.
ilu.hh
The incomplete LU factorization kernels.
vbvector.hh
???
Dune
Definition allocator.hh:11
Dune::Amg
Definition novlpschwarz.hh:256
Dune::NonoverlappingSchwarzOperator::PIS
C::PIS PIS
Definition novlpschwarz.hh:74
Dune::NonoverlappingSchwarzOperator::communication_type
C communication_type
The type of the communication object.
Definition novlpschwarz.hh:72
Dune::NonoverlappingSchwarzOperator::RIMap
std::multimap< int, std::pair< int, RILIterator > > RIMap
Definition novlpschwarz.hh:82
Dune::NonoverlappingSchwarzOperator::RI
C::RI RI
Definition novlpschwarz.hh:75
Dune::NonoverlappingSchwarzOperator::novlp_op_apply
void novlp_op_apply(const X &x, Y &y, field_type alpha) const
Definition novlpschwarz.hh:126
Dune::NonoverlappingSchwarzOperator::NonoverlappingSchwarzOperator
NonoverlappingSchwarzOperator(std::shared_ptr< const matrix_type > A, const communication_type &com)
Definition novlpschwarz.hh:96
Dune::NonoverlappingSchwarzOperator::category
SolverCategory::Category category() const override
Category of the linear operator (see SolverCategory::Category).
Definition novlpschwarz.hh:235
Dune::NonoverlappingSchwarzOperator::ColIterator
M::ConstColIterator ColIterator
Definition novlpschwarz.hh:79
Dune::NonoverlappingSchwarzOperator::domain_type
X domain_type
The type of the domain.
Definition novlpschwarz.hh:66
Dune::NonoverlappingSchwarzOperator::getmat
const M & getmat() const override
get reference to matrix
Definition novlpschwarz.hh:121
Dune::NonoverlappingSchwarzOperator::RIMapit
RIMap::iterator RIMapit
Definition novlpschwarz.hh:83
Dune::NonoverlappingSchwarzOperator::RILIterator
RIL::const_iterator RILIterator
Definition novlpschwarz.hh:78
Dune::NonoverlappingSchwarzOperator::MM
std::multimap< int, int > MM
Definition novlpschwarz.hh:81
Dune::NonoverlappingSchwarzOperator::range_type
Y range_type
The type of the range.
Definition novlpschwarz.hh:68
Dune::NonoverlappingSchwarzOperator::matrix_type
M matrix_type
The type of the matrix we operate on.
Definition novlpschwarz.hh:64
Dune::NonoverlappingSchwarzOperator::RowIterator
M::ConstRowIterator RowIterator
Definition novlpschwarz.hh:80
Dune::NonoverlappingSchwarzOperator::getCommunication
const communication_type & getCommunication() const
Get the object responsible for communication.
Definition novlpschwarz.hh:241
Dune::NonoverlappingSchwarzOperator::RIL
RI::RemoteIndexList RIL
Definition novlpschwarz.hh:76
Dune::NonoverlappingSchwarzOperator::field_type
X::field_type field_type
The field type of the range.
Definition novlpschwarz.hh:70
Dune::NonoverlappingSchwarzOperator::NonoverlappingSchwarzOperator
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition novlpschwarz.hh:92
Dune::NonoverlappingSchwarzOperator::applyscaleadd
void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition novlpschwarz.hh:109
Dune::NonoverlappingSchwarzOperator::apply
void apply(const X &x, Y &y) const override
apply operator to x:
Definition novlpschwarz.hh:101
Dune::NonoverlappingSchwarzOperator::RIIterator
RI::const_iterator RIIterator
Definition novlpschwarz.hh:77
Dune::Amg::ConstructionTraits
Traits class for generically constructing non default constructable types.
Definition construction.hh:39
Dune::NonoverlappingBlockPreconditioner::NonoverlappingBlockPreconditioner
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition novlpschwarz.hh:295
Dune::NonoverlappingBlockPreconditioner::post
void post(domain_type &x) override
Clean up.
Definition novlpschwarz.hh:346
Dune::NonoverlappingBlockPreconditioner::range_type
P::range_type range_type
The range type of the preconditioner.
Definition novlpschwarz.hh:284
Dune::NonoverlappingBlockPreconditioner::apply
void apply(domain_type &v, const range_type &d) override
Apply the preconditioner.
Definition novlpschwarz.hh:325
Dune::NonoverlappingBlockPreconditioner::NonoverlappingBlockPreconditioner
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition novlpschwarz.hh:306
Dune::NonoverlappingBlockPreconditioner::communication_type
C communication_type
The type of the communication object.
Definition novlpschwarz.hh:286
Dune::NonoverlappingBlockPreconditioner::pre
void pre(domain_type &x, range_type &b) override
Prepare the preconditioner.
Definition novlpschwarz.hh:315
Dune::NonoverlappingBlockPreconditioner::apply
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition novlpschwarz.hh:335
Dune::NonoverlappingBlockPreconditioner::category
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category).
Definition novlpschwarz.hh:352
Dune::NonoverlappingBlockPreconditioner::domain_type
P::domain_type domain_type
The domain type of the preconditioner.
Definition novlpschwarz.hh:282
Dune::AssembledLinearOperator
A linear operator exporting itself in matrix form.
Definition operators.hh:110
Dune::OwnerOverlapCopyAttributeSet::owner
@ owner
Definition owneroverlapcopy.hh:61
Dune::OwnerOverlapCopyAttributeSet::copy
@ copy
Definition owneroverlapcopy.hh:61
Dune::OwnerOverlapCopyAttributeSet::overlap
@ overlap
Definition owneroverlapcopy.hh:61
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:33
Dune::SolverCategory::Category
Category
Definition solvercategory.hh:23
Dune::SolverCategory::nonoverlapping
@ nonoverlapping
Category for non-overlapping solvers.
Definition solvercategory.hh:27