amgparallelhelpers.hh 36.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
Thomas Fetzer's avatar
Thomas Fetzer committed
19
20
/*!
 * \file
21
 * \ingroup Linear
Thomas Fetzer's avatar
Thomas Fetzer committed
22
 * \brief Provides a helper class for nonoverlapping
23
 *        decomposition using the ISTL AMG.
Thomas Fetzer's avatar
Thomas Fetzer committed
24
 */
25
26
27
#ifndef DUMUX_AMGPARALLELHELPERS_HH
#define DUMUX_AMGPARALLELHELPERS_HH

Timo Koch's avatar
Timo Koch committed
28
29
#include <dune/geometry/dimension.hh>
#include <dune/grid/common/datahandleif.hh>
Timo Koch's avatar
Timo Koch committed
30
#include <dune/grid/common/partitionset.hh>
31
32
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
Timo Koch's avatar
Timo Koch committed
33
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
34
#include <dumux/linear/amgtraits.hh>
35

36
namespace Dumux {
37

Thomas Fetzer's avatar
Thomas Fetzer committed
38
/*!
39
 * \ingroup Linear
Thomas Fetzer's avatar
Thomas Fetzer committed
40
41
42
 * \brief A parallel helper class providing a nonoverlapping
 *        decomposition of all degrees of freedom
 */
43
44
45
46
// operator that resets result to zero at constrained DOFS
template<class TypeTag>
class ParallelISTLHelper
{
47
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
48
    using AmgTraits = Dumux::AmgTraits<TypeTag>;
49
    using DofMapper = typename AmgTraits::DofMapper;
50

51
    enum { dofCodim = AmgTraits::dofCodim };
52
53
54
55

    class BaseGatherScatter
    {
    public:
56
57
        BaseGatherScatter(const DofMapper& mapper)
        : mapper_(mapper) {}
58
59

        template<class EntityType>
60
        int index(const EntityType& e) const
61
        { return mapper_.index(e); }
62
63

    private:
64
        const DofMapper& mapper_;
65
66
    };

67
68
    /*!
     * \brief GatherScatter implementation that makes a right hand side in the box model consistent.
69
70
71
72
73
74
75
     */
    template<class V>
    class ConsistencyBoxGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<ConsistencyBoxGatherScatter<V>,typename V::block_type>
    {
    public:
76
        using DataType = typename V::block_type;
77

78
79
        ConsistencyBoxGatherScatter(V& container, const DofMapper& mapper)
        : BaseGatherScatter(mapper), container_(container)
80
        {}
81

82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

        bool fixedsize(int dim, int codim) const
        {
            return true;
        }

        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
        }

        template<class MessageBuffer, class EntityType>
        void gather (MessageBuffer& buff, const EntityType& e) const
        {
101
            buff.write(container_[this->index(e)]);
102
103
104
105
106
107
108
        }

        template<class MessageBuffer, class EntityType>
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
        {
            typename V::block_type block;
            buff.read(block);
109
            container_[this->index(e)]+=block;
110
111
112
113
        }
    private:
        V& container_;
    };
114

115
116

    /**
117
     * \brief Writes 1<<24 to each data item (of the container) that is gathered or scattered
118
119
120
121
122
     * and is neither interior nor border.
     *
     * Can be used to mark ghost cells.
     */
    class GhostGatherScatter
123
124
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<GhostGatherScatter,std::size_t>
125
126
    {
    public:
127
        using DataType = std::size_t;
128

129
130
        GhostGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
        : BaseGatherScatter(mapper), ranks_(ranks)
131
132
        {}

133
134
135
136
137
138

        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

139
        bool fixedsize(int dim, int codim) const
140
141
        {
            return true;
142
        }
143

144
145
146
147
        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
148
149
        }

150
        template<class MessageBuffer, class EntityType>
151
        void gather (MessageBuffer& buff, const EntityType& e) const
152
        {
153
            std::size_t& data= ranks_[this->index(e)];
154
155
156
157
158
159
            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
                data = (1<<24);
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
160
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
161
162
        {
            std::size_t x;
163
            std::size_t& data = ranks_[this->index(e)];
164
165
166
167
168
169
170
171
172
            buff.read(x);
            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
                data= (1<<24);
        }
    private:
        std::vector<std::size_t>& ranks_;
    };

