parallelhelpers.hh 30.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// -*- 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 3 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/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup Linear
 * \brief Provides a helper class for nonoverlapping
 *        decomposition.
 */
25
26
27
28
#ifndef DUMUX_LINEAR_PARALLELHELPERS_HH
#define DUMUX_LINEAR_PARALLELHELPERS_HH

#if HAVE_MPI
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

#include <dune/common/version.hh>
#include <dune/geometry/dimension.hh>
#include <dune/grid/common/datahandleif.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>

namespace Dumux {

/*!
 * \ingroup Linear
 * \brief A parallel helper class providing a nonoverlapping
 *        decomposition of all degrees of freedom
 */
// operator that resets result to zero at constrained DOFS
45
template<class LinearSolverTraits>
46
47
class ParallelISTLHelper
{
48
    using GridView = typename LinearSolverTraits::GridView;
49
50
51
    using DofMapper = typename LinearSolverTraits::DofMapper;
    enum { dofCodim = LinearSolverTraits::dofCodim };

52
53
54
    // TODO: this is some large number (replace by limits?)
    static constexpr std::size_t ghostMarker_ = 1<<24;

55
56
57
58
59
60
61
62
63
64
65
66
67
    class BaseGatherScatter
    {
    public:
        BaseGatherScatter(const DofMapper& mapper)
        : mapper_(mapper) {}

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

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

68
69
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
        //! returns true if size per entity of given dim and codim is a constant
70
        bool fixedSize(int dim, int codim) const
71
        { return true; }
72
73
74
75
76
#else
        //! returns true if size per entity of given dim and codim is a constant
        bool fixedsize(int dim, int codim) const
        { return true; }
#endif
77
78

        template<class EntityType>
79
        std::size_t size(EntityType& e) const
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
        { return 1; }

        template<class EntityType>
        bool isNeitherInteriorNorBorderEntity(EntityType& e) const
        { return e.partitionType() != Dune::InteriorEntity && e.partitionType() != Dune::BorderEntity; }

    private:
        const DofMapper& mapper_;
    };

    /*!
     * \brief GatherScatter implementation that makes a right hand side in the box model consistent.
     */
    template<class V>
    class ConsistencyBoxGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<ConsistencyBoxGatherScatter<V>,typename V::block_type>
    {
    public:
        using DataType = typename V::block_type;
        using BaseGatherScatter::contains;
101
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
102
        using BaseGatherScatter::fixedSize;
103
104
105
#else
        using BaseGatherScatter::fixedsize;
#endif
106
107
108
109
110
111
112
113
114
115
116
117
118
        using BaseGatherScatter::size;

        ConsistencyBoxGatherScatter(V& container, const DofMapper& mapper)
        : BaseGatherScatter(mapper), container_(container)
        {}

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, const EntityType& e) const
        {
            buff.write(container_[this->index(e)]);
        }

        template<class MessageBuffer, class EntityType>
119
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
120
121
122
123
124
125
126
127
128
129
130
        {
            typename V::block_type block;
            buff.read(block);
            container_[this->index(e)] += block;
        }
    private:
        V& container_;
    };


    /*!
131
     * \brief Writes ghostMarker_ to each data item (of the container) that is gathered or scattered
132
133
134
135
136
137
138
139
140
141
142
     * and is neither interior nor border.
     *
     * Can be used to mark ghost cells.
     */
    class GhostGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<GhostGatherScatter,std::size_t>
    {
    public:
        using DataType = std::size_t;
        using BaseGatherScatter::contains;
143
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
144
        using BaseGatherScatter::fixedSize;
145
146
147
#else
        using BaseGatherScatter::fixedsize;
#endif
148
149
150
151
152
153
154
155
156
        using BaseGatherScatter::size;

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

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, const EntityType& e) const
        {
157
            auto& data = ranks_[this->index(e)];
158
            if (this->isNeitherInteriorNorBorderEntity(e))
159
                data = ghostMarker_;
160
161
162
163
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
164
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
165
166
        {
            std::size_t x;
167
168
            buff.read(x);
            auto& data = ranks_[this->index(e)];
169
            if (this->isNeitherInteriorNorBorderEntity(e))
170
                data = ghostMarker_;
171
172
173
174
175
176
        }
    private:
        std::vector<std::size_t>& ranks_;
    };

