boxmodel.hh 31.9 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
Andreas Lauser's avatar
Andreas Lauser committed
3
/*****************************************************************************
4
5
 *   Copyright (C) 2008-2010 by Andreas Lauser                               *
 *   Copyright (C) 2008-2010 by Bernd Flemisch                               *
Andreas Lauser's avatar
Andreas Lauser committed
6
 *   Institute for Modelling Hydraulic and Environmental Systems             *
Andreas Lauser's avatar
Andreas Lauser committed
7
8
9
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
10
 *   This program is free software: you can redistribute it and/or modify    *
Andreas Lauser's avatar
Andreas Lauser committed
11
 *   it under the terms of the GNU General Public License as published by    *
12
13
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
Andreas Lauser's avatar
Andreas Lauser committed
14
 *                                                                           *
15
16
17
18
19
20
21
 *   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/>.   *
Andreas Lauser's avatar
Andreas Lauser committed
22
 *****************************************************************************/
Bernd Flemisch's avatar
Bernd Flemisch committed
23
24
25
26
27
/*!
 * \file
 *
 * \brief Base class for models using box discretization
 */
28
29
#ifndef DUMUX_BOX_MODEL_HH
#define DUMUX_BOX_MODEL_HH
Andreas Lauser's avatar
Andreas Lauser committed
30
31

#include "boxproperties.hh"
Andreas Lauser's avatar
Andreas Lauser committed
32
#include "boxpropertydefaults.hh"
Andreas Lauser's avatar
Andreas Lauser committed
33

34
#include "boxelementvolumevariables.hh"
Andreas Lauser's avatar
Andreas Lauser committed
35
36
37
#include "boxlocaljacobian.hh"
#include "boxlocalresidual.hh"

38
39
#include <dumux/parallel/vertexhandles.hh>

40
#include <dune/grid/common/geometry.hh>
41

Andreas Lauser's avatar
Andreas Lauser committed
42
43
44
45
46
47
48
49
50
51
52
53
namespace Dumux
{

/*!
 * \ingroup BoxModel
 *
 * \brief The base class for the vertex centered finite volume
 *        discretization scheme.
 */
template<class TypeTag>
class BoxModel
{
54
55
56
57
58
59
60
61
62
63
64
65
    typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
Andreas Lauser's avatar
Andreas Lauser committed
66

Andreas Lauser's avatar
Andreas Lauser committed
67
    enum {
68
        numEq = GET_PROP_VALUE(TypeTag, NumEq),
69
        dim = GridView::dimension
Andreas Lauser's avatar
Andreas Lauser committed
70
71
    };

72
73
74
75
76
    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, LocalJacobian) LocalJacobian;
    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
    typedef typename GET_PROP_TYPE(TypeTag, NewtonController) NewtonController;
Andreas Lauser's avatar
Andreas Lauser committed
77

78
    typedef typename GridView::ctype CoordScalar;
Andreas Lauser's avatar
Andreas Lauser committed
79
80
81
82
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::template Codim<dim>::Entity Vertex;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
83
84
85
86
    typedef typename GridView::IntersectionIterator IntersectionIterator;

    typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
    typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
Andreas Lauser's avatar
Andreas Lauser committed
87

88
89
90
    // copying a model is not a good idea
    BoxModel(const BoxModel &);

Andreas Lauser's avatar
Andreas Lauser committed
91
92
93
94
95
public:
    /*!
     * \brief The constructor.
     */
    BoxModel()
96
97
98
    {
        enableHints_ = GET_PARAM(TypeTag, bool, EnableHints);
    }
Andreas Lauser's avatar
Andreas Lauser committed
99
100
101
102
103
104

    ~BoxModel()
    { delete jacAsm_;  }

