DUNE-FEM (unstable)

novlpschwarz.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
6#define DUNE_ISTL_NOVLPSCHWARZ_HH
7
8#include <iostream> // for input/output to shell
9#include <fstream> // for input/output to files
10#include <vector> // STL vector class
11#include <sstream>
12
13#include <cmath> // Yes, we do some math here
14
15#include <dune/common/timer.hh>
16
17#include <dune/common/hybridutilities.hh>
18
19#include "io.hh"
20#include "bvector.hh"
21#include "vbvector.hh"
22#include "bcrsmatrix.hh"
23#include "io.hh"
24#include "gsetc.hh"
25#include "ilu.hh"
26#include "operators.hh"
27#include "solvers.hh"
28#include "preconditioners.hh"
29#include "scalarproducts.hh"
30#include "owneroverlapcopy.hh"
31
32namespace Dune {
33
59 template<class M, class X, class Y, class C>
61 {
62 public:
64 typedef M matrix_type;
66 typedef X domain_type;
68 typedef Y range_type;
70 typedef typename X::field_type field_type;
73
74 typedef typename C::PIS PIS;
75 typedef typename C::RI RI;
76 typedef typename RI::RemoteIndexList RIL;
77 typedef typename RI::const_iterator RIIterator;
78 typedef typename RIL::const_iterator RILIterator;
79 typedef typename M::ConstColIterator ColIterator;
80 typedef typename M::ConstRowIterator RowIterator;
81 typedef std::multimap<int,int> MM;
82 typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
83 typedef typename RIMap::iterator RIMapit;
84
93 : _A_(stackobject_to_shared_ptr(A)), communication(com), buildcomm(true)
94 {}
95
96 NonoverlappingSchwarzOperator (std::shared_ptr<const matrix_type> A, const communication_type& com)
97 : _A_(A), communication(com), buildcomm(true)
98 {}
99
101 void apply (const X& x, Y& y) const override
102 {
103 y = 0;
104 novlp_op_apply(x,y,1);
105 communication.addOwnerCopyToOwnerCopy(y,y);
106 }
107
109 void applyscaleadd (field_type alpha, const X& x, Y& y) const override
110 {
111 // only apply communication to alpha*A*x to make it consistent,
112 // y already has to be consistent.
113 Y y1(y);
114 y = 0;
115 novlp_op_apply(x,y,alpha);
116 communication.addOwnerCopyToOwnerCopy(y,y);
117 y += y1;
118 }
119
121 const M& getmat () const override
122 {
123 return *_A_;
124 }
125
126 void novlp_op_apply (const X& x, Y& y, field_type alpha) const
127 {
128 //get index sets
129 const PIS& pis=communication.indexSet();
130 const RI& ri = communication.remoteIndices();
131
132 // at the beginning make a multimap "bordercontribution".
133 // process has i and j as border dofs but is not the owner
134 // => only contribute to Ax if i,j is in bordercontribution
135 if (buildcomm) {
136
137 // set up mask vector
138 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
139 mask.resize(x.size());
140 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
141 mask[i] = 1;
142 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
143 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
144 mask[i->local().local()] = 0;
145 else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
146 mask[i->local().local()] = 2;
147 }
148
149 for (MM::iterator iter = bordercontribution.begin();
150 iter != bordercontribution.end(); ++iter)
151 bordercontribution.erase(iter);
152 std::map<int,int> owner; //key: local index i, value: process, that owns i
153 RIMap rimap;
154
155 // for each local index make multimap rimap:
156 // key: local index i, data: pair of process that knows i and pointer to RI entry
157 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i)
158 if (mask[i.index()] == 0)
159 for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
160 RIL& ril = *(remote->second.first);
161 for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
162 if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
163 if (rindex->localIndexPair().local().local() == i.index()) {
164 rimap.insert
165 (std::make_pair(i.index(),
166 std::pair<int,RILIterator>(remote->first, rindex)));
167 if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
168 owner.insert(std::make_pair(i.index(),remote->first));
169 }
170 }
171
172 int iowner = 0;
173 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
174 if (mask[i.index()] == 0) {
175 std::map<int,int>::iterator it = owner.find(i.index());
176 iowner = it->second;
177 std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
178 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
179 if (mask[j.index()] == 0) {
180 bool flag = true;
181 for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
182 std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
183 for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
184 if (foundj->second.first == foundi->second.first)
185 if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
186 || foundj->second.first == iowner
187 || foundj->second.first < communication.communicator().rank()) {
188 flag = false;
189 continue;
190 }
191 if (!flag)
192 continue;
193 }
194 // don´t contribute to Ax if
195 // 1. the owner of j has i as interior/border dof
196 // 2. iowner has j as interior/border dof
197 // 3. there is another process with smaller rank that has i and j
198 // as interor/border dofs
199 // if the owner of j does not have i as interior/border dof,
200 // it will not be taken into account
201 if (flag)
202 bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
203 }
204 }
205 }
206 }
207 buildcomm = false;
208 }
209
210 //compute alpha*A*x nonoverlapping case
211 for (RowIterator i = _A_->begin(); i != _A_->end(); ++i) {
212 if (mask[i.index()] == 0) {
213 //dof doesn't belong to process but is border (not ghost)
214 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j) {
215 if (mask[j.index()] == 1) //j is owner => then sum entries
216 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
217 else if (mask[j.index()] == 0) {
218 std::pair<MM::iterator, MM::iterator> itp =
219 bordercontribution.equal_range(i.index());
220 for (MM::iterator it = itp.first; it != itp.second; ++it)
221 if ((*it).second == (int)j.index())
222 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
223 }
224 }
225 }
226 else if (mask[i.index()] == 1) {
227 for (ColIterator j = (*_A_)[i.index()].begin(); j != (*_A_)[i.index()].end(); ++j)
228 if (mask[j.index()] != 2)
229 Impl::asMatrix(*j).usmv(alpha,x[j.index()],y[i.index()]);
230 }
231 }
232 }
233
236 {
238 }
239
242 {
243 return communication;
244 }
245 private:
246 std::shared_ptr<const matrix_type> _A_;
247 const communication_type& communication;
248 mutable bool buildcomm;
249 mutable std::vector<double> mask;
250 mutable std::multimap<int,int> bordercontribution;
251 };
252
255 namespace Amg
256 {
257 template<class T> struct ConstructionTraits;
258 }
259
274 template<class C, class P>
276 : public Preconditioner<typename P::domain_type,typename P::range_type> {
278 using X = typename P::domain_type;
279 using Y = typename P::range_type;
280 public:
282 typedef typename P::domain_type domain_type;
284 typedef typename P::range_type range_type;
287
296 : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
297 { }
298
306 NonoverlappingBlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
307 : _preconditioner(p), _communication(c)
308 { }
309
315 void pre (domain_type& x, range_type& b) override
316 {
317 _preconditioner->pre(x,b);
318 }
319
325 void apply (domain_type& v, const range_type& d) override
326 {
327 // block preconditioner equivalent to WrappedPreconditioner from
328 // pdelab/backend/ovlpistsolverbackend.hh,
329 // but not to BlockPreconditioner from schwarz.hh
330 _preconditioner->apply(v,d);
331 _communication.addOwnerCopyToOwnerCopy(v,v);
332 }
333
334 template<bool forward>
335 void apply (X& v, const Y& d)
336 {
337 _preconditioner->template apply<forward>(v,d);
338 _communication.addOwnerCopyToOwnerCopy(v,v);
339 }
340
346 void post (domain_type& x) override
347 {
348 _preconditioner->post(x);
349 }
350
353 {
355 }
356
357 private:
359 std::shared_ptr<P> _preconditioner;
360
362 const communication_type& _communication;
363 };
364
367} // end namespace
368
369#endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A linear operator exporting itself in matrix form.
Definition: operators.hh:110
Nonoverlapping parallel preconditioner.
Definition: novlpschwarz.hh:276
NonoverlappingBlockPreconditioner(P &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:295
void post(domain_type &x) override
Clean up.
Definition: novlpschwarz.hh:346
P::range_type range_type
The range type of the preconditioner.
Definition: novlpschwarz.hh:284
void apply(domain_type &v, const range_type &d) override
Apply the preconditioner.
Definition: novlpschwarz.hh:325
NonoverlappingBlockPreconditioner(const std::shared_ptr< P > &p, const communication_type &c)
Constructor.
Definition: novlpschwarz.hh:306
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:286
void pre(domain_type &x, range_type &b) override
Prepare the preconditioner.
Definition: novlpschwarz.hh:315
void apply(X &v, const Y &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: novlpschwarz.hh:335
SolverCategory::Category category() const override
Category of the preconditioner (see SolverCategory::Category)
Definition: novlpschwarz.hh:352
P::domain_type domain_type
The domain type of the preconditioner.
Definition: novlpschwarz.hh:282
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:61
C communication_type
The type of the communication object.
Definition: novlpschwarz.hh:72
SolverCategory::Category category() const override
Category of the linear operator (see SolverCategory::Category)
Definition: novlpschwarz.hh:235
X domain_type
The type of the domain.
Definition: novlpschwarz.hh:66
const M & getmat() const override
get reference to matrix
Definition: novlpschwarz.hh:121
Y range_type
The type of the range.
Definition: novlpschwarz.hh:68
M matrix_type
The type of the matrix we operate on.
Definition: novlpschwarz.hh:64
const communication_type & getCommunication() const
Get the object responsible for communication.
Definition: novlpschwarz.hh:241
X::field_type field_type
The field type of the range.
Definition: novlpschwarz.hh:70
NonoverlappingSchwarzOperator(const matrix_type &A, const communication_type &com)
constructor: just store a reference to a matrix.
Definition: novlpschwarz.hh:92
void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition: novlpschwarz.hh:109
void apply(const X &x, Y &y) const override
apply operator to x:
Definition: novlpschwarz.hh:101
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Some generic functions for pretty printing vectors and matrices.
Classes providing communication interfaces for overlapping Schwarz methods.
Define base class for scalar product and norm.
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
The incomplete LU factorization kernels.
Dune namespace.
Definition: alignedallocator.hh:13
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:72
Define general, extensible interface for operators. The available implementation wraps a matrix.
Define general preconditioner interface.
Implementations of the inverse operator interface.
Traits class for generically constructing non default constructable types.
Definition: construction.hh:39
Category
Definition: solvercategory.hh:23
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:27
A simple timing class.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)