    /*!
177
     * \brief GatherScatter handle that sets ghostMarker_ for data items neither associated to
178
179
180
181
182
183
184
185
186
187
188
     * the interior or border and take the minimum when scattering.
     *
     * Used to compute an owner rank for each unknown.
     */
    class InteriorBorderGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<InteriorBorderGatherScatter,std::size_t>
    {
    public:
        using DataType = std::size_t;
        using BaseGatherScatter::contains;
189
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
190
        using BaseGatherScatter::fixedSize;
191
192
193
#else
        using BaseGatherScatter::fixedsize;
#endif
194
195
196
197
198
199
200
201
202
        using BaseGatherScatter::size;

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

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, const EntityType& e) const
        {
203
            auto& data = ranks_[this->index(e)];
204
            if (this->isNeitherInteriorNorBorderEntity(e))
205
                data = ghostMarker_;
206
207
208
209
210
211
212
213
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
        {
            std::size_t x;
            buff.read(x);
214
215
216
            auto& data = ranks_[this->index(e)];
            using std::min;
            data = this->isNeitherInteriorNorBorderEntity(e) ? x : min(data,x);
217
218
219
220
221
222
223
224
225
226
227
228
229
        }
    private:
        std::vector<std::size_t>& ranks_;
    };

    /*!
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
     *
     */
    struct NeighbourGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<NeighbourGatherScatter,int>
    {
230
        using DataType = int;
231
        using BaseGatherScatter::contains;
232
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
233
        using BaseGatherScatter::fixedSize;
234
235
236
#else
        using BaseGatherScatter::fixedsize;
#endif
237
238
239
240
241
242
243
244
245
246
247
248
249
        using BaseGatherScatter::size;

        NeighbourGatherScatter(const DofMapper& mapper, int rank, std::set<int>& neighbours)
        : BaseGatherScatter(mapper), rank_(rank), neighbours_(neighbours)
        {}

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, const EntityType& e) const
        {
            buff.write(rank_);
        }

        template<class MessageBuffer, class EntityType>
250
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
251
        {
252
            int x;
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
            buff.read(x);
            neighbours_.insert(x);
        }
    private:
        int rank_;
        std::set<int>& neighbours_;
    };


    /*!
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
     *
     */
    struct SharedGatherScatter
        : public BaseGatherScatter,
          public Dune::CommDataHandleIF<SharedGatherScatter,int>
    {
        using DataType = int;
        using BaseGatherScatter::contains;
272
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
273
        using BaseGatherScatter::fixedSize;
274
275
276
#else
        using BaseGatherScatter::fixedsize;
#endif
277
278
279
280
281
282
283
284
285
286
287
288
289
290
        using BaseGatherScatter::size;

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

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, EntityType& e) const
        {
            int data = true;
            buff.write(data);
        }

        template<class MessageBuffer, class EntityType>
291
        void scatter(MessageBuffer& buff, const EntityType &e, std::size_t n)
292
293
294
        {
            int x;
            buff.read(x);
295
            auto& data = shared_[this->index(e)];
296
297
298
299
300
301
302
303
304
305
            data = data || x;
        }
    private:
        std::vector<int>& shared_;
    };

    /*!
     * \brief GatherScatter handle for finding out about neighbouring processor ranks.
     *
     */
306
    template<typename GlobalIndex>
307
308
    struct GlobalIndexGatherScatter
        : public BaseGatherScatter,
309
          public Dune::CommDataHandleIF<GlobalIndexGatherScatter<GlobalIndex>, GlobalIndex>