    /**
173
     * \brief GatherScatter handle that sets 1<<24 for data items neither associated to
174
175
176
177
178
     * the interior or border and take the minimum when scattering.
     *
     * Used to compute an owner rank for each unknown.
     */
    class InteriorBorderGatherScatter
179
        : public BaseGatherScatter,
180
          public Dune::CommDataHandleIF<InteriorBorderGatherScatter,std::size_t>
181
182
    {
    public:
183
        using DataType = std::size_t;
184

185
186
        InteriorBorderGatherScatter(std::vector<std::size_t>& ranks, const DofMapper& mapper)
        : BaseGatherScatter(mapper), ranks_(ranks)
187
188
        {}

189
190
191
192
193
194

        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

195
        bool fixedsize(int dim, int codim) const
196
197
198
199
        {
            return true;

        }
200
201
202
203
204
205
206

        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
        }

207
        template<class MessageBuffer, class EntityType>
208
        void gather (MessageBuffer& buff, const EntityType& e) const
209
210
        {

211
            std::size_t& data = ranks_[this->index(e)];
212
213
214
215
216
217
            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
                data = (1<<24);
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
218
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
219
        {
220
            using std::min;
221
            std::size_t x;
222
            std::size_t& data = ranks_[this->index(e)];
223
224
225
226
            buff.read(x);
            if (e.partitionType()!=Dune::InteriorEntity && e.partitionType()!=Dune::BorderEntity)
                data = x;
            else
227
                data = min(data,x);
228
229
230
231
232
233
        }
    private:
        std::vector<std::size_t>& ranks_;
    };

    /**
234
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
235
236
237
     *
     */
    struct NeighbourGatherScatter
238
239
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<NeighbourGatherScatter,int>
240
    {
241
        using DataType = int;
242

243
244
        NeighbourGatherScatter(const DofMapper& mapper, int rank_, std::set<int>& neighbours_)
        : BaseGatherScatter(mapper), myrank(rank_), neighbours(neighbours_)
245
246
        {}

247
248
249
250
251
252

        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

253
        bool fixedsize(int dim, int codim) const
254
255
        {
            return true;
256
        }
257

258
259
260
261
        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
262
263
        }

264
        template<class MessageBuffer, class EntityType>
265
        void gather (MessageBuffer& buff, const EntityType &e) const
266
267
268
269
270
        {
            buff.write(myrank);
        }

        template<class MessageBuffer, class EntityType>
271
        void scatter (MessageBuffer& buff, const EntityType &e, size_t n)
272
273
274
275
276
277
278
279
280
281
282
        {
            int x;
            buff.read(x);
            neighbours.insert(x);
        }
        int myrank;
        std::set<int>& neighbours;
    };


    /**
283
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
284
285
286
     *
     */
    struct SharedGatherScatter
287
288
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<SharedGatherScatter,int>
289
    {
290
        using DataType = int;
291

292
293
        SharedGatherScatter(std::vector<int>& shared, const DofMapper& mapper)
        : BaseGatherScatter(mapper), shared_(shared)
294
295
        {}

296
297
298
299
300
        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

301
        bool fixedsize(int dim, int codim) const
302
303
304
305
306
        {
            return true;

        }

307
308
309
310
311
312
        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
        }

313
        template<class MessageBuffer, class EntityType>
314
        void gather (MessageBuffer& buff, EntityType& e) const
315
        {
316
            int data=true;
317
318
319
320
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
321
        void scatter (MessageBuffer& buff, const EntityType &e, size_t n)
322
        {
323
            int x;
324
            buff.read(x);
325
            int& data= shared_[this->index(e)];
326
327
328
329
330
331
332
            data = data || x;
        }
    private:
        std::vector<int>& shared_;
    };

    /**
333
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
334
335
336
337
     *
     */
    template<typename GI>
    struct GlobalIndexGatherScatter