    /*!
     * \brief Apply the initial conditions to the model.
105
     *
106
     * \param problem The object representing the problem which needs to
107
     *             be simulated.
Andreas Lauser's avatar
Andreas Lauser committed
108
     */
109
    void init(Problem &problem)
Andreas Lauser's avatar
Andreas Lauser committed
110
    {
111
        problemPtr_ = &problem;
112

113
114
        updateBoundaryIndices_();

115
        int nDofs = asImp_().numDofs();
116
117
        uCur_.resize(nDofs);
        uPrev_.resize(nDofs);
Andreas Lauser's avatar
Andreas Lauser committed
118
119
120
121
122
        boxVolume_.resize(nDofs);

        localJacobian_.init(problem_());
        jacAsm_ = new JacobianAssembler();
        jacAsm_->init(problem_());
123

124
        asImp_().applyInitialSolution_();
Andreas Lauser's avatar
Andreas Lauser committed
125

Andreas Lauser's avatar
Andreas Lauser committed
126
        // resize the hint vectors
127
        if (enableHints_) {
Andreas Lauser's avatar
Andreas Lauser committed
128
            int nVerts = gridView_().size(dim);
129
130
            curHints_.resize(nVerts);
            prevHints_.resize(nVerts);
Andreas Lauser's avatar
Andreas Lauser committed
131
            hintsUsable_.resize(nVerts);
Andreas Lauser's avatar
Andreas Lauser committed
132
            std::fill(hintsUsable_.begin(),
Andreas Lauser's avatar
Andreas Lauser committed
133
134
135
136
                      hintsUsable_.end(),
                      false);
        }

Andreas Lauser's avatar
Andreas Lauser committed
137
138
        // also set the solution of the "previous" time step to the
        // initial solution.
139
        uPrev_ = uCur_;
Andreas Lauser's avatar
Andreas Lauser committed
140
141
    }

142
    void setHints(const Element &element,
143
144
                  ElementVolumeVariables &prevVolVars,
                  ElementVolumeVariables &curVolVars) const
Andreas Lauser's avatar
Andreas Lauser committed
145
    {
146
        if (!enableHints_)
Andreas Lauser's avatar
Andreas Lauser committed
147
148
            return;

149
        int n = element.template count<dim>();
150
151
152
        prevVolVars.resize(n);
        curVolVars.resize(n);
        for (int i = 0; i < n; ++i) {
153
            int globalIdx = vertexMapper().map(element, i, dim);
154
155
156
157
158
159
160
161
162
163
164

            if (!hintsUsable_[globalIdx]) {
                curVolVars[i].setHint(NULL);
                prevVolVars[i].setHint(NULL);
            }
            else {
                curVolVars[i].setHint(&curHints_[globalIdx]);
                prevVolVars[i].setHint(&prevHints_[globalIdx]);
            }
        }
    };
Andreas Lauser's avatar
Andreas Lauser committed
165

166
    void setHints(const Element &element,
167
                  ElementVolumeVariables &curVolVars) const
Andreas Lauser's avatar
Andreas Lauser committed
168
    {
169
        if (!enableHints_)
Andreas Lauser's avatar
Andreas Lauser committed
170
171
            return;

172
        int n = element.template count<dim>();
173
174
        curVolVars.resize(n);
        for (int i = 0; i < n; ++i) {
175
            int globalIdx = vertexMapper().map(element, i, dim);
176

177
178
179
180
181
            if (!hintsUsable_[globalIdx])
                curVolVars[i].setHint(NULL);
            else
                curVolVars[i].setHint(&curHints_[globalIdx]);
        }
Andreas Lauser's avatar
Andreas Lauser committed
182
183
    };

184
    void updatePrevHints()
Andreas Lauser's avatar
Andreas Lauser committed
185
    {
186
        if (!enableHints_)
Andreas Lauser's avatar
Andreas Lauser committed
187
188
            return;

189
        prevHints_ = curHints_;
Andreas Lauser's avatar
Andreas Lauser committed
190
191
    };

192
193
    void updateCurHints(const Element &element,
                        const ElementVolumeVariables &elemVolVars) const
Andreas Lauser's avatar
Andreas Lauser committed
194
    {
195
        if (!enableHints_)
Andreas Lauser's avatar
Andreas Lauser committed
196
            return;
197

198
        for (unsigned int i = 0; i < elemVolVars.size(); ++i) {
199
200
            int globalIdx = vertexMapper().map(element, i, dim);
            curHints_[globalIdx] = elemVolVars[i];
201
            if (!hintsUsable_[globalIdx])
202
                prevHints_[globalIdx] = elemVolVars[i];
Andreas Lauser's avatar
Andreas Lauser committed
203
204
205
206
207
            hintsUsable_[globalIdx] = true;
        }
    };


Andreas Lauser's avatar
Andreas Lauser committed
208
209
210
    /*!
     * \brief Compute the global residual for an arbitrary solution
     *        vector.
211
     *
212
     * \param residual Stores the result
213
     * \param u The solution for which the residual ought to be calculated
Andreas Lauser's avatar
Andreas Lauser committed
214
     */
215
    Scalar globalResidual(SolutionVector &residual,
216
                          const SolutionVector &u)
Andreas Lauser's avatar
Andreas Lauser committed
217
    {
218
219
        SolutionVector tmp(curSol());
        curSol() = u;
220
        Scalar res = globalResidual(residual);
221
        curSol() = tmp;
Andreas Lauser's avatar
Andreas Lauser committed
222
223
224
225
226
227
        return res;
    }

    /*!
     * \brief Compute the global residual for the current solution
     *        vector.
228
     *
229
     * \param residual Stores the result
Andreas Lauser's avatar
Andreas Lauser committed
230
     */
231
    Scalar globalResidual(SolutionVector &residual)
Andreas Lauser's avatar
Andreas Lauser committed
232
    {
233
        residual = 0;
234

Andreas Lauser's avatar
Andreas Lauser committed
235
236
237
        ElementIterator elemIt = gridView_().template begin<0>();
        const ElementIterator elemEndIt = gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
238
            localResidual().eval(*elemIt);
Andreas Lauser's avatar
Andreas Lauser committed
239

240
241
            for (int i = 0; i < elemIt->template count<dim>(); ++i) {
                int globalI = vertexMapper().map(*elemIt, i, dim);
242
                residual[globalI] += localResidual().residual(i);
Andreas Lauser's avatar
Andreas Lauser committed
243
244
245
            }
        };

246
        // calculate the square norm of the residual
247
        Scalar result2 = residual.two_norm2();
248
249
        if (gridView_().comm().size() > 1)
            result2 = gridView_().comm().sum(result2);
250
251