310
    {
311
        using DataType = GlobalIndex;
312
        using BaseGatherScatter::contains;
313
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
314
        using BaseGatherScatter::fixedSize;
315
316
317
#else
        using BaseGatherScatter::fixedsize;
#endif
318
319
        using BaseGatherScatter::size;

320
321
        GlobalIndexGatherScatter(std::vector<GlobalIndex>& globalIndices, const DofMapper& mapper)
        : BaseGatherScatter(mapper), globalIndices_(globalIndices)
322
323
324
325
326
        {}

        template<class MessageBuffer, class EntityType>
        void gather(MessageBuffer& buff, const EntityType& e) const
        {
327
            buff.write(globalIndices_[this->index(e)]);
328
329
330
        }

        template<class MessageBuffer, class EntityType>
331
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
332
333
334
        {
            DataType x;
            buff.read(x);
335
336
337
            using std::min;
            auto& data = globalIndices_[this->index(e)];
            data = min(data, x);
338
339
        }
    private:
340
        std::vector<GlobalIndex>& globalIndices_;
341
342
343
344
    };

public:

345
346
347
348
    ParallelISTLHelper(const GridView& gridView, const DofMapper& mapper)
        : gridView_(gridView), mapper_(mapper), initialized_(false)
    {}

349
350
351
352
    // \brief Initializes the markers for ghosts and owners with the correct size and values.
    //
    void initGhostsAndOwners()
    {
353
354
        const auto rank = gridView_.comm().rank();
        isOwned_.resize(mapper_.size(), rank);
355
        // find out about ghosts
356
        GhostGatherScatter ggs(isOwned_, mapper_);
357
358
359
360

        if (gridView_.comm().size() > 1)
            gridView_.communicate(ggs, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication);

361
        isGhost_ = isOwned_;
362
363

        // partition interior/border
364
        InteriorBorderGatherScatter dh(isOwned_, mapper_);
365
366
367
368
369

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

        // convert vector into mask vector
370
371
        for (auto& v : isOwned_)
            v = (v == rank) ? 1 : 0;
372
373
374
375

        initialized_ = true;
    }

376
377
378
    bool isGhost(std::size_t i) const
    { return isGhost_[i] == ghostMarker_; }

379
380
381
382
383
384
385
386
387
388
389
    // \brief Make a vector of the box model consistent.
    template<class B, class A>
    void makeNonOverlappingConsistent(Dune::BlockVector<B,A>& v)
    {
        ConsistencyBoxGatherScatter<Dune::BlockVector<B,A> > gs(v, mapper_);
        if (gridView_.comm().size() > 1)
            gridView_.communicate(gs, Dune::InteriorBorder_InteriorBorder_Interface,
                                  Dune::ForwardCommunication);
    }

    /*!
390
     * \brief Creates a parallel index set
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
     *
     * \tparam Comm The type of the OwnerOverlapCopyCommunication
     * communicators.
     */
    template<class Comm>
    void createParallelIndexSet(Comm& comm)
    {
        if (!initialized_)
        {
            // This is the first time this function is called.
            // Therefore we need to initialize the marker vectors for ghosts and
            // owned dofs
            initGhostsAndOwners();
        }

406
407
408
409
410
411
        if (gridView_.comm().size() <= 1)
        {
            comm.remoteIndices().template rebuild<false>();
            return;
        }

412
        // First find out which dofs we share with other processors
413
        std::vector<int> isShared(mapper_.size(), false);
414

415
416
        SharedGatherScatter sgs(isShared, mapper_);
        gridView_.communicate(sgs, Dune::All_All_Interface, Dune::ForwardCommunication);
417
418
419
420

        // Count shared dofs that we own
        using GlobalIndex = typename Comm::ParallelIndexSet::GlobalIndex;
        GlobalIndex count = 0;
421
422
        for (std::size_t i = 0; i < isShared.size(); ++i)
            if (isShared[i] && isOwned_[i] == 1)
423
424
425
                ++count;

        std::vector<GlobalIndex> counts(gridView_.comm().size());
426
        gridView_.comm().allgather(&count, 1, counts.data());
427
428

        // Compute start index start_p = \sum_{i=0}^{i<p} counts_i
429
430
        const int rank = gridView_.comm().rank();
        auto start = std::accumulate(counts.begin(), counts.begin() + rank, GlobalIndex(0));
431

432
        std::vector<GlobalIndex> globalIndices(mapper_.size(), std::numeric_limits<GlobalIndex>::max());
433

434
        for (std::size_t i = 0; i < globalIndices.size(); ++i)
435
        {
436
            if (isOwned_[i] == 1 && isShared[i])
437
            {
438
                globalIndices[i] = start; // GlobalIndex does not implement postfix ++
439
440
441
442
443
                ++start;
            }
        }

        // publish global indices for the shared DOFS to other processors.
444
445
        GlobalIndexGatherScatter<GlobalIndex> gigs(globalIndices, mapper_);
        gridView_.communicate(gigs, Dune::All_All_Interface, Dune::ForwardCommunication);
446

447
        resizeIndexSet_(comm, globalIndices);
448
449
450

        // Compute neighbours using communication
        std::set<int> neighbours;
451
452
        NeighbourGatherScatter ngs(mapper_, gridView_.comm().rank(), neighbours);
        gridView_.communicate(ngs, Dune::All_All_Interface, Dune::ForwardCommunication);
453
454
        comm.remoteIndices().setNeighbours(neighbours);

455
        comm.remoteIndices().template rebuild<false>();
456
457
458
459
460
461
462
463
464
465
    }

    //! Return the dofMapper
    const DofMapper& dofMapper() const
    { return mapper_; }

    //! Return the gridView
    const GridView& gridView() const
    { return gridView_; }