338
339
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<GlobalIndexGatherScatter<GI>, GI>
340
    {
341
        using DataType = GI;
342
343
        GlobalIndexGatherScatter(std::vector<GI>& gindices, const DofMapper& mapper)
        : BaseGatherScatter(mapper), gindices_(gindices)
344
345
        {}

346
347
348
349
350
        bool contains(int dim, int codim) const
        {
            return dofCodim==codim;
        }

351
        bool fixedsize(int dim, int codim) const
352
353
        {
            return true;
354
        }
355

356
357
358
359
        template<class EntityType>
        size_t size (EntityType& e) const
        {
            return 1;
360
361
        }

362
        template<class MessageBuffer, class EntityType>
363
        void gather (MessageBuffer& buff, const EntityType& e) const
364
        {
365
            buff.write(gindices_[this->index(e)]);
366
367
368
        }

        template<class MessageBuffer, class EntityType>
369
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
370
        {
371
            using std::min;
372
373
            DataType x;
            buff.read(x);
374
            gindices_[this->index(e)] = min(gindices_[this->index(e)], x);
375
376
377
378
379
380
381
        }
    private:
        std::vector<GI>& gindices_;
    };

public:

382
383
    ParallelISTLHelper (const GridView& gridView, const DofMapper& mapper, int verbose=1)
        : gridView_(gridView), mapper_(mapper), verbose_(verbose), initialized_(false)
384
    {}
385

386
387
    // \brief Initializes the markers for ghosts and owners with the correct size and values.
    //
388
389
390
391
392
    void initGhostsAndOwners()
    {
        owner_.resize(mapper_.size(),
                      gridView_.comm().rank());
        isGhost_.resize(mapper_.size(),0.0);
393
        // find out about ghosts
394
        GhostGatherScatter ggs(owner_, mapper_);
395

396
397
        if (gridView_.comm().size()>1)
            gridView_.communicate(ggs, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);
398
399
400
401

        isGhost_ = owner_;

        // partition interior/border
402
        InteriorBorderGatherScatter dh(owner_, mapper_);
403

404
405
        if (gridView_.comm().size()>1)
            gridView_.communicate(dh, Dune::InteriorBorder_InteriorBorder_Interface, Dune::ForwardCommunication);
406
407

        // convert vector into mask vector
408
409
        for (auto& v : owner_)
            v = (v == gridView_.comm().rank()) ? 1.0 : 0.0;
410
411

        initialized_=true;
412
413
414
415
416
417
418
419
    }

    // keep only DOFs assigned to this processor
    template<typename W>
    void mask (W& w) const
    {
        auto v1=w.begin();

420
        for(auto v2=owner_.begin(), vend=owner_.end(); v2!=vend; ++v1, ++v2)
421
422
423
424
425
426
427
428
429
430
431
432
433
434
            v1*=v2;
    }

    // access to mask vector
    double mask (std::size_t i) const
    {
        return owner_[i];
    }

    // access to ghost vector
    double ghost (std::size_t i) const
    {
        return isGhost_[i];
    }
435

436
437
438
439
440
    // \brief Make a vector of the box model consistent.
    template<typename B, typename A>
    void makeNonOverlappingConsistent(Dune::BlockVector<B,A>& v)
    {
#if HAVE_MPI
441
442
443
444
        ConsistencyBoxGatherScatter<Dune::BlockVector<B,A> > gs(v, mapper_);
        if (gridView_.comm().size() > 1)
            gridView_.communicate(gs, Dune::InteriorBorder_InteriorBorder_Interface,
                                  Dune::ForwardCommunication);
445
446
#endif
    }
447

448
449
450
451

#if HAVE_MPI

    /**
452
     * \brief Creates a matrix suitable for parallel AMG and the parallel information
453
454
     *
     *
455
456
457
458
     * \tparam MatrixType The type of the ISTL matrix used.
     * \tparam Comm The type of the OwnerOverlapCopyCommunication
     * \param m The local matrix.
     * \param c The parallel information object providing index set, interfaces and
459
460
461
462
463
     * communicators.
     */
    template<typename MatrixType, typename Comm>
    void createIndexSetAndProjectForAMG(MatrixType& m, Comm& c);
