boxproblem.hh 14.8 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
3
4
5
6
7
/*****************************************************************************
 *   Copyright (C) 2009 by Andreas Lauser                                    *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
8
 *   This program is free software: you can redistribute it and/or modify    *
9
 *   it under the terms of the GNU General Public License as published by    *
10
11
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
12
 *                                                                           *
13
14
15
16
17
18
19
 *   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/>.   *
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 *****************************************************************************/
/*!
 * \file
 * \brief Base class for all problems which use the box scheme
 */
#ifndef DUMUX_BOX_PROBLEM_HH
#define DUMUX_BOX_PROBLEM_HH

#include "boxproperties.hh"

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

namespace Dumux
{
/*!
Andreas Lauser's avatar
Andreas Lauser committed
36
 * \ingroup BoxModel
Andreas Lauser's avatar
Andreas Lauser committed
37
38
39
40
41
42
 * \brief Base class for all problems which use the box scheme.
 *
 * \note All quantities are specified assuming a threedimensional
 *       world. Problems discretized using 2D grids are assumed to be
 *       extruded by \f$1 m\f$ and 1D grids are assumed to have a
 *       cross section of \f$1m \times 1m\f$.
43
 */
Andreas Lauser's avatar
Andreas Lauser committed
44
template<class TypeTag>
45
46
47
class BoxProblem
{
private:
Andreas Lauser's avatar
Andreas Lauser committed
48
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Implementation;
49
50
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;

Andreas Lauser's avatar
Andreas Lauser committed
51
    typedef Dumux::VtkMultiWriter<GridView> VtkMultiWriter;
52

Andreas Lauser's avatar
Andreas Lauser committed
53
54
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
55

Andreas Lauser's avatar
Andreas Lauser committed
56
57
58
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
59

Andreas Lauser's avatar
Andreas Lauser committed
60
61
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
62
63
64
65
66
67

    enum {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

Andreas Lauser's avatar
Andreas Lauser committed
68
69
70
71
72
    typedef typename GridView::template Codim<0>::Entity Element;

    typedef typename GridView::Grid::ctype CoordScalar;
    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
73

74
75
76
    // copying a problem is not a good idea
    BoxProblem(const BoxProblem &);

77
public:
78
79
80
    /*!
     * \brief Constructor
     *
Felix Bode's avatar
Felix Bode committed
81
     * \param timeManager The TimeManager which is used by the simulation
82
     * \param gridView The simulation's idea about physical space
83
     */
Andreas Lauser's avatar
Andreas Lauser committed
84
    BoxProblem(TimeManager &timeManager, const GridView &gridView)
85
86
87
        : gridView_(gridView),
          bboxMin_(std::numeric_limits<double>::max()),
          bboxMax_(-std::numeric_limits<double>::max()),
Andreas Lauser's avatar
Andreas Lauser committed
88
89
90
          elementMapper_(gridView),
          vertexMapper_(gridView),
          timeManager_(&timeManager),
91
          newtonMethod_(asImp_()),
Andreas Lauser's avatar
Andreas Lauser committed
92
          resultWriter_(asImp_().name())
93
94
95
96
97
98
99
100
101
102
103
    {
        // calculate the bounding box of the grid view
        VertexIterator vIt = gridView.template begin<dim>();
        const VertexIterator vEndIt = gridView.template end<dim>();
        for (; vIt!=vEndIt; ++vIt) {
            for (int i=0; i<dim; i++) {
                bboxMin_[i] = std::min(bboxMin_[i], vIt->geometry().corner(0)[i]);
                bboxMax_[i] = std::max(bboxMax_[i], vIt->geometry().corner(0)[i]);
            }
        }
        for (int i=0; i<dim; i++) {
Bernd Flemisch's avatar
Bernd Flemisch committed
104
105
106
107
108
            if (gridView.comm().size() > 1)
            {
            	bboxMin_[i] = gridView.comm().min(bboxMin_[i]);
            	bboxMax_[i] = gridView.comm().max(bboxMax_[i]);
            }
109
110
111
112
113
114
115
116
117
118
        }
    }

    ~BoxProblem()
    {
    };