466
private:
467
468
469
470
471
472
473
474
475
476
    template<class Comm, class GlobalIndices>
    void resizeIndexSet_(Comm& comm, const GlobalIndices& globalIndices) const
    {
        comm.indexSet().beginResize();

        for (std::size_t localIndex = 0; localIndex < globalIndices.size(); ++localIndex)
        {
            const auto globalIndex = globalIndices[localIndex];
            if (globalIndex != std::numeric_limits<typename GlobalIndices::value_type>::max())
            {
477
478
479
480
                const bool isOwned = isOwned_[localIndex] > 0;
                const auto attr = getAttribute_(comm, isOwned, isGhost(localIndex));
                using LocalIndex = typename Comm::ParallelIndexSet::LocalIndex;
                comm.indexSet().add(globalIndex, LocalIndex{localIndex, attr});
481
482
483
484
485
486
            }
        }

        comm.indexSet().endResize();
    }

487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
    template<class Comm>
    Dune::OwnerOverlapCopyAttributeSet::AttributeSet
    getAttribute_(const Comm& comm, const bool isOwned, const bool isGhost) const
    {
        if (isOwned)
            return Dune::OwnerOverlapCopyAttributeSet::owner;
#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 7)
        else if (isGhost && (comm.category() == static_cast<int>(Dune::SolverCategory::nonoverlapping)) )
#else
        else if (isGhost && (comm.getSolverCategory() == static_cast<int>(Dune::SolverCategory::nonoverlapping)) )
#endif
            return Dune::OwnerOverlapCopyAttributeSet::overlap;
        else
            return Dune::OwnerOverlapCopyAttributeSet::copy;
    }

503
504
    const GridView gridView_; //!< the grid view
    const DofMapper& mapper_; //!< the dof mapper
505
    std::vector<std::size_t> isOwned_; //!< vector to identify unique decomposition
506
507
508
509
510
511
512
513
514
    std::vector<std::size_t> isGhost_; //!< vector to identify ghost dofs
    bool initialized_; //!< whether isGhost and owner arrays are initialized

}; // class ParallelISTLHelper

/*!
 * \ingroup Linear
 * \brief Helper class for adding up matrix entries on border.
 */
