DUNE PDELab (git)

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
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.
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.
Classes providing communication interfaces for overlapping Schwarz methods.
Define general preconditioner interface.
Define base class for scalar product and norm.
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)