#endif
464

465
private:
466
467
468
469
470
471
    const GridView gridView_; //!< the grid view
    const DofMapper& mapper_; //!< the dof mapper
    std::vector<std::size_t> owner_; //!< vector to identify unique decomposition
    std::vector<std::size_t> isGhost_; //!< vector to identify ghost dofs
    int verbose_; //!< verbosity
    bool initialized_; //!< whether isGhost and owner arrays are initialized
472
473

}; // class ParallelISTLHelper
474

475
476
477
478
479
/*!
 * \ingroup Linear
 * \brief Helper class for adding up matrix entries on border.
 * \tparam GridOperator The grid operator to work on.
 * \tparam MatrixType The MatrixType.
480
481
482
483
 */
template<class TypeTag>
class EntityExchanger
{
484
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
485
    using AmgTraits = Dumux::AmgTraits<TypeTag>;
486
    enum { numEq = AmgTraits::numEq };
487
488
    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,numEq,numEq> >;
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
489
490
    enum { dim = GridView::dimension };
    enum { dofCodim = AmgTraits::dofCodim };
491
492
493
494
495
496
497
    using Grid = typename GridView::Traits::Grid;
    using BlockType = typename Matrix::block_type;
    using IDS = typename Grid::Traits::GlobalIdSet;
    using IdType = typename IDS::IdType;
    using RowIterator = typename Matrix::RowIterator;
    using ColIterator = typename Matrix::ColIterator;
    using DofMapper = typename AmgTraits::DofMapper;
498
499
500

public:
    /*! \brief Constructor. Sets up the local to global relations.
501
      \param[in] gridView The gridView on which we are operating
502
      \param[in] mapper The local dof mapper
503
    */
504
505
    EntityExchanger(const GridView& gridView, const DofMapper& mapper)
    : gridView_(gridView), mapper_(mapper)
506
507
508
509
    {
        gid2Index_.clear();
        index2GID_.clear();

510
        for (const auto& entity : entities(gridView_, Dune::Codim<dofCodim>()))
511
        {
512
            if (entity.partitionType() == Dune::BorderEntity)
513
            {
514
515
                int localIdx = mapper_.index(entity);
                IdType dofIdxGlobal = gridView_.grid().globalIdSet().id(entity);
516

517
                std::pair<IdType, int> g2iPair(dofIdxGlobal, localIdx);
518
519
                gid2Index_.insert(g2iPair);

520
                std::pair<int, IdType> i2gPair(localIdx, dofIdxGlobal);
521
522
523
524
525
526
                index2GID_.insert(i2gPair);
            }
        }
    }