    /*!
     * \brief Called by the Dumux::TimeManager in order to
     *        initialize the problem.
119
120
121
     *
     * If you overload this method don't forget to call
     * ParentType::init()
122
123
124
125
     */
    void init()
    {
        // set the initial condition of the model
Andreas Lauser's avatar
Andreas Lauser committed
126
        model().init(asImp_());
127
128
    }

Andreas Lauser's avatar
Andreas Lauser committed
129
130
131
132
133
134
    /*!
     * \brief If model coupling is used, this updates the parameters
     *        required to calculate the coupling fluxes between the
     *        sub-models.
     *
     * By default it does nothing
135
136
137
     *
     * \param element The DUNE Codim<0> entity for which the coupling
     *                parameters should be computed.
Andreas Lauser's avatar
Andreas Lauser committed
138
139
140
     */
    void updateCouplingParams(const Element &element) const
    {}
141

Andreas Lauser's avatar
Andreas Lauser committed
142
143
144
145
146
    /*!
     * \name Simulation steering
     */
    // \{

147
148
149
    /*!
     * \brief Called by the time manager before the time integration.
     */
150
    void preTimeStep()
151
152
153
154
155
156
157
158
    {}

    /*!
     * \brief Called by Dumux::TimeManager in order to do a time
     *        integration on the model.
     */
    void timeIntegration()
    {
Andreas Lauser's avatar
Andreas Lauser committed
159
160
161
162
163
164
        const int maxFails = 10;
        for (int i = 0; i < maxFails; ++i) {
            if (i > 0 && gridView().comm().rank() == 0)
                std::cout << "Newton solver did not converge. Retrying with time step of "
                          << timeManager().timeStepSize() << "sec\n";

165
            if (model_.update(newtonMethod_, newtonCtl_))
Andreas Lauser's avatar
Andreas Lauser committed
166
                return;
167

Andreas Lauser's avatar
Andreas Lauser committed
168
169
170
171
172
            // update failed
            Scalar dt = timeManager().timeStepSize();
            Scalar nextDt = dt / 2;
            timeManager().setTimeStepSize(nextDt);
        }
173

Andreas Lauser's avatar
Andreas Lauser committed
174
175
176
177
178
        DUNE_THROW(Dune::MathError,
                   "Newton solver didn't converge after "
                   << maxFails
                   << " timestep divisions. dt="
                   << timeManager().timeStepSize());
179
180
    }

181
    /*!
182
     * \brief Returns the newton method object
183
184
185
186
187
188
189
190
191
192
     */
    NewtonMethod &newtonMethod()
    { return newtonMethod_; }

    /*!
     * \copydoc newtonMethod()
     */
    const NewtonMethod &newtonMethod() const
    { return newtonMethod_; }

193
194
195
196
197
198
199
200
201
202
203
204
    /*!
     * \brief Returns the newton contoller object
     */
    NewtonController &newtonController()
    { return newtonCtl_; }

    /*!
     * \copydoc newtonController()
     */
    const NewtonController &newtonController() const
    { return newtonCtl_; }

205
206
207
208
    /*!
     * \brief Called by Dumux::TimeManager whenever a solution for a
     *        timestep has been computed and the simulation time has
     *        been updated.
209
210
     *
     * \param dt The current time step size
211
     */
212
213
    Scalar nextTimeStepSize(Scalar dt)
    { return newtonCtl_.suggestTimeStepSize(dt); };
214
215
216
217
218
219
220
221
222

    /*!
     * \brief Returns true if a restart file should be written to
     *        disk.
     *
     * The default behaviour is to write one restart file every 5 time
     * steps. This file is intented to be overwritten by the
     * implementation.
     */
223
    bool shouldWriteRestartFile() const
224
    {
225
        return timeManager().timeStepIndex() > 0 &&
Andreas Lauser's avatar
Andreas Lauser committed
226
            (timeManager().timeStepIndex() % 10 == 0);
227
228
229
230
231
232
233
234
235
236
    }

    /*!
     * \brief Returns true if the current solution should be written to
     *        disk (i.e. as a VTK file)
     *
     * The default behaviour is to write out every the solution for
     * very time step. This file is intented to be overwritten by the
     * implementation.
     */
237
    bool shouldWriteOutput() const
Andreas Lauser's avatar
Andreas Lauser committed
238
    { return true; }
239
240