        // add up the residuals on the process borders
252
        if (gridView_().comm().size() > 1) {
Andreas Lauser's avatar
Andreas Lauser committed
253
            VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
254
                sumHandle(residual, vertexMapper());
255
            gridView_().communicate(sumHandle,
Andreas Lauser's avatar
Andreas Lauser committed
256
257
                                    Dune::InteriorBorder_InteriorBorder_Interface,
                                    Dune::ForwardCommunication);
258
259
260
        }

        return std::sqrt(result2);
Andreas Lauser's avatar
Andreas Lauser committed
261
262
    }

Andreas Lauser's avatar
Andreas Lauser committed
263
    /*!
264
265
266
     * \brief Compute the integral over the domain of the storage
     *        terms of all conservation quantities.
     *
267
     * \param storage Stores the result
Andreas Lauser's avatar
Andreas Lauser committed
268
     */
269
    void globalStorage(PrimaryVariables &storage)
Andreas Lauser's avatar
Andreas Lauser committed
270
    {
271
        storage = 0;
272

Andreas Lauser's avatar
Andreas Lauser committed
273
274
275
        ElementIterator elemIt = gridView_().template begin<0>();
        const ElementIterator elemEndIt = gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
276
            localResidual().evalStorage(*elemIt);
Andreas Lauser's avatar
Andreas Lauser committed
277
278

            for (int i = 0; i < elemIt->template count<dim>(); ++i)
279
                storage += localResidual().storageTerm()[i];
Andreas Lauser's avatar
Andreas Lauser committed
280
281
        };

282
        if (gridView_().comm().size() > 1)
283
            storage = gridView_().comm().sum(storage);
Andreas Lauser's avatar
Andreas Lauser committed
284
285
    }

Andreas Lauser's avatar
Andreas Lauser committed
286
    /*!
Felix Bode's avatar
Felix Bode committed
287
     * \brief Returns the volume \f$\mathrm{[m^3]}\f$ of a given control volume.
288
289
290
     *
     * \param globalIdx The global index of the control volume's
     *                  associated vertex
Andreas Lauser's avatar
Andreas Lauser committed
291
292
293
294
295
     */
    Scalar boxVolume(int globalIdx) const
    { return boxVolume_[globalIdx][0]; }

    /*!
296
297
298
299
300
301
302
303
304
305
306
307
308
     * \brief Reference to the current solution as a block vector.
     */
    const SolutionVector &curSol() const
    { return uCur_; }

    /*!
     * \brief Reference to the current solution as a block vector.
     */
    SolutionVector &curSol()
    { return uCur_; }

    /*!
     * \brief Reference to the previous solution as a block vector.
Andreas Lauser's avatar
Andreas Lauser committed
309
     */
310
311
    const SolutionVector &prevSol() const
    { return uPrev_; }
Andreas Lauser's avatar
Andreas Lauser committed
312
313

    /*!
314
     * \brief Reference to the previous solution as a block vector.
Andreas Lauser's avatar
Andreas Lauser committed
315
     */
316
317
    SolutionVector &prevSol()
    { return uPrev_; }
Andreas Lauser's avatar
Andreas Lauser committed
318
319
320
321
322
323
324
325

    /*!
     * \brief Returns the operator assembler for the global jacobian of
     *        the problem.
     */
    JacobianAssembler &jacobianAssembler()
    { return *jacAsm_; }

326
327
328
329
330
331
    /*!
     * \copydoc jacobianAssembler()
     */
    const JacobianAssembler &jacobianAssembler() const
    { return *jacAsm_; }

Andreas Lauser's avatar
Andreas Lauser committed
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
    /*!
     * \brief Returns the local jacobian which calculates the local
     *        stiffness matrix for an arbitrary element.
     *
     * The local stiffness matrices of the element are used by
     * the jacobian assembler to produce a global linerization of the
     * problem.
     */
    LocalJacobian &localJacobian()
    { return localJacobian_; }
    /*!
     * \copydoc localJacobian()
     */
    const LocalJacobian &localJacobian() const
    { return localJacobian_; }

    /*!
     * \brief Returns the local residual function.
     */
    LocalResidual &localResidual()
    { return localJacobian().localResidual(); }
    /*!
     * \copydoc localResidual()
     */
    const LocalResidual &localResidual() const
    { return localJacobian().localResidual(); }