    /**
527
     * \brief A DataHandle class to exchange matrix sparsity patterns.
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
     *
     *  We look at a 2D example with a nonoverlapping grid,
     *  two processes and no ghosts with Q1 discretization.
     *  Process 0 has the left part of the domain
     *  with three cells and eight vertices (1-8),
     *  Process 1 the right part with three cells
     *  and eight vertices (2,4,7-12).
     *  <pre>
     *  1 _ 2        2 _ 9 _ 10
     *  |   |        |   |   |
     *  3 _ 4 _ 7    4 _ 7 _ 11
     *  |   |   |        |   |
     *  5 _ 6 _ 8        8 _ 12
     *  </pre>
     *  If we look at vertex 7 and the corresponding entries in the matrix for P0,
     *  there will be entries for (7,4) and (7,8), but not for (7,2).
     *  The MatPatternExchange class will find these entries and returns a vector "sparsity",
     *  that contains all missing connections.
     */
    class MatPatternExchange
548
549
550
551
        : public Dune::CommDataHandleIF<MatPatternExchange, IdType>
    {
        using RowIterator = typename Matrix::RowIterator;
        using ColIterator = typename Matrix::ColIterator;
552
553
    public:
        //! Export type of data for message buffer
554
555
        using DataType = IdType;

556
557
558
559
560
561
        /** \brief Constructor
            \param[in] mapper The local dof mapper.
            \param[in] g2i Global to local index map.
            \param[in] i2g Local to global index map.
            \param[in] A Matrix to operate on.
            \param[in] helper parallel istl helper.
562
563
564
565
566
567
568
569
        */
        MatPatternExchange (const DofMapper& mapper,
                            const std::map<IdType,int>& g2i,
                            const std::map<int,IdType>& i2g, Matrix& A,
                            const ParallelISTLHelper<TypeTag>& helper)
            : mapper_(mapper), gid2Index_(g2i), index2GID_(i2g),
              sparsity_(A.N()), A_(A), helper_(helper)
        {}
570

571
        /** \brief Returns true if data for given valid codim should be communicated
572
573
574
575
576
577
         */
        bool contains (int dim, int codim) const
        {
            return (codim==dofCodim);
        }

578
        /** \brief Returns true if size of data per entity of given dim and codim is a constant
579
580
581
582
583
584
         */
        bool fixedsize (int dim, int codim) const
        {
            return false;
        }

585
        /** \brief How many objects of type DataType have to be sent for a given entity
586
587
588
589
         */
        template<class EntityType>
        size_t size (EntityType& e) const
        {
590
            int i = mapper_.index(e);
591
592
593
594
595
596
597
598
599
600
601
            int n = 0;
            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
            {
                typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
                if (it != index2GID_.end())
                    n++;
            }

            return n;
        }

602
        /** \brief Pack data from user to message buffer
603
604
605
606
         */
        template<class MessageBuffer, class EntityType>
        void gather (MessageBuffer& buff, const EntityType& e) const
        {
607
            int i = mapper_.index(e);
608
609
610
611
612
613
614
615
616
            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
            {
                typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
                if (it != index2GID_.end())
                    buff.write(it->second);
            }

        }

617
        /** \brief Unpack data from message buffer to user
618
619
620
621
         */
        template<class MessageBuffer, class EntityType>
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
        {
622
            int i = mapper_.index(e);
623
624
625
626
627
628
629
630
631
632
633
634
635
636
            for (size_t k = 0; k < n; k++)
            {
                IdType id;
                buff.read(id);
                // only add entries corresponding to border entities
                typename std::map<IdType,int>::const_iterator it = gid2Index_.find(id);
                if (it != gid2Index_.end()
                    && sparsity_[i].find(it->second) == sparsity_[i].end()
                    && helper_.ghost(it->second) != 1<<24)
                    sparsity_[i].insert(it->second);
            }
        }

        /**
637
         * \brief Get the communicated sparsity pattern
638
639
640
641
642
643
644
645
         * @return the vector with the sparsity pattern
         */
        std::vector<std::set<int> >& sparsity ()
        {
            return sparsity_;
        }

    private:
646
        const DofMapper& mapper_;
647
648
649
650
651
        const std::map<IdType,int>& gid2Index_;
        const std::map<int,IdType>& index2GID_;
        std::vector<std::set<int> > sparsity_;
        Matrix& A_;
        const ParallelISTLHelper<TypeTag>& helper_;
652
653

    }; // class MatPatternExchange
654
655
656
657
658
659

    //! Local matrix blocks associated with the global id set
    struct MatEntry
    {
        IdType first;
        BlockType second;
660
        MatEntry (const IdType& f, const BlockType& s) : first(f), second(s) {}
661
662
663
664
665
        MatEntry () {}
    };