    /*!
241
242
     * \brief Called by the time manager after the time integration to
     *        do some post processing on the solution.
243
     */
244
    void postTimeStep()
Andreas Lauser's avatar
Andreas Lauser committed
245
    { }
246

247
    /*!
248
249
250
251
252
     * \brief Called by the time manager after everything which can be
     *        done about the current time step is finished and the
     *        model should be prepared to do the next time integration.
     */
    void advanceTimeLevel()
253
    {
254
        model_.advanceTimeLevel();
255
256
257
    }

    /*!
258
     * \brief Called when the end of an simulation episode is reached.
259
260
     *
     * Typically a new episode should be started in this method.
261
262
263
264
265
266
267
     */
    void episodeEnd()
    {
        std::cerr << "The end of an episode is reached, but the problem "
                  << "does not override the episodeEnd() method. "
                  << "Doing nothing!\n";
    };
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    // \}

    /*!
     * \brief The problem name.
     *
     * This is used as a prefix for files generated by the simulation.
     * It could be either overwritten by the problem files, or simply
     * declared over the setName() function in the application file.
     */
    const char *name() const
    {
        return simname_.c_str();
    }

    /*!
     * \brief Set the problem name.
     *
285
286
287
288
289
     * This static method sets the simulation name, which should be
     * called before the application problem is declared! If not, the
     * default name "sim" will be used.
     *
     * \param newName The problem's name
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
     */
    static void setName(const char *newName)
    {
        simname_ = newName;
    }

    /*!
     * \brief The GridView which used by the problem.
     */
    const GridView &gridView() const
    { return gridView_; }

    /*!
     * \brief The coordinate of the corner of the GridView's bounding
     *        box with the smallest values.
     */
    const GlobalPosition &bboxMin() const
    { return bboxMin_; }

    /*!
     * \brief The coordinate of the corner of the GridView's bounding
     *        box with the largest values.
     */
    const GlobalPosition &bboxMax() const
    { return bboxMax_; }

Andreas Lauser's avatar
Andreas Lauser committed
316
317
318
319
320
321
322
323
324
325
326
327
    /*!
     * \brief Returns the mapper for vertices to indices.
     */
    const VertexMapper &vertexMapper() const
    { return vertexMapper_; }

    /*!
     * \brief Returns the mapper for elements to indices.
     */
    const ElementMapper &elementMapper() const
    { return elementMapper_; }

328
329
330
331
    /*!
     * \brief Returns TimeManager object used by the simulation
     */
    TimeManager &timeManager()
Andreas Lauser's avatar
Andreas Lauser committed
332
    { return *timeManager_; }
333
334
335
336
337

    /*!
     * \copydoc timeManager()
     */
    const TimeManager &timeManager() const
Andreas Lauser's avatar
Andreas Lauser committed
338
    { return *timeManager_; }
339
340
341
342
343

    /*!
     * \brief Returns numerical model used for the problem.
     */
    Model &model()
344
    { return model_; }
345
346
347
348
349

    /*!
     * \copydoc model()
     */
    const Model &model() const
350
    { return model_; }
351
352
353
354
355
356
357
358
    // \}