Andreas Lauser's avatar
Andreas Lauser committed
359
360
361
    /*!
     * \brief Returns the relative weight of a primary variable for
     *        calculating relative errors.
362
363
364
     *
     * \param vertIdx The global index of the control volume
     * \param pvIdx The index of the primary variable
Andreas Lauser's avatar
Andreas Lauser committed
365
366
     */
    Scalar primaryVarWeight(int vertIdx, int pvIdx) const
367
    {
368
        return 1.0/std::max(std::abs(this->prevSol()[vertIdx][pvIdx]), 1.0);
369
    }
Andreas Lauser's avatar
Andreas Lauser committed
370
371
372
373
374

    /*!
     * \brief Returns the relative error between two vectors of
     *        primary variables.
     *
375
376
377
378
379
     * \param vertexIdx The global index of the control volume's
     *                  associated vertex
     * \param pv1 The first vector of primary variables
     * \param pv2 The second vector of primary variables
     *
Andreas Lauser's avatar
Andreas Lauser committed
380
381
382
383
384
385
386
387
388
389
390
391
     * \todo The vertexIdx argument is pretty hacky. it is required by
     *       models with pseudo primary variables (i.e. the primary
     *       variable switching models). the clean solution would be
     *       to access the pseudo primary variables from the primary
     *       variables.
     */
    Scalar relativeErrorVertex(int vertexIdx,
                               const PrimaryVariables &pv1,
                               const PrimaryVariables &pv2)
    {
        Scalar result = 0.0;
        for (int j = 0; j < numEq; ++j) {
392
393
394
            //Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
            //Scalar eqErr = std::abs(pv1[j] - pv2[j])*weight;
            Scalar eqErr = std::abs(pv1[j] - pv2[j]);
395
            eqErr /= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
Andreas Lauser's avatar
Andreas Lauser committed
396
397
398
399
400

            result = std::max(result, eqErr);
        }
        return result;
    }
401

Andreas Lauser's avatar
Andreas Lauser committed
402
403
    /*!
     * \brief Try to progress the model to the next timestep.
404
405
     *
     * \param solver The non-linear solver
406
     * \param controller The controller which specifies the behaviour
407
     *                   of the non-linear solver
Andreas Lauser's avatar
Andreas Lauser committed
408
409
410
411
412
     */
    bool update(NewtonMethod &solver,
                NewtonController &controller)
    {
#if HAVE_VALGRIND
413
414
        for (size_t i = 0; i < curSol().size(); ++i)
            Valgrind::CheckDefined(curSol()[i]);
Andreas Lauser's avatar
Andreas Lauser committed
415
416
417
418
419
#endif // HAVE_VALGRIND

        asImp_().updateBegin();

        bool converged = solver.execute(controller);
420
        if (converged) {
Andreas Lauser's avatar
Andreas Lauser committed
421
            asImp_().updateSuccessful();
422
        }
Andreas Lauser's avatar
Andreas Lauser committed
423
424
        else
            asImp_().updateFailed();
425

Andreas Lauser's avatar
Andreas Lauser committed
426
#if HAVE_VALGRIND
427
428
        for (size_t i = 0; i < curSol().size(); ++i) {
            Valgrind::CheckDefined(curSol()[i]);
Andreas Lauser's avatar
Andreas Lauser committed
429
430
        }
#endif // HAVE_VALGRIND
431

Andreas Lauser's avatar
Andreas Lauser committed
432
433
434
435
436
437
438
439
440
441
        return converged;
    }


    /*!
     * \brief Called by the update() method before it tries to
     *        apply the newton method. This is primary a hook
     *        which the actual model can overload.
     */
    void updateBegin()
442
    { }
Andreas Lauser's avatar
Andreas Lauser committed
443
444
445
446
447
448
449
450


    /*!
     * \brief Called by the update() method if it was
     *        successful. This is primary a hook which the actual
     *        model can overload.
     */
    void updateSuccessful()
451
    { };
Andreas Lauser's avatar
Andreas Lauser committed
452
453
454
455
456
457
458
459
460
461
462

    /*!
     * \brief Called by the update() method if it was
     *        unsuccessful. This is primary a hook which the actual
     *        model can overload.
     */
    void updateFailed()
    {
        // Reset the current solution to the one of the
        // previous time step so that we can start the next
        // update at a physically meaningful solution.
463
464
        uCur_ = uPrev_;
        curHints_ = prevHints_;
Andreas Lauser's avatar
Andreas Lauser committed
465

466
        jacAsm_->reassembleAll();
Andreas Lauser's avatar
Andreas Lauser committed
467
468
    };

469
    /*!
470
471
472
     * \brief Called by the problem if a time integration was
     *        successful, post processing of the solution is done and
     *        the result has been written to disk.
473
     *
474
     * This should prepare the model for the next time integration.
475
476
477
478
     */
    void advanceTimeLevel()
    {
        // make the current solution the previous one.
479
480
        uPrev_ = uCur_;
        prevHints_ = curHints_;
Andreas Lauser's avatar
Andreas Lauser committed
481

482
        updatePrevHints();
483
484
    }