    //! A DataHandle class to exchange matrix entries
    class MatEntryExchange
666
667
668
669
        : public Dune::CommDataHandleIF<MatEntryExchange,MatEntry>
    {
        using RowIterator = typename Matrix::RowIterator;
        using ColIterator = typename Matrix::ColIterator;
670
671
    public:
        //! Export type of data for message buffer
672
673
        using DataType = MatEntry;

674
675
676
677
678
        /** \brief Constructor
            \param[in] mapper The local dof mapper.
            \param[in] g2i Global to local index map.
            \param[in] i2g Local to global index map.
            \param[in] A Matrix to operate on.
679
680
681
682
683
684
        */
        MatEntryExchange (const DofMapper& mapper, const std::map<IdType,int>& g2i,
                          const std::map<int,IdType>& i2g,
                          Matrix& A)
            : mapper_(mapper), gid2Index_(g2i), index2GID_(i2g), A_(A)
        {}
685

686
        /** \brief Returns true if data for given valid codim should be communicated
687
688
689
690
691
692
         */
        bool contains (int dim, int codim) const
        {
            return (codim==dofCodim);
        }

693
        /** \brief Returns true if size of data per entity of given dim and codim is a constant
694
695
696
697
698
699
         */
        bool fixedsize (int dim, int codim) const
        {
            return false;
        }

700
        /** \brief How many objects of type DataType have to be sent for a given entity
701
702
703
704
         */
        template<class EntityType>
        size_t size (EntityType& e) const
        {
705
            int i = mapper_.index(e);
706
707
708
709
710
711
712
713
714
715
716
            int n = 0;
            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
            {
                typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
                if (it != index2GID_.end())
                    n++;
            }

            return n;
        }

717
        /** \brief Pack data from user to message buffer
718
719
720
721
         */
        template<class MessageBuffer, class EntityType>
        void gather (MessageBuffer& buff, const EntityType& e) const
        {
722
            int i = mapper_.index(e);
723
724
725
726
727
728
729
730
731
            for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
            {
                typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
                if (it != index2GID_.end())
                    buff.write(MatEntry(it->second,*j));
            }

        }

732
        /** \brief Unpack data from message buffer to user
733
734
735
736
         */
        template<class MessageBuffer, class EntityType>
        void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
        {
737
            int i = mapper_.index(e);
738
739
740
741
742
743
744
745
746
747
748
749
750
            for (size_t k = 0; k < n; k++)
            {
                MatEntry m;
                buff.read(m);
                // only add entries corresponding to border entities
                typename std::map<IdType,int>::const_iterator it = gid2Index_.find(m.first);
                if (it != gid2Index_.end())
                    if (A_[i].find(it->second) != A_[i].end())
                        A_[i][it->second] += m.second;
            }
        }

    private:
751
        const DofMapper& mapper_;
752
753
754
        const std::map<IdType,int>& gid2Index_;
        const std::map<int,IdType>& index2GID_;
        Matrix& A_;
755
756

    }; // class MatEntryExchange
757

758
759
760
    /** \brief communicates values for the sparsity pattern of the new matrix.
        \param A Matrix to operate on.
        \param helper ParallelelISTLHelper.
761
    */
762
    void getExtendedMatrix (Matrix& A, const ParallelISTLHelper<TypeTag>& helper)
763
    {
764
765
        if (gridView_.comm().size() > 1)
        {
766
767
768
            Matrix tmp(A);
            std::size_t nnz=0;
            // get entries from other processes
769
770
771
            MatPatternExchange datahandle(mapper_, gid2Index_, index2GID_, A, helper);
            gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface,
                                              Dune::ForwardCommunication);
772
773
            std::vector<std::set<int> >& sparsity = datahandle.sparsity();
            // add own entries, count number of nonzeros
774
775
776
777
            for (RowIterator i = A.begin(); i != A.end(); ++i)
            {
                for (ColIterator j = A[i.index()].begin(); j != A[i.index()].end(); ++j)
                {
778
                    if (sparsity[i.index()].find(j.index()) == sparsity[i.index()].end())
779
                    {
780
                        sparsity[i.index()].insert(j.index());
781
                    }
782
783
784
                }
                nnz += sparsity[i.index()].size();
            }
785

786
787
788
            A.setSize(tmp.N(), tmp.N(), nnz);
            A.setBuildMode(Matrix::row_wise);
            typename Matrix::CreateIterator citer = A.createbegin();
789
            using Iter = typename std::vector<std::set<int> >::const_iterator;
790
791
            for (Iter i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer)
            {
792
                using SIter = typename std::set<int>::const_iterator;
793
794
795
                for (SIter si = i->begin(), send = i->end(); si!=send; ++si)
                    citer.insert(*si);
            }
796

797
798
799
            // set matrix old values
            A = 0;
            for (RowIterator i = tmp.begin(); i != tmp.end(); ++i)
800
                for (ColIterator j = tmp[i.index()].begin(); j != tmp[i.index()].end(); ++j)
801
802
803
804
                    A[i.index()][j.index()] = tmp[i.index()][j.index()];
        }
    }