515
516
template<class Matrix, class GridView, class DofMapper, int dofCodim>
class ParallelMatrixHelper
517
{
518
    static constexpr int dim = GridView::dimension;
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
    using Grid = typename GridView::Traits::Grid;
    using BlockType = typename Matrix::block_type;
    using IDS = typename Grid::Traits::GlobalIdSet;
    using IdType = typename IDS::IdType;

    /*!
     * \brief A DataHandle class to exchange matrix sparsity patterns.
     *
     *  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).
542
     *  The MatrixPatternExchange class will find these entries and fills a vector "sparsity",
543
544
     *  that contains all missing connections.
     */
545
    template<class IsGhostFunc>
546
    struct MatrixPatternExchange
547
    : public Dune::CommDataHandleIF<MatrixPatternExchange<IsGhostFunc>, IdType>
548
549
550
551
    {
        //! Export type of data for message buffer
        using DataType = IdType;

552
553
        MatrixPatternExchange(const DofMapper& mapper,
                              const std::map<IdType,int>& globalToLocal,
554
555
556
557
                              const std::map<int,IdType>& localToGlobal,
                              Matrix& A,
                              std::vector<std::set<int>>& sparsity,
                              const IsGhostFunc& isGhost)
558
        : mapper_(mapper), idToIndex_(globalToLocal), indexToID_(localToGlobal)
559
560
561
562
        , sparsity_(sparsity), A_(A), isGhost_(isGhost)
        {
            sparsity_.resize(A.N());
        }
563
564
565
566
567

        /*!
         * \brief Returns true if data for given valid codim should be communicated
         */
        bool contains (int dim, int codim) const
568
        { return (codim == dofCodim); }
569

570
571
572
573
574
575
576
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
        //! returns true if size per entity of given dim and codim is a constant
        bool fixedSize(int dim, int codim) const
        { return false; }
#else
        //! returns true if size per entity of given dim and codim is a constant
        bool fixedsize(int dim, int codim) const
577
        { return false; }
578
#endif
579
580
581
582
583

        /*!
         * \brief How many objects of type DataType have to be sent for a given entity
         */
        template<class EntityType>
584
        std::size_t size(EntityType& e) const
585
        {
586
587
588
589
            const auto rowIdx = mapper_.index(e);
            std::size_t  n = 0;
            for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
                if (indexToID_.count(colIt.index()))
590
591
592
593
594
595
596
597
598
                    n++;

            return n;
        }

        /*!
         * \brief Pack data from user to message buffer
         */
        template<class MessageBuffer, class EntityType>
599
        void gather(MessageBuffer& buff, const EntityType& e) const
600
        {
601
602
            const auto rowIdx = mapper_.index(e);
            for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
603
            {
604
605
                auto it = indexToID_.find(colIt.index());
                if (it != indexToID_.end())
606
607
608
609
610
611
612
613
                    buff.write(it->second);
            }
        }

        /*!
         * \brief Unpack data from message buffer to user
         */
        template<class MessageBuffer, class EntityType>
614
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
615
        {
616
617
            const auto rowIdx = mapper_.index(e);
            for (std::size_t k = 0; k < n; k++)
618
619
620
621
            {
                IdType id;
                buff.read(id);
                // only add entries corresponding to border entities
622
623
624
625
                const auto it = idToIndex_.find(id);
                if (it != idToIndex_.end())
                {
                    const auto colIdx = it->second;
626
                    if (!sparsity_[rowIdx].count(colIdx) && !isGhost_(colIdx))
627
628
                        sparsity_[rowIdx].insert(colIdx);
                }
629
630
631
632
633
            }
        }

    private:
        const DofMapper& mapper_;
634
635
        const std::map<IdType,int>& idToIndex_;
        const std::map<int,IdType>& indexToID_;
636
        std::vector<std::set<int>>& sparsity_;
637
        Matrix& A_;
638
        const IsGhostFunc& isGhost_;
639

640
    }; // class MatrixPatternExchange
641
642

    //! Local matrix blocks associated with the global id set
643
    struct MatrixEntry
644
645
646
    {
        IdType first;
        BlockType second;
647
648
        MatrixEntry (const IdType& f, const BlockType& s) : first(f), second(s) {}
        MatrixEntry () {}
649
650
651
    };

    //! A DataHandle class to exchange matrix entries
652
    struct MatrixEntryExchange
653
    : public Dune::CommDataHandleIF<MatrixEntryExchange, MatrixEntry>
654
655
    {
        //! Export type of data for message buffer
656
        using DataType = MatrixEntry;
657

658
659
660
661
662
        MatrixEntryExchange(const DofMapper& mapper,
                            const std::map<IdType,int>& globalToLocal,
                            const std::map<int,IdType>& localToGlobal,
                            Matrix& A)
        : mapper_(mapper), idToIndex_(globalToLocal), indexToID_(localToGlobal), A_(A)
663
664
665
666
667
        {}

        /*!
         * \brief Returns true if data for given valid codim should be communicated
         */
668
669
        bool contains(int dim, int codim) const
        { return (codim == dofCodim); }
670

671
672
#if DUNE_VERSION_GT(DUNE_GRID,2,7)
        //! returns true if size per entity of given dim and codim is a constant
673
        bool fixedSize(int dim, int codim) const
674
        { return false; }
675
676
677
678
679
#else
        //! returns true if size per entity of given dim and codim is a constant
        bool fixedsize(int dim, int codim) const
        { return false; }
#endif
680
681
682
683
684

        /*!
         * \brief How many objects of type DataType have to be sent for a given entity
         */
        template<class EntityType>
685
        std::size_t size(EntityType& e) const
686
        {
687
688
689
690
            const auto rowIdx = mapper_.index(e);
            std::size_t n = 0;
            for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
                if (indexToID_.count(colIt.index()))
691
692
693
694
695
696
697
698
699
                    n++;

            return n;
        }

        /*!
         * \brief Pack data from user to message buffer
         */
        template<class MessageBuffer, class EntityType>
700
        void gather(MessageBuffer& buff, const EntityType& e) const
701
        {
702
703
            const auto rowIdx = mapper_.index(e);
            for (auto colIt = A_[rowIdx].begin(); colIt != A_[rowIdx].end(); ++colIt)
704
            {
705
706
707
                auto it = indexToID_.find(colIt.index());
                if (it != indexToID_.end())
                    buff.write(MatrixEntry(it->second,*colIt));
708
709
710
711
712
713
714
            }
        }

        /*!
         * \brief Unpack data from message buffer to user
         */
        template<class MessageBuffer, class EntityType>
715
        void scatter(MessageBuffer& buff, const EntityType& e, std::size_t n)
716
        {
717
718
            const auto rowIdx = mapper_.index(e);
            for (std::size_t k = 0; k < n; k++)
719
            {
720
                MatrixEntry m;
721
722
                buff.read(m);
                // only add entries corresponding to border entities
723
724
725
726
                auto it = idToIndex_.find(m.first);
                if (it != idToIndex_.end())
                    if (A_[rowIdx].find(it->second) != A_[rowIdx].end())
                        A_[rowIdx][it->second] += m.second;
727
728
729
730
731
            }
        }

    private:
        const DofMapper& mapper_;
732
733
        const std::map<IdType,int>& idToIndex_;
        const std::map<int,IdType>& indexToID_;
734
735
        Matrix& A_;

736
737
738
739
    }; // class MatrixEntryExchange

public:

740
    ParallelMatrixHelper(const GridView& gridView, const DofMapper& mapper)
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
    : gridView_(gridView), mapper_(mapper)
    {
        idToIndex_.clear();
        indexToID_.clear();

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

                idToIndex_.emplace(dofIdxGlobal, localIdx);
                indexToID_.emplace(localIdx, dofIdxGlobal);
            }
        }
    }