Andreas Lauser's avatar
Andreas Lauser committed
485
486
    /*!
     * \brief Serializes the current state of the model.
487
488
489
490
     *
     * \tparam Restarter The type of the serializer class
     *
     * \param res The serializer object
Andreas Lauser's avatar
Andreas Lauser committed
491
492
493
494
495
496
497
     */
    template <class Restarter>
    void serialize(Restarter &res)
    { res.template serializeEntities<dim>(asImp_(), this->gridView_()); }

    /*!
     * \brief Deserializes the state of the model.
498
499
500
501
     *
     * \tparam Restarter The type of the serializer class
     *
     * \param res The serializer object
Andreas Lauser's avatar
Andreas Lauser committed
502
503
504
505
506
     */
    template <class Restarter>
    void deserialize(Restarter &res)
    {
        res.template deserializeEntities<dim>(asImp_(), this->gridView_());
507
        prevSol() = curSol();
Andreas Lauser's avatar
Andreas Lauser committed
508
509
510
511
512
    }

    /*!
     * \brief Write the current solution for a vertex to a restart
     *        file.
513
514
515
     *
     * \param outstream The stream into which the vertex data should
     *                  be serialized to
516
     * \param vertex The DUNE Codim<dim> entity which's data should be
517
     *             serialized
Andreas Lauser's avatar
Andreas Lauser committed
518
519
     */
    void serializeEntity(std::ostream &outstream,
520
                         const Vertex &vertex)
Andreas Lauser's avatar
Andreas Lauser committed
521
    {
522
        int vertIdx = dofMapper().map(vertex);
Andreas Lauser's avatar
Andreas Lauser committed
523
524
525
526
527
528
529
530
531

        // write phase state
        if (!outstream.good()) {
            DUNE_THROW(Dune::IOError,
                       "Could not serialize vertex "
                       << vertIdx);
        }

        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
532
            outstream << curSol()[vertIdx][eqIdx] << " ";
Andreas Lauser's avatar
Andreas Lauser committed
533
534
535
536
537
538
        }
    };

    /*!
     * \brief Reads the current solution variables for a vertex from a
     *        restart file.
539
540
541
     *
     * \param instream The stream from which the vertex data should
     *                  be deserialized from
542
     * \param vertex The DUNE Codim<dim> entity which's data should be
543
     *             deserialized
Andreas Lauser's avatar
Andreas Lauser committed
544
545
     */
    void deserializeEntity(std::istream &instream,
546
                           const Vertex &vertex)
Andreas Lauser's avatar
Andreas Lauser committed
547
    {
548
        int vertIdx = dofMapper().map(vertex);
Andreas Lauser's avatar
Andreas Lauser committed
549
550
551
552
553
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            if (!instream.good())
                DUNE_THROW(Dune::IOError,
                           "Could not deserialize vertex "
                           << vertIdx);
554
            instream >> curSol()[vertIdx][eqIdx];
Andreas Lauser's avatar
Andreas Lauser committed
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
        }
    };

    /*!
     * \brief Returns the number of global degrees of freedoms (DOFs)
     */
    size_t numDofs() const
    { return gridView_().size(dim); }

    /*!
     * \brief Mapper for the entities where degrees of freedoms are
     *        defined to indices.
     *
     * This usually means a mapper for vertices.
     */
    const DofMapper &dofMapper() const
    { return problem_().vertexMapper(); };

    /*!
     * \brief Mapper for vertices to indices.
     */
    const VertexMapper &vertexMapper() const
    { return problem_().vertexMapper(); };

    /*!
     * \brief Mapper for elements to indices.
     */
    const ElementMapper &elementMapper() const
    { return problem_().elementMapper(); };

585
    /*!
586
     * \brief Resets the Jacobian matrix assembler, so that the
587
588
     *        boundary types can be altered.
     */
Andreas Lauser's avatar
Andreas Lauser committed
589
590
591
592
593
594
595
    void resetJacobianAssembler ()
    {
        delete jacAsm_;
        jacAsm_ = new JacobianAssembler;
        jacAsm_->init(problem_());
    }

596
597
598
599
600
601
602
    /*!
     * \brief Update the weights of all primary variables within an
     *        element given the complete set of volume variables
     *
     * \param element The DUNE codim 0 entity
     * \param volVars All volume variables for the element
     */
603
604
    void updatePVWeights(const Element &element,
                         const ElementVolumeVariables &volVars) const
605
606
    { };