805
806
    /** \brief Sums up the entries corresponding to border vertices.
        \param A Matrix to operate on.
807
808
809
    */
    void sumEntries (Matrix& A)
    {
810
        if (gridView_.comm().size() > 1)
811
        {
812
813
814
            MatEntryExchange datahandle(mapper_, gid2Index_, index2GID_, A);
            gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface,
                                              Dune::ForwardCommunication);
815
816
817
        }
    }

818
819
#if HAVE_MPI
    /**
820
821
     * \brief Extends the sparsity pattern of the discretization matrix for AMG.
     * \param A A reference to the matrix to change.
822
823
824
825
     */
    void getExtendedMatrix (Matrix& A) const;
#endif

826
private:
827
828
829
830
831
832
    const GridView gridView_;
    const DofMapper& mapper_;
    std::map<IdType, int> gid2Index_;
    std::map<int, IdType> index2GID_;

}; // class EntityExchanger
833

834
// methods that only exist if MPI is available
835
#if HAVE_MPI
836
837
template<class TypeTag>
void EntityExchanger<TypeTag>::getExtendedMatrix (Matrix& A) const
838
{
839
840
    if (gridView_.comm().size() > 1)
    {
841
        Matrix tmp(A);
842
        std::size_t nnz = 0;
843
        // get entries from other processes
844
845
846
        MatPatternExchange datahandle(mapper_, gid2Index_, index2GID_, A, *this);
        gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface,
                              Dune::ForwardCommunication);
847
        std::vector<std::set<int> >& sparsity = datahandle.sparsity();
848

849
        // add own entries, count number of nonzeros
850
851
852
853
        for (RowIterator i = A.begin(); i != A.end(); ++i)
        {
            for (ColIterator j = A[i.index()].begin(); j != A[i.index()].end(); ++j)
            {
854
855
856
857
858
                if (sparsity[i.index()].find(j.index()) == sparsity[i.index()].end())
                    sparsity[i.index()].insert(j.index());
            }
            nnz += sparsity[i.index()].size();
        }
859

860
861
862
        A.setSize(tmp.N(), tmp.N(), nnz);
        A.setBuildMode(Matrix::row_wise);
        typename Matrix::CreateIterator citer = A.createbegin();
863
        using Iter = typename std::vector<std::set<int> >::const_iterator;
864
        for (Iter i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer){
865
            using SIter = typename std::set<int>::const_iterator;
866
867
868
            for (SIter si = i->begin(), send = i->end(); si!=send; ++si)
                citer.insert(*si);
        }
869

870
871
872
        // set matrix old values
        A = 0;
        for (RowIterator i = tmp.begin(); i != tmp.end(); ++i)
873
            for (ColIterator j = tmp[i.index()].begin(); j != tmp[i.index()].end(); ++j)
874
                A[i.index()][j.index()] = tmp[i.index()][j.index()];
875

876
        sumEntries(A);
877
    }
878
879

} // EntityExchanger::getExtendedMatrix
880

