multidomainproblem.hh 14.3 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
// -*- 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/>.   *
 *****************************************************************************/
/*!
 * \file
 * \brief Base class for problems which involve two sub problems
 */
23

24
25
#ifndef DUMUX_MULTIDOMAIN_PROBLEM_HH
#define DUMUX_MULTIDOMAIN_PROBLEM_HH
26

27
28
#include <dune/common/deprecated.hh>

29
30
31
32
33
#include "multidomainmodel.hh"
#include "multidomainnewtoncontroller.hh"
#include "multidomainpropertydefaults.hh"
#include "subdomainpropertydefaults.hh"
#include "multidomainassembler.hh"
34
35

#include <dumux/io/restart.hh>
36

37

38
39
40
41
namespace Dumux
{

/*!
Thomas Fetzer's avatar
Thomas Fetzer committed
42
43
 * \ingroup ImplicitBaseProblems
 * \ingroup MultidomainModel
44
 * \brief Base class for problems which involve two sub problems (multidomain problems)s
45
46
 */
template<class TypeTag>
47
class MultiDomainProblem
48
{
49

50
51
52
53
54
55
56
57
58
59
private:
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;

    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;

60
61
    typedef typename GET_PROP_TYPE(TypeTag, SubDomain1TypeTag) SubDomain1TypeTag;
    typedef typename GET_PROP_TYPE(TypeTag, SubDomain2TypeTag) SubDomain2TypeTag;
62

63
64
    typedef typename GET_PROP_TYPE(SubDomain1TypeTag, LocalResidual) LocalResidual1;
    typedef typename GET_PROP_TYPE(SubDomain2TypeTag, LocalResidual) LocalResidual2;
65

66
67
    typedef typename GET_PROP_TYPE(SubDomain1TypeTag, Problem) SubDomainProblem1;
    typedef typename GET_PROP_TYPE(SubDomain2TypeTag, Problem) SubDomainProblem2;
68

69
70
    typedef typename GET_PROP_TYPE(SubDomain1TypeTag, GridView) SubDomainGridView1;
    typedef typename GET_PROP_TYPE(SubDomain2TypeTag, GridView) SubDomainGridView2;
71

72
     typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MultiDomainGrid;
73

74
75
76
77
    typedef typename MultiDomainGrid::LeafGridView MultiDomainGridView;
    typedef typename MultiDomainGrid::Traits::template Codim<0>::Entity MultiDomainElement;
    typedef typename MultiDomainGrid::SubDomainGrid SubDomainGrid;
    typedef typename SubDomainGrid::template Codim<0>::EntityPointer SubDomainElementPointer;
78

79
    typedef Dune::MultiDomainMCMGMapper<MultiDomainGridView, Dune::MCMGVertexLayout> VertexMapper;
80

81
public:
82
    /*!
83
      * \brief Base class for the multi domain problem
84
      *
85
      * \param mdGrid The multi domain grid
86
87
88
      * \param timeManager The TimeManager which is used by the simulation
      *
      */
89
    MultiDomainProblem(MultiDomainGrid &mdGrid,
90
                       TimeManager &timeManager)
91
        : timeManager_(timeManager)
92
93
94
95
96
97
98
99
100
101
102
        , newtonMethod_(asImp_())
        , newtonCtl_(asImp_())
        , mdGrid_(mdGrid)
        , mdGridView_(mdGrid.leafGridView())
        , mdVertexMapper_(mdGrid_.leafGridView())
        , sdID1_(0)
        , sdID2_(1)
        , sdGrid1_(mdGrid.subDomain(sdID1_))
        , sdGrid2_(mdGrid.subDomain(sdID2_))
        , sdProblem1_(timeManager, sdGrid1_.leafGridView())
        , sdProblem2_(timeManager, sdGrid2_.leafGridView())
103
    {}
104

105
    //! \copydoc Dumux::ImplicitProblem::init()
106
107
108
    void init()
    {
        // initialize the sub-problems
109
110
        sdProblem1().init();
        sdProblem2().init();
111
112
113
114
115
116
117
118

        // set the initial condition of the model
        model().init(asImp_());

        // intialize lagrange multipliers
        asImp_().initMortarElements();
    }

119
    //! \copydoc Dumux::ImplicitProblem::serialize()
120
121
122
    template <class Restarter>
    void serialize(Restarter &res)
    {
Thomas Fetzer's avatar
Thomas Fetzer committed
123
        this->model().serialize(res);
124
125
    }

126
    //! \copydoc Dumux::ImplicitProblem::serialize()
127
128
129
130
131
132
133
134
135
136
137
138
    void serialize()
    {
        typedef Dumux::Restart Restarter;
        Restarter res;
        res.serializeBegin(this->asImp_());
        std::cerr << "Serialize to file '" << res.fileName() << "'\n";

        this->timeManager().serialize(res);
        this->asImp_().serialize(res);
        res.serializeEnd();
    }

139
    //! \copydoc Dumux::ImplicitProblem::restart()
140
141
142
143
144
145
146
147
148
149
150
151
152
    void restart(Scalar tRestart)
    {
        typedef Dumux::Restart Restarter;

        Restarter res;

        res.deserializeBegin(this->asImp_(), tRestart);
        std::cout << "Deserialize from file '" << res.fileName() << "'\n";
        this->timeManager().deserialize(res);
        this->asImp_().deserialize(res);
        res.deserializeEnd();
    }

153
    //! \copydoc Dumux::ImplicitProblem::deserialize()
154
155
156
    template <class Restarter>
    void deserialize(Restarter &res)
    {
157
        this->model().deserialize(res);
158
159
160
161
162
163
164
    }

    /*!
     * \name Simulation control
     */
    // \{

165
166
167
168
169
170
171
    /*!
     * \brief Set the initial time step and the time where the simulation ends
     *        and starts simulation
     *
     * \param dtInitial Initial time step
     * \param tEnd Time, when simulation ends
     */
172
    DUNE_DEPRECATED_MSG("Is deprecated and will be removed without replacement.")
173
174
175
176
177
178
    bool simulate(Scalar dtInitial, Scalar tEnd)
    {
        timeManager_.setEndTime(tEnd);
        timeManager_.setTimeStepSize(dtInitial);
        timeManager_.runSimulation(asImp_());
        return true;
179
    }
180
181
182
183
184
185
186

    /*!
     * \brief Called by the time manager before the time integration. Calls preTimeStep()
     *        of the subproblems.
     */
    void preTimeStep()
    {
187
188
        asImp_().sdProblem1().preTimeStep();
        asImp_().sdProblem2().preTimeStep();
189
190
    }

191
    //! \copydoc Dumux::ImplicitProblem::timeIntegration()
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    void timeIntegration()
    {
        // TODO: should be called from the group Implicit
        const int maxFails =
                GET_PARAM_FROM_GROUP(TypeTag, int, Newton, MaxTimeStepDivisions);
        for (int i = 0; i < maxFails; ++i)
        {
            if (model_.update(newtonMethod_, newtonCtl_))
                return;

            // update failed
            Scalar dt = timeManager().timeStepSize();
            Scalar nextDt = dt / 2;
            timeManager().setTimeStepSize(nextDt);

            std::cout << "Newton solver did not converge. Retrying with time step of "
                      << timeManager().timeStepSize() << "sec\n";
        }

        DUNE_THROW(Dune::MathError,
                   "Newton solver didn't converge after "
                   << maxFails
                   << " timestep divisions. dt="
                   << timeManager().timeStepSize());
    }

    /*!
     * \brief Called by the time manager after the time integration to
     *        do some post processing on the solution. Calls postTimeStep()
     *        of the subproblems.
     */
    void postTimeStep()
    {
225
226
        asImp_().sdProblem1().postTimeStep();
        asImp_().sdProblem2().postTimeStep();
227
228
    }

229
    //! \copydoc Dumux::ImplicitProblem::nextTimeStepSize()
230
231
232
    Scalar nextTimeStepSize(const Scalar dt)
    {
        return newtonCtl_.suggestTimeStepSize(dt);
233
    }
234
235
236
237
238
239
240

    /*!
     * \brief This method is called by the model if the update to the
     *        next time step failed completely.
     */
    void updateSuccessful()
    {
241
        model_.updateSuccessful();
242
    }
243

244
    //! \copydoc Dumux::ImplicitProblem::shouldWriteOutput()
245
246
247
    bool shouldWriteOutput() const
    { return true; }

248
    //! \copydoc Dumux::ImplicitProblem::shouldWriteRestartFile()
249
250
251
    bool shouldWriteRestartFile() const
    { return false; }

252
    //! \copydoc Dumux::ImplicitProblem::episodeEnd()
253
254
255
256
257
258
259
    void episodeEnd()
    {
        std::cerr << "The end of an episode is reached, but the problem "
                  << "does not override the episodeEnd() method. "
                  << "Doing nothing!\n";
    }

260
    //! \copydoc Dumux::ImplicitProblem::advanceTimeLevel()
261
262
    void advanceTimeLevel()
    {
263
264
        asImp_().sdProblem1().advanceTimeLevel();
        asImp_().sdProblem2().advanceTimeLevel();
265
266
267
268

        model_.advanceTimeLevel();
    }

269
    //! \copydoc Dumux::ImplicitProblem::writeOutput()
270
271
272
273
    void writeOutput()
    {
        // write the current result to disk
        if (asImp_().shouldWriteOutput()) {
274
275
            asImp_().sdProblem1().writeOutput();
            asImp_().sdProblem2().writeOutput();
276
277
278
279
        }
    }


280
281
    // \}

282
    //! \copydoc Dumux::ImplicitProblem::name()
283
    const char *name() const
284
    { return simname_.c_str(); }
285

286
    //! \copydoc Dumux::ImplicitProblem::setName()
287
    static void setName(const char *newName)
288
    { simname_ = newName; }
289

290
    //! \copydoc Dumux::ImplicitProblem::timeManager()
291
292
293
    TimeManager &timeManager()
    { return timeManager_; }

294
    //! \copydoc Dumux::ImplicitProblem::timeManager()
295
296
297
    const TimeManager &timeManager() const
    { return timeManager_; }

298
    //! \copydoc Dumux::ImplicitProblem::newtonController()
299
300
301
    NewtonController &newtonController()
    { return newtonCtl_; }

302
    //! \copydoc Dumux::ImplicitProblem::newtonController()
303
304
305
    const NewtonController &newtonController() const
    { return newtonCtl_; }

306
    //! \copydoc Dumux::ImplicitProblem::model()
307
308
309
    Model &model()
    { return model_; }

310
    //! \copydoc Dumux::ImplicitProblem::model()
311
312
313
314
    const Model &model() const
    { return model_; }
    // \}

315
316
317
    /*!
     * \brief Returns the ID of the first domain
     */
318
    const typename MultiDomainGrid::SubDomainIndex sdID1() const
319
    { return sdID1_; }
320

321
322
323
    /*!
     * \brief Returns the ID of the second domain
     */
324
    const typename MultiDomainGrid::SubDomainIndex sdID2() const
325
    { return sdID2_; }
326

327
328
329
    /*!
     * \brief Returns a reference to subproblem1
     */
330
331
    SubDomainProblem1& sdProblem1()
    { return sdProblem1_; }
332

333
334
335
336
337
    /*!
     * \brief Returns a const reference to subproblem1
     */
    const SubDomainProblem1& sdProblem1() const
    { return sdProblem1_; }
338
339

    /*!
340
     * \brief Returns a reference to subproblem2
341
     */
342
343
    SubDomainProblem2& sdProblem2()
    { return sdProblem2_; }
344

345
346
347
348
349
350
    /*!
     * \brief Returns a const reference to subproblem2
     */
    const SubDomainProblem2& sdProblem2() const
    { return sdProblem2_; }

351
    /*!
352
     * \brief Returns a reference to the localresidual1
353
     */
354
    LocalResidual1& localResidual1()
355
    { return sdProblem1().model().localResidual(); }
356
357

    /*!
358
     * \brief Returns a reference to the localresidual2
359
     */
360
    LocalResidual2& localResidual2()
361
    { return sdProblem2().model().localResidual(); }
362
363
364
365

    /*!
     * \brief Returns a reference to the multidomain grid
     */
366
    MultiDomainGrid& mdGrid()
367
    { return mdGrid_; }
368
369
370
371

    /*!
     * \brief Returns a const reference to the multidomain grid
     */
372
    const MultiDomainGrid& mdGrid() const
373
    { return mdGrid_; }
374
375
376
377

    /*!
     * \brief Returns the multidomain gridview
     */
378
    const MultiDomainGridView& mdGridView() const
379
    { return mdGridView_; }
380
381
382
383

    /*!
     * \brief Returns the multidomain gridview
     */
384
    const MultiDomainGridView& gridView() const
385
    { return mdGridView_; }
386
387

    /*!
388
     * \brief Provides a vertex mapper for the multidomain
389
     */
390
391
    VertexMapper& mdVertexMapper()
    { return mdVertexMapper_; }
392
393
394
395

    /*!
     * \brief Returns a const reference to the subdomain1 grid
     */
396
    const SubDomainGrid& sdGrid1() const
397
    { return sdGrid1_; }
398
399
400
401

    /*!
     * \brief Returns a const reference to the subdomain2 grid
     */
402
    const SubDomainGrid& sdGrid2() const
403
    { return sdGrid2_; }
404

405
406
407
408
    /*!
     * \brief Returns the gridview of subdomain1
     */
    const SubDomainGridView1 sdGridView1() const
409
    { return sdGrid1_.leafGridView(); }
410
411
412
413
414

    /*!
     * \brief Returns the gridview of subdomain2
     */
    const SubDomainGridView2 sdGridView2() const
415
    { return sdGrid2_.leafGridView(); }
416

417
418
    /*!
     * \brief Returns a pointer to the subdomain1 element
419
420
     *
     * \param mdElement1 The multi domain element1
421
     */
422
    SubDomainElementPointer sdElementPointer1(const MultiDomainElement& mdElement1)
423
    { return sdGrid1_.subDomainEntityPointer(mdElement1); }
424
425
426

    /*!
     * \brief Returns a pointer to the subdomain2 element
427
     *
428
     * \param mdElement2 The multi domain element2
429
     */
430
    SubDomainElementPointer sdElementPointer2(const MultiDomainElement& mdElement2)
431
    { return sdGrid2_.subDomainEntityPointer(mdElement2); }
432
433


434
protected:
435
    void initMortarElements()
436
    {}
437

438

439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }

    //! \copydoc asImp_()
    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }

private:
    // a string for the name of the current simulation, which could be
    // set by means of an program argument, for example.
    static std::string simname_;

    TimeManager &timeManager_;
    NewtonMethod newtonMethod_;
    NewtonController newtonCtl_;
454

455
    Model model_;
456

457
    MultiDomainGrid &mdGrid_;
458
    const MultiDomainGridView mdGridView_;
459
460
    VertexMapper mdVertexMapper_;

461
462
    typename MultiDomainGrid::SubDomainIndex sdID1_;
    typename MultiDomainGrid::SubDomainIndex sdID2_;
463

464
465
    const SubDomainGrid& sdGrid1_;
    const SubDomainGrid& sdGrid2_;
466

467
468
    SubDomainProblem1 sdProblem1_;
    SubDomainProblem2 sdProblem2_;
469
};
470

471
472
473
// definition of the static class member simname_,
// which is necessary because it is of type string.
template <class TypeTag>
474
std::string MultiDomainProblem<TypeTag>::simname_ = "simCoupled";
475

476
} // namespace Dumux
477

478
#endif // DUMUX_MULTIDOMAIN_PROBLEM_HH