Andreas Lauser's avatar
Andreas Lauser committed
607
608
609
610
    /*!
     * \brief Add the vector fields for analysing the convergence of
     *        the newton method to the a VTK multi writer.
     *
611
612
613
     * \tparam MultiWriter The type of the VTK multi writer
     *
     * \param writer  The VTK multi writer object on which the fields should be added.
Bernd Flemisch's avatar
Bernd Flemisch committed
614
615
     * \param u       The solution function
     * \param deltaU  The delta of the solution function before and after the Newton update
Andreas Lauser's avatar
Andreas Lauser committed
616
617
618
     */
    template <class MultiWriter>
    void addConvergenceVtkFields(MultiWriter &writer,
619
                                 const SolutionVector &u,
Andreas Lauser's avatar
Andreas Lauser committed
620
621
                                 const SolutionVector &deltaU)
    {
622
        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
Andreas Lauser's avatar
Andreas Lauser committed
623

624
625
        SolutionVector residual(u);
        asImp_().globalResidual(residual, u);
Andreas Lauser's avatar
Andreas Lauser committed
626
627
628
629
630
631
632
633

        // create the required scalar fields
        unsigned numVertices = this->gridView_().size(dim);

        // global defect of the two auxiliary equations
        ScalarField* def[numEq];
        ScalarField* delta[numEq];
        ScalarField* x[numEq];
634
635
636
637
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            x[eqIdx] = writer.allocateManagedBuffer(numVertices);
            delta[eqIdx] = writer.allocateManagedBuffer(numVertices);
            def[eqIdx] = writer.allocateManagedBuffer(numVertices);
Andreas Lauser's avatar
Andreas Lauser committed
638
639
640
641
642
643
644
        }

        VertexIterator vIt = this->gridView_().template begin<dim>();
        VertexIterator vEndIt = this->gridView_().template end<dim>();
        for (; vIt != vEndIt; ++ vIt)
        {
            int globalIdx = vertexMapper().map(*vIt);
645
646
647
648
            for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                (*x[eqIdx])[globalIdx] = u[globalIdx][eqIdx];
                (*delta[eqIdx])[globalIdx] = - deltaU[globalIdx][eqIdx];
                (*def[eqIdx])[globalIdx] = residual[globalIdx][eqIdx];
Andreas Lauser's avatar
Andreas Lauser committed
649
650
651
            }
        }

652
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
653
            std::ostringstream oss;
654
655
656
657
658
659
            oss.str(""); oss << "x_" << eqIdx;
            writer.attachVertexData(*x[eqIdx], oss.str());
            oss.str(""); oss << "delta_" << eqIdx;
            writer.attachVertexData(*delta[eqIdx], oss.str());
            oss.str(""); oss << "defect_" << eqIdx;
            writer.attachVertexData(*def[eqIdx], oss.str());
Andreas Lauser's avatar
Andreas Lauser committed
660
661
662
663
664
665
        }

        asImp_().addOutputVtkFields(u, writer);
    }

    /*!
666
     * \brief Add the quantities of a time step which ought to be written to disk.
Andreas Lauser's avatar
Andreas Lauser committed
667
     *
668
     * This should be overwritten by the actual model if any secondary
Andreas Lauser's avatar
Andreas Lauser committed
669
670
671
     * variables should be written out. Read: This should _always_ be
     * overwritten by well behaved models!
     *
672
673
     * \tparam MultiWriter The type of the VTK multi writer
     *
674
675
     * \param sol The global vector of primary variable values.
     * \param writer The VTK multi writer where the fields should be added.
Andreas Lauser's avatar
Andreas Lauser committed
676
677
     */
    template <class MultiWriter>
678
    void addOutputVtkFields(const SolutionVector &sol,
Andreas Lauser's avatar
Andreas Lauser committed
679
680
                            MultiWriter &writer)
    {
681
        typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
Andreas Lauser's avatar
Andreas Lauser committed
682
683
684

        // create the required scalar fields
        unsigned numVertices = this->gridView_().size(dim);
685

Andreas Lauser's avatar
Andreas Lauser committed
686
687
        // global defect of the two auxiliary equations
        ScalarField* x[numEq];
688
689
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            x[eqIdx] = writer.allocateManagedBuffer(numVertices);
Andreas Lauser's avatar
Andreas Lauser committed
690
691
692
693
694
695
696
        }

        VertexIterator vIt = this->gridView_().template begin<dim>();
        VertexIterator vEndIt = this->gridView_().template end<dim>();
        for (; vIt != vEndIt; ++ vIt)
        {
            int globalIdx = vertexMapper().map(*vIt);
697
698
            for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                (*x[eqIdx])[globalIdx] = sol[globalIdx][eqIdx];
Andreas Lauser's avatar
Andreas Lauser committed
699
700
701
            }
        }

702
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
703
            std::ostringstream oss;
704
705
            oss << "primaryVar_" << eqIdx;
            writer.attachVertexData(*x[eqIdx], oss.str());
706
        }
Andreas Lauser's avatar
Andreas Lauser committed
707
708
    }

Bernd Flemisch's avatar
Bernd Flemisch committed
709
710
711
712
713
    /*!
     * \brief Reference to the grid view of the spatial domain.
     */
    const GridView &gridView() const
    { return problem_().gridView(); }