881
template<class TypeTag>
882
883
template<typename MatrixType, typename Comm>
void ParallelISTLHelper<TypeTag>::createIndexSetAndProjectForAMG(MatrixType& m, Comm& comm)
884
{
885
886
    if(!initialized_)
    {
887
888
889
890
891
        // This is the first time this function is called.
        // Therefore we need to initialize the marker vectors for ghosts and
        // owned dofs
        initGhostsAndOwners();
    }
892
893

    // First find out which dofs we share with other processors
894
    std::vector<int> sharedDofs(mapper_.size(), false);
895

896
    SharedGatherScatter sgs(sharedDofs, mapper_);
897

898
899
900
    if (gridView_.comm().size() > 1)
        gridView_.communicate(sgs, Dune::All_All_Interface,
                              Dune::ForwardCommunication);
901
902

    // Count shared dofs that we own
903
904
905
    using GlobalIndex = typename Comm::ParallelIndexSet::GlobalIndex;
    GlobalIndex count = 0;
    auto owned = owner_.begin();
906

907
    for (auto v=sharedDofs.begin(), vend=sharedDofs.end(); v != vend; ++v, ++owned)
908
        if(*v && *owned==1)
909
910
            ++count;

911
912
    Dune::dverb << gridView_.comm().rank() << ": shared count is " << count.touint()
                << std::endl;
913

914
915
    std::vector<GlobalIndex> counts(gridView_.comm().size());
    gridView_.comm().allgather(&count, 1, &(counts[0]));
916
917

    // Compute start index start_p = \sum_{i=0}^{i<p} counts_i
918
919
    GlobalIndex start = 0;
    for (int i = 0; i < gridView_.comm().rank(); ++i)
920
        start = start + counts[i];
921

922
923
    std::vector<GlobalIndex> scalarIndices(mapper_.size(),
                                           std::numeric_limits<GlobalIndex>::max());
924

925
926
    auto shared = sharedDofs.begin();
    auto index = scalarIndices.begin();
927

928
929
930
931
    for (auto i=owner_.begin(), iend=owner_.end(); i!=iend; ++i, ++shared, ++index)
    {
        if(*i==1 && *shared)
        {
932
            *index=start;
933
            ++start;
934
        }
935
    }
936
937

    // publish global indices for the shared DOFS to other processors.
938
939
940
941
    GlobalIndexGatherScatter<GlobalIndex> gigs(scalarIndices, mapper_);
    if (gridView_.comm().size()>1)
        gridView_.communicate(gigs, Dune::All_All_Interface,
                              Dune::ForwardCommunication);
942
943
944


    // Setup the index set
945
946
947
    comm.indexSet().beginResize();
    index = scalarIndices.begin();
    auto ghost = isGhost_.begin();
948

949
    for (auto i=owner_.begin(), iend=owner_.end(); i!=iend; ++i, ++ghost, ++index)
950
951
    {
        Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
952
953
        if (*index!=std::numeric_limits<GlobalIndex>::max())
        {
954
            // global index exist in index set
955
956
            if (*i>0)
            {
957
958
959
                // This dof is managed by us.
                attr = Dune::OwnerOverlapCopyAttributeSet::owner;
            }
960
961
962
            else if ( *ghost==(1<<24) && ( comm.getSolverCategory() ==
                                           static_cast<int>(Dune::SolverCategory::nonoverlapping)) )
            {
963
964
965
                //use attribute overlap for ghosts in novlp grids
                attr = Dune::OwnerOverlapCopyAttributeSet::overlap;
            }
966
967
            else
            {
968
969
                attr = Dune::OwnerOverlapCopyAttributeSet::copy;
            }
970
            comm.indexSet().add(*index, typename Comm::ParallelIndexSet::LocalIndex(i-owner_.begin(), attr));
971
        }
972
    }
973
    comm.indexSet().endResize();
974
975
976

    // Compute neighbours using communication
    std::set<int> neighbours;
977
    NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(),
978
979
                               neighbours);

980
981
982
    if (gridView_.comm().size() > 1)
        gridView_.communicate(ngs, Dune::All_All_Interface,
                              Dune::ForwardCommunication);
983

984
985
986
987
    comm.remoteIndices().setNeighbours(neighbours);
    comm.remoteIndices().template rebuild<false>();

} // ParallelISTLHelper::createIndexSetAndProjectForAMG
988

989
#endif // HAVE_MPI
990
991
992
993
994
995
996
997
998
999
1000

/*!
 * \brief Prepare the linear algebra member variables.
 *
 * At compile time, correct constructor calls have to be chosen,
 * depending on whether the setting is parallel or sequential.
 * Since several template parameters are present, this cannot be solved
 * by a full function template specialization. Instead, class template
 * specialization has to be used.
 * This adapts example 4 from http://www.gotw.ca/publications/mill17.htm.
 *
For faster browsing, not all history is shown. View entire blame