    /*!
     * \name Restart mechanism
     */
    // \{

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
359
     * \brief This method writes the complete state of the simulation
360
361
362
363
364
365
366
367
368
     *        to the harddisk.
     *
     * The file will start with the prefix returned by the name()
     * method, has the current time of the simulation clock in it's
     * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
     * file.)  See Dumux::Restart for details.
     */
    void serialize()
    {
Andreas Lauser's avatar
Andreas Lauser committed
369
        typedef Dumux::Restart Restarter;
370
        Restarter res;
Andreas Lauser's avatar
Andreas Lauser committed
371
372
        res.serializeBegin(asImp_());
        std::cerr << "Serialize to file '" << res.fileName() << "'\n";
373

Andreas Lauser's avatar
Andreas Lauser committed
374
375
376
377
        timeManager().serialize(res);
        asImp_().serialize(res);
        res.serializeEnd();
    }
378

Andreas Lauser's avatar
Andreas Lauser committed
379
380
381
382
383
384
385
386
    /*!
     * \brief This method writes the complete state of the problem
     *        to the harddisk.
     *
     * The file will start with the prefix returned by the name()
     * method, has the current time of the simulation clock in it's
     * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
     * file.)  See Dumux::Restart for details.
387
388
389
390
     *
     * \tparam Restarter The serializer type
     *
     * \param res The serializer object
Andreas Lauser's avatar
Andreas Lauser committed
391
392
393
394
     */
    template <class Restarter>
    void serialize(Restarter &res)
    {
395
396
        resultWriter_.serialize(res);
        model().serialize(res);
Andreas Lauser's avatar
Andreas Lauser committed
397
    }
398

Andreas Lauser's avatar
Andreas Lauser committed
399
400
401
    /*!
     * \brief Load a previously saved state of the whole simulation
     *        from disk.
402
     *
403
404
     * \param tRestart The simulation time on which the program was
     *                 written to disk.
Andreas Lauser's avatar
Andreas Lauser committed
405
406
407
408
     */
    void restart(Scalar tRestart)
    {
        typedef Dumux::Restart Restarter;
409

Andreas Lauser's avatar
Andreas Lauser committed
410
        Restarter res;
411

Andreas Lauser's avatar
Andreas Lauser committed
412
413
414
415
416
        res.deserializeBegin(asImp_(), tRestart);
        std::cout << "Deserialize from file '" << res.fileName() << "'\n";
        timeManager().deserialize(res);
        asImp_().deserialize(res);
        res.deserializeEnd();
417
418
419
420
421
422
423
    }

    /*!
     * \brief This method restores the complete state of the problem
     *        from disk.
     *
     * It is the inverse of the serialize() method.
424
425
426
427
     *
     * \tparam Restarter The deserializer type
     *
     * \param res The desrializer object
428
     */
Andreas Lauser's avatar
Andreas Lauser committed
429
430
    template <class Restarter>
    void deserialize(Restarter &res)
431
432
433
434
435
436
437
    {
        resultWriter_.deserialize(res);
        model().deserialize(res);
    };

    // \}

Andreas Lauser's avatar
Andreas Lauser committed
438
439
440
441
442
    /*!
     * \brief Write the relavant secondar variables of the current
     *        solution into an VTK output file.
     */
    void writeOutput()
443
444
    {
        // write the current result to disk
445
        if (asImp_().shouldWriteOutput()) {
446
            if (gridView().comm().rank() == 0)
Andreas Lauser's avatar
Andreas Lauser committed
447
                std::cout << "Writing result file for \"" << asImp_().name() << "\"\n";
448
449

            // calculate the time _after_ the time was updated
Andreas Lauser's avatar
Andreas Lauser committed
450
451
452
            Scalar t = timeManager().time() + timeManager().timeStepSize();
            resultWriter_.beginTimestep(t, gridView());
            model().addOutputVtkFields(model().curSol(), resultWriter_);
453
454
455
456
            resultWriter_.endTimestep();
        }
    }

Andreas Lauser's avatar
Andreas Lauser committed
457
458
459
460
461
462
463
464
465
protected:
    //! Returns the implementation of the problem (i.e. static polymorphism)
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }

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

466
467
468
469
private:
    static std::string simname_; // a string for the name of the current simulation,
                                  // which could be set by means of an program argument,
                                  // for example.
Andreas Lauser's avatar
Andreas Lauser committed
470
471
472
473
    const GridView gridView_;

    GlobalPosition bboxMin_;
    GlobalPosition bboxMax_;
474

Andreas Lauser's avatar
Andreas Lauser committed
475
476
    ElementMapper elementMapper_;
    VertexMapper vertexMapper_;
477

Andreas Lauser's avatar
Andreas Lauser committed
478
    TimeManager *timeManager_;
479

480
    Model model_;
481

482
    NewtonMethod newtonMethod_;
483
484
    NewtonController newtonCtl_;

Andreas Lauser's avatar
Andreas Lauser committed
485
    VtkMultiWriter resultWriter_;
486
487
488
};
// definition of the static class member simname_,
// which is necessary because it is of type string.
Andreas Lauser's avatar
Andreas Lauser committed
489
490
template <class TypeTag>
std::string BoxProblem<TypeTag>::simname_="sim"; //initialized with default "sim"
491
492
493
494

}

#endif