758
759
760
761

    /*!
     * \brief communicates values for the sparsity pattern of the new matrix.
     * \param A Matrix to operate on.
Melanie Lipp's avatar
Melanie Lipp committed
762
     * \param isGhost Function returning if is ghost.
763
     */
764
765
    template<class IsGhostFunc>
    void extendMatrix(Matrix& A, const IsGhostFunc& isGhost)
766
    {
767
768
769
770
771
772
        if (gridView_.comm().size() <= 1)
            return;

        Matrix tmp(A);
        std::size_t nnz = 0;
        // get entries from other processes
773
        std::vector<std::set<int>> sparsity;
774
        MatrixPatternExchange<IsGhostFunc> datahandle(mapper_, idToIndex_, indexToID_, A, sparsity, isGhost);
775
776
777
778
        gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface,
                                          Dune::ForwardCommunication);
        // add own entries, count number of nonzeros
        for (auto rowIt = A.begin(); rowIt != A.end(); ++rowIt)
779
        {
780
781
782
            for (auto colIt = A[rowIt.index()].begin(); colIt != A[rowIt.index()].end(); ++colIt)
                if (!sparsity[rowIt.index()].count(colIt.index()))
                    sparsity[rowIt.index()].insert(colIt.index());
783

784
785
            nnz += sparsity[rowIt.index()].size();
        }