Andreas Lauser's avatar
Andreas Lauser committed
714

715
716
717
718
719
720
721
722
723
724
725
726
727
728
    /*!
     * \brief Returns true if the vertex with 'globalVertIdx' is
     *        located on the grid's boundary.
     *
     * \param globalVertIdx The global index of the control volume's
     *                      associated vertex
     */
    bool onBoundary(int globalVertIdx) const
    { return boundaryIndices_[globalVertIdx]; }

    /*!
     * \brief Returns true if a vertex is located on the grid's
     *        boundary.
     *
729
     * \param element A DUNE Codim<0> entity which contains the control
730
     *             volume's associated vertex.
731
     * \param vIdx The local vertex index inside element
732
     */
733
734
    bool onBoundary(const Element &element, int vIdx) const
    { return onBoundary(vertexMapper().map(element, vIdx, dim)); }
735

736
737
738
739
740
741
742
743
744
745
746
747
748
749
    /*!
     * \brief Fill the fluid state according to the primary variables. 
     * 
     * Taking the information from the primary variables, 
     * the fluid state is filled with every information that is 
     * necessary to evaluate the model's local residual. 
     * 
     * \param primaryVariables The primary variables of the model. 
     * \param problem The problem at hand. 
     * \param element The current element. 
     * \param elementGeometry The finite volume element geometry. 
     * \param scvIdx The index of the subcontrol volume. 
     * \param fluidState The fluid state to fill. 
     */
750
751
752
753
754
755
756
757
    template <class FluidState>
    static void completeFluidState(const PrimaryVariables& primaryVariables,
                                   const Problem& problem,
                                   const Element& element,
                                   const FVElementGeometry& elementGeometry,
                                   int scvIdx,
                                   FluidState& fluidState)
    {
Andreas Lauser's avatar
Andreas Lauser committed
758
759
        VolumeVariables::completeFluidState(primaryVariables, problem, element,
                                            elementGeometry, scvIdx, fluidState);
760
    }
Andreas Lauser's avatar
Andreas Lauser committed
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
protected:
    /*!
     * \brief A reference to the problem on which the model is applied.
     */
    Problem &problem_()
    { return *problemPtr_; }
    /*!
     * \copydoc problem_()
     */
    const Problem &problem_() const
    { return *problemPtr_; }

    /*!
     * \brief Reference to the grid view of the spatial domain.
     */
    const GridView &gridView_() const
    { return problem_().gridView(); }

779
780
781
    /*!
     * \brief Reference to the local residal object
     */
Andreas Lauser's avatar
Andreas Lauser committed
782
783
784
    LocalResidual &localResidual_()
    { return localJacobian_.localResidual(); }

785
    /*!
786
     * \brief Applies the initial solution for all vertices of the grid.
787
     */
Andreas Lauser's avatar
Andreas Lauser committed
788
789
790
    void applyInitialSolution_()
    {
        // first set the whole domain to zero
791
        uCur_ = Scalar(0.0);
Andreas Lauser's avatar
Andreas Lauser committed
792
        boxVolume_ = Scalar(0.0);
793

794
        FVElementGeometry fvGeometry;
795

Andreas Lauser's avatar
Andreas Lauser committed
796
797
798
799
800
        // iterate through leaf grid and evaluate initial
        // condition at the center of each sub control volume
        //
        // TODO: the initial condition needs to be unique for
        // each vertex. we should think about the API...
801
802
803
        ElementIterator eIt = gridView_().template begin<0>();
        const ElementIterator &eEndIt = gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt) {
Andreas Lauser's avatar
Andreas Lauser committed
804
            // deal with the current element
805
            fvGeometry.update(gridView_(), *eIt);
Andreas Lauser's avatar
Andreas Lauser committed
806
807

            // loop over all element vertices, i.e. sub control volumes
808
            for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; scvIdx++)
Andreas Lauser's avatar
Andreas Lauser committed
809
810
            {
                // map the local vertex index to the global one
811
                int globalIdx = vertexMapper().map(*eIt,
Andreas Lauser's avatar
Andreas Lauser committed
812
813
814
815
816
                                                   scvIdx,
                                                   dim);

                // let the problem do the dirty work of nailing down
                // the initial solution.
817
                PrimaryVariables initVal;
Andreas Lauser's avatar
Andreas Lauser committed
818
819
                Valgrind::SetUndefined(initVal);
                problem_().initial(initVal,
820
821
                                   *eIt,
                                   fvGeometry,
Andreas Lauser's avatar
Andreas Lauser committed
822
823
824
825
                                   scvIdx);
                Valgrind::CheckDefined(initVal);

                // add up the initial values of all sub-control
826
                // volumes. If the initial values disagree for
Andreas Lauser's avatar
Andreas Lauser committed
827
828
                // different sub control volumes, the initial value
                // will be the arithmetic mean.
829
830
                initVal *= fvGeometry.subContVol[scvIdx].volume;
                boxVolume_[globalIdx] += fvGeometry.subContVol[scvIdx].volume;
831
832
                uCur_[globalIdx] += initVal;
                Valgrind::CheckDefined(uCur_[globalIdx]);
Andreas Lauser's avatar
Andreas Lauser committed
833
834
            }
        }
835
836
837

        // add up the primary variables and the volumes of the boxes
        // which cross process borders
838
        if (gridView_().comm().size() > 1) {
Andreas Lauser's avatar
Andreas Lauser committed
839
            VertexHandleSum<Dune::FieldVector<Scalar, 1>,
840
841
                Dune::BlockVector<Dune::FieldVector<Scalar, 1> >,
                VertexMapper> sumVolumeHandle(boxVolume_, vertexMapper());
842
            gridView_().communicate(sumVolumeHandle,
Andreas Lauser's avatar
Andreas Lauser committed
843
844
                                    Dune::InteriorBorder_InteriorBorder_Interface,
                                    Dune::ForwardCommunication);
845

Andreas Lauser's avatar
Andreas Lauser committed
846
            VertexHandleSum<PrimaryVariables, SolutionVector, VertexMapper>
847
                sumPVHandle(uCur_, vertexMapper());
848
            gridView_().communicate(sumPVHandle,
Andreas Lauser's avatar
Andreas Lauser committed
849
850
                                    Dune::InteriorBorder_InteriorBorder_Interface,
                                    Dune::ForwardCommunication);
851
852
        }

Andreas Lauser's avatar
Andreas Lauser committed
853
854
855
        // divide all primary variables by the volume of their boxes
        int n = gridView_().size(dim);
        for (int i = 0; i < n; ++i) {
856
            uCur_[i] /= boxVolume(i);
Andreas Lauser's avatar
Andreas Lauser committed
857
858
859
        }
    }

860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
    /*!
     * \brief Find all indices of boundary vertices.
     *
     * For this we need to loop over all intersections (which is slow
     * in general). If the DUNE grid interface would provide a
     * onBoundary() method for entities this could be done in a much
     * nicer way (actually this would not be necessary)
     */
    void updateBoundaryIndices_()
    {
        boundaryIndices_.resize(numDofs());
        std::fill(boundaryIndices_.begin(), boundaryIndices_.end(), false);

        ElementIterator eIt = gridView_().template begin<0>();
        ElementIterator eEndIt = gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt) {
            Dune::GeometryType geoType = eIt->geometry().type();
            const ReferenceElement &refElem = ReferenceElements::general(geoType);

            IntersectionIterator isIt = gridView_().ibegin(*eIt);
            IntersectionIterator isEndIt = gridView_().iend(*eIt);
            for (; isIt != isEndIt; ++isIt) {
                if (!isIt->boundary())
                    continue;
                // add all vertices on the intersection to the set of
                // boundary vertices
                int faceIdx = isIt->indexInInside();
                int numFaceVerts = refElem.size(faceIdx, 1, dim);
                for (int faceVertIdx = 0;
                     faceVertIdx < numFaceVerts;
                     ++faceVertIdx)
                {
                    int elemVertIdx = refElem.subEntity(faceIdx,
                                                        1,
                                                        faceVertIdx,
                                                        dim);
                    int globalVertIdx = vertexMapper().map(*eIt, elemVertIdx, dim);
                    boundaryIndices_[globalVertIdx] = true;
                }
            }
        }
    }

Andreas Lauser's avatar
Andreas Lauser committed
903
904
905
    // the hint cache for the previous and the current volume
    // variables
    mutable std::vector<bool> hintsUsable_;
906
907
    mutable std::vector<VolumeVariables> curHints_;
    mutable std::vector<VolumeVariables> prevHints_;
Andreas Lauser's avatar
Andreas Lauser committed
908

Andreas Lauser's avatar
Andreas Lauser committed
909
910
911
912
913
    // the problem we want to solve. defines the constitutive
    // relations, matxerial laws, etc.
    Problem *problemPtr_;

    // calculates the local jacobian matrix for a given element
914
    LocalJacobian localJacobian_;
Andreas Lauser's avatar
Andreas Lauser committed
915
916
917
918
    // Linearizes the problem at the current time step using the
    // local jacobian
    JacobianAssembler *jacAsm_;

919
920
921
    // the set of all indices of vertices on the boundary
    std::vector<bool> boundaryIndices_;

Andreas Lauser's avatar
Andreas Lauser committed
922
923
    // cur is the current iterative solution, prev the converged
    // solution of the previous time step
924
925
    SolutionVector uCur_;
    SolutionVector uPrev_;
926

927
    Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
928
929

private:
930
931
932
933
934
935
936
937
938
939
940
    /*!
     * \brief Returns whether messages should be printed
     */
    bool verbose_() const
    { return gridView_().comm().rank() == 0; };

    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }

941
    bool enableHints_;
Andreas Lauser's avatar
Andreas Lauser committed
942
943
944
945
};
}

#endif