786

787
788
789
790
791
792
793
        A.setSize(tmp.N(), tmp.N(), nnz);
        A.setBuildMode(Matrix::row_wise);
        auto citer = A.createbegin();
        for (auto i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer)
        {
            for (auto si = i->begin(), send = i->end(); si!=send; ++si)
                citer.insert(*si);
794
        }
795
796
797
798
799
800

        // set matrix old values
        A = 0;
        for (auto rowIt = tmp.begin(); rowIt != tmp.end(); ++rowIt)
            for (auto colIt = tmp[rowIt.index()].begin(); colIt != tmp[rowIt.index()].end(); ++colIt)
                A[rowIt.index()][colIt.index()] = tmp[rowIt.index()][colIt.index()];
801
802
803
804
805
806
    }

    /*!
     * \brief Sums up the entries corresponding to border vertices.
     * \param A Matrix to operate on.
     */
807
    void sumEntries(Matrix& A)
808
    {
809
810
811
812
813
814
        if (gridView_.comm().size() <= 1)
            return;

        MatrixEntryExchange datahandle(mapper_, idToIndex_, indexToID_, A);
        gridView_.communicate(datahandle, Dune::InteriorBorder_InteriorBorder_Interface,
                                          Dune::ForwardCommunication);
815
816
    }

817
818
819
820
821
private:
    const GridView gridView_;
    const DofMapper& mapper_;
    std::map<IdType, int> idToIndex_;
    std::map<int, IdType> indexToID_;
822

823
}; // class EntityExchanger
824

825
826
827
/*!
 * \brief Prepare linear algebra variables for parallel solvers
 */
828
829
template<class LinearSolverTraits, class ParallelTraits,
         class Matrix, class Vector, class ParallelHelper>
830
void prepareLinearAlgebraParallel(Matrix& A, Vector& b,
831
832
833
                                  std::shared_ptr<typename ParallelTraits::Comm>& comm,
                                  std::shared_ptr<typename ParallelTraits::LinearOperator>& fop,
                                  std::shared_ptr<typename ParallelTraits::ScalarProduct>& sp,
834
                                  ParallelHelper& pHelper)
835
{
836
    if constexpr (ParallelTraits::isNonOverlapping)
837
838
    {
        // extend the matrix pattern such that it is usable for a parallel solver
839
        // and make right-hand side consistent
840
        using GridView = typename LinearSolverTraits::GridView;
841
842
843
844
845
        using DofMapper = typename LinearSolverTraits::DofMapper;
        static constexpr int dofCodim = LinearSolverTraits::dofCodim;
        ParallelMatrixHelper<Matrix, GridView, DofMapper, dofCodim> matrixHelper(pHelper.gridView(), pHelper.dofMapper());
        matrixHelper.extendMatrix(A, [&pHelper](auto idx){ return pHelper.isGhost(idx); });
        matrixHelper.sumEntries(A);
846

847
        pHelper.makeNonOverlappingConsistent(b);
848
849

        // create commicator, operator, scalar product
850
        using Traits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
851
852
853
854
855
856
857
858
859
        const auto category = Dune::SolverCategory::nonoverlapping;
        comm = std::make_shared<typename Traits::Comm>(pHelper.gridView().comm(), category);
        pHelper.createParallelIndexSet(*comm);
        fop = std::make_shared<typename Traits::LinearOperator>(A, *comm);
        sp = std::make_shared<typename Traits::ScalarProduct>(*comm);
    }
    else
    {
        // create commicator, operator, scalar product
860
        using Traits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
861
862
863
864
865
866
        const auto category = Dune::SolverCategory::overlapping;
        comm = std::make_shared<typename Traits::Comm>(pHelper.gridView().comm(), category);
        pHelper.createParallelIndexSet(*comm);
        fop = std::make_shared<typename Traits::LinearOperator>(A, *comm);
        sp = std::make_shared<typename Traits::ScalarProduct>(*comm);
    }
867
}
868
869
870

} // end namespace Dumux

871
#endif // HAVE_MPI
872
#endif // DUMUX_PARALLELHELPERS_HH