boxlocaljacobian.hh 18.4 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
2
3
4
5
6
7
8
9
10
11
12
13
// $Id$
/*****************************************************************************
 *   Copyright (C) 2007 by Peter Bastian                                     *
 *   Institute of Parallel and Distributed System                            *
 *   Department Simulation of Large Systems                                  *
 *   University of Stuttgart, Germany                                        *
 *                                                                           *
 *   Copyright (C) 2008-2009 by Andreas Lauser                               *
 *   Copyright (C) 2007-2009 by Bernd Flemisch                               *
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
14
 *   This program is free software: you can redistribute it and/or modify    *
Andreas Lauser's avatar
Andreas Lauser committed
15
 *   it under the terms of the GNU General Public License as published by    *
16
17
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
Andreas Lauser's avatar
Andreas Lauser committed
18
 *                                                                           *
19
20
21
22
23
24
25
 *   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
26
27
28
 *****************************************************************************/
/*!
 * \file
29
 * \brief Caculates the Jacobian of the local residual for box models
Andreas Lauser's avatar
Andreas Lauser committed
30
31
32
33
34
35
 */
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH

#include <dune/istl/matrix.hh>

Andreas Lauser's avatar
Andreas Lauser committed
36
#include "boxelementboundarytypes.hh"
Andreas Lauser's avatar
Andreas Lauser committed
37
38
39
40
41

namespace Dumux
{
/*!
 * \ingroup BoxModel
42
 * \brief Caculates the Jacobian of the local residual for box models
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
 *
 * The default behaviour is to use numeric differentiation, i.e.
 * forward or backward differences (2nd order), or central
 * differences (3rd order). The method used is determined by the
 * "NumericDifferenceMethod" property:
 *
 * - if the value of this property is smaller than 0, backward
 *   differences are used, i.e.:
 *   \f[
 \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
 *   \f]
 *
 * - if the value of this property is 0, central
 *   differences are used, i.e.:
 *   \f[
 \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
 *   \f]
 *
 * - if the value of this property is larger than 0, forward
 *   differences are used, i.e.:
 *   \f[
 \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
 *   \f]
 *
 * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
 * is the value of a sub-control volume's primary variable at the
 * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
Andreas Lauser's avatar
Andreas Lauser committed
70
71
72
73
74
75
76
77
78
79
80
 */
template<class TypeTag>
class BoxLocalJacobian
{
private:
    typedef BoxLocalJacobian<TypeTag> ThisType;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
81
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
Andreas Lauser's avatar
Andreas Lauser committed
82
83

    enum {
84
        numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
Andreas Lauser's avatar
Andreas Lauser committed
85

86
        dim = GridView::dimension,
87
88
89
90
        dimWorld = GridView::dimensionworld,

        Red = JacobianAssembler::Red,
        Yellow = JacobianAssembler::Yellow,
91
92
93
94
        Green = JacobianAssembler::Green,

        numDiffMethod = GET_PROP_VALUE(TypeTag,
                                       PTAG(NumericDifferenceMethod))
Andreas Lauser's avatar
Andreas Lauser committed
95
96
    };

97

Andreas Lauser's avatar
Andreas Lauser committed
98
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
99
    typedef typename GridView::Grid::ctype CoordScalar;
Andreas Lauser's avatar
Andreas Lauser committed
100
101
102
103

    typedef Dune::FieldVector<Scalar, dim>       LocalPosition;
    typedef Dune::FieldVector<Scalar, dimWorld>  GlobalPosition;

104
105
106
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename Element::EntityPointer ElementPointer;
Andreas Lauser's avatar
Andreas Lauser committed
107

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

111
112
113
    typedef typename GridView::IntersectionIterator IntersectionIterator;
    typedef typename Element::Geometry Geometry;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
Andreas Lauser's avatar
Andreas Lauser committed
114
115
116

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
117
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
Andreas Lauser's avatar
Andreas Lauser committed
118

119
120
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
Andreas Lauser's avatar
Andreas Lauser committed
121
122
123
124
125
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;

    typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
    typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;

126
127
128
    // copying a local jacobian is not a good idea
    BoxLocalJacobian(const BoxLocalJacobian &);

Andreas Lauser's avatar
Andreas Lauser committed
129
130
public:
    BoxLocalJacobian()
131
    { Valgrind::SetUndefined(problemPtr_); }
Andreas Lauser's avatar
Andreas Lauser committed
132

133
134
135
136
137
138
139
140
141
    
    /*!
     * \brief Initialize the local Jacobian object.
     *
     * At this point we can assume that everything has been allocated,
     * although some objects may not yet be completely initialized.
     *
     * \param prob The problem which we want to simulate.
     */
Andreas Lauser's avatar
Andreas Lauser committed
142
143
144
145
146
147
148
    void init(Problem &prob)
    {
        problemPtr_ = &prob;
        localResidual_.init(prob);
        // assume quadrilinears as elements with most vertices
        A_.setSize(2<<dim, 2<<dim);
    }
149

Andreas Lauser's avatar
Andreas Lauser committed
150
    /*!
151
152
153
154
     * \brief Assemble an element's local Jacobian matrix of the
     *        defect.
     *
     * \param element The DUNE Codim<0> entity which we look at.
Andreas Lauser's avatar
Andreas Lauser committed
155
156
157
158
159
160
     */
    void assemble(const Element &element)
    {
        // set the current grid element and update the element's
        // finite volume geometry
        elemPtr_ = &element;
161
        fvElemGeom_.update(gridView_(), element);
162
        reset_();
163

Andreas Lauser's avatar
Andreas Lauser committed
164
165
166
167
168
169
170
171
172
173
174
175
176
        bcTypes_.update(problem_(), elem_(), fvElemGeom_);

        // this is pretty much a HACK because the internal state of
        // the problem is not supposed to be changed during the
        // evaluation of the residual. (Reasons: It is a violation of
        // abstraction, makes everything more prone to errors and is
        // not thread save.) The real solution are context objects!
        problem_().updateCouplingParams(elem_());

        int numVertices = fvElemGeom_.numVertices;

        // update the secondary variables for the element at the last
        // and the current time levels
177
178
179
180
/*        prevVolVars_.resize(numVertices);
        for (int i = 0; i < numVertices; ++i)
            prevVolVars_[i].setNoHint();
*/
181
        prevVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
182
183
184
                            elem_(),
                            fvElemGeom_,
                            true /* isOldSol? */);
185
186
187
188
189

/*        curVolVars_.resize(numVertices);
        for (int i = 0; i < numVertices; ++i)
            curVolVars_[i].setHint(prevVolVars_[i]);
*/
190
        curVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
191
192
193
                           elem_(),
                           fvElemGeom_,
                           false /* isOldSol? */);
194
195
196
197
198
199
200
201
        // calculate the local residual
        localResidual().eval(elem_(),
                             fvElemGeom_,
                             prevVolVars_,
                             curVolVars_,
                             bcTypes_);
        residual_ = localResidual().residual();
        
Andreas Lauser's avatar
Andreas Lauser committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
        // calculate the local jacobian matrix
        ElementSolutionVector partialDeriv(numVertices);
        for (int j = 0; j < numVertices; j++) {
            for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
                asImp_().evalPartialDerivative_(partialDeriv,
                                                j,
                                                pvIdx);

                // update the local stiffness matrix with the current partial
                // derivatives
                updateLocalJacobian_(j,
                                     pvIdx,
                                     partialDeriv);
            }
        }
    }

    /*!
220
221
     * \brief Returns a reference to the object which calculates the
     *        local residual.
Andreas Lauser's avatar
Andreas Lauser committed
222
223
224
225
226
     */
    const LocalResidual &localResidual() const
    { return localResidual_; }

    /*!
227
228
     * \brief Returns a reference to the object which calculates the
     *        local residual.
Andreas Lauser's avatar
Andreas Lauser committed
229
     */
230
    LocalResidual &localResidual()
Andreas Lauser's avatar
Andreas Lauser committed
231
232
233
    { return localResidual_; }

    /*!
234
     * \brief Returns the Jacobian of the equations at vertex i to the
Andreas Lauser's avatar
Andreas Lauser committed
235
     *        primary variables at vertex j.
236
237
238
239
240
     *
     * \param i The local vertex (or sub-contol volume) index on which
     *          the equations are defined
     * \param j The local vertex (or sub-contol volume) index which holds
     *          primary variables
Andreas Lauser's avatar
Andreas Lauser committed
241
242
243
244
     */
    const MatrixBlock &mat(int i, int j) const
    { return A_[i][j]; }

245
246
247
248
249
250
251
252
253
    /*!
     * \brief Returns the residual of the equations at vertex i.
     *
     * \param i The local vertex (or sub-contol volume) index on which
     *          the equations are defined
     */
    const PrimaryVariables &residual(int i) const
    { return residual_[i]; }

Andreas Lauser's avatar
Andreas Lauser committed
254
255
256
257
258
259
260
261
262
263
protected:
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }

    /*!
     * \brief Returns a reference to the problem.
     */
    const Problem &problem_() const
264
    {
Andreas Lauser's avatar
Andreas Lauser committed
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
        Valgrind::CheckDefined(problemPtr_);
        return *problemPtr_;
    };

    /*!
     * \brief Returns a reference to the grid view.
     */
    const GridView &gridView_() const
    { return problem_().gridView(); }

    /*!
     * \brief Returns a reference to the element.
     */
    const Element &elem_() const
    {
        Valgrind::CheckDefined(elemPtr_);
        return *elemPtr_;
    };

    /*!
     * \brief Returns a reference to the model.
     */
    const Model &model_() const
    { return problem_().model(); };

290
291
292
293
294
295
    /*!
     * \brief Returns a reference to the jacobian assembler.
     */
    const JacobianAssembler &jacAsm_() const
    { return model_().jacobianAssembler(); }

Andreas Lauser's avatar
Andreas Lauser committed
296
297
298
299
300
301
    /*!
     * \brief Returns a reference to the vertex mapper.
     */
    const VertexMapper &vertexMapper_() const
    { return problem_().vertexMapper(); };

302
303
304
305
306
307
308
309
310
311
312
313
314
    /*!
     * \brief Reset the local jacobian matrix to 0
     */
    void reset_()
    {
        int n = elem_().template count<dim>();
        for (int i = 0; i < n; ++ i) {
            for (int j = 0; j < n; ++ j) {
                A_[i][j] = 0.0;
            }
        }
    }

Andreas Lauser's avatar
Andreas Lauser committed
315
316
317
318
319
320
    /*!
     * \brief Compute the partial derivatives to a primary variable at
     *        an degree of freedom.
     *
     * This method can be overwritten by the implementation if a
     * better scheme than central differences ought to be used.
321
322
323
324
325
326
327
328
329
330
331
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
     *
     * The default implementation of this method uses numeric
     * differentiation, i.e. forward or backward differences (2nd
     * order), or central differences (3rd order). The method used is
     * determined by the "NumericDifferenceMethod" property:
     *
     * - if the value of this property is smaller than 0, backward
     *   differences are used, i.e.:
     *   \f[
         \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
     *   \f]
     *
     * - if the value of this property is 0, central
     *   differences are used, i.e.:
     *   \f[
           \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
     *   \f]
     *
     * - if the value of this property is larger than 0, forward
     *   differences are used, i.e.:
     *   \f[
           \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
     *   \f]
     *
     * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$
     * is the value of a sub-control volume's primary variable at the
     * evaluation point and \f$\epsilon\f$ is a small value larger than 0.
     *
     * \param dest The vector storing the partial derivatives of all
     *              equations
     * \param scvIdx The sub-control volume index of the current
     *               finite element for which the partial derivative
     *               ought to be calculated
     * \param pvIdx The index of the primary variable at the scvIdx'
     *              sub-control volume of the current finite element
     *              for which the partial derivative ought to be
     *              calculated
Andreas Lauser's avatar
Andreas Lauser committed
358
359
360
361
362
363
     */
    void evalPartialDerivative_(ElementSolutionVector &dest,
                                int scvIdx,
                                int pvIdx)
    {
        int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
364
365
        PrimaryVariables priVars(model_().curSol()[globalIdx]);
        VolumeVariables origVolVars(curVolVars_[scvIdx]);
Andreas Lauser's avatar
Andreas Lauser committed
366

367
        curVolVars_[scvIdx].setEvalPoint(&origVolVars);
Andreas Lauser's avatar
Andreas Lauser committed
368
        Scalar eps = asImp_().numericEpsilon_(scvIdx, pvIdx);
369
370
        Scalar delta = 0;

371
372
373
374
        if (numDiffMethod >= 0) { 
            // we are not using backward differences, i.e. we need to
            // calculate f(x - \epsilon)

375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
            // deflect primary variables
            priVars[pvIdx] += eps;
            delta += eps;

            // calculate the residual
            curVolVars_[scvIdx].update(priVars,
                                       problem_(),
                                       elem_(),
                                       fvElemGeom_,
                                       scvIdx,
                                       false);
            localResidual().eval(elem_(),
                                 fvElemGeom_,
                                 prevVolVars_,
                                 curVolVars_,
                                 bcTypes_);
            
            // store the residual
            dest = localResidual().residual();
        }
        else {
396
397
398
            // we are using backward differences, i.e. we don't need
            // to calculate f(x - \epsilon) and we can recycle the
            // (already calculated) residual f(x)
399
            dest = residual_;
400
401
402
        }


403
404
405
406
        if (numDiffMethod <= 0) { 
            // we are not using forward differences, i.e. we don't
            // need to calculate f(x + \epsilon)

407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
            // deflect the primary variables
            priVars[pvIdx] -= delta + eps;
            delta += eps;

            // calculate residual again
            curVolVars_[scvIdx].update(priVars,
                                       problem_(),
                                       elem_(),
                                       fvElemGeom_,
                                       scvIdx,
                                       false);
            localResidual().eval(elem_(),
                                 fvElemGeom_,
                                 prevVolVars_,
                                 curVolVars_,
                                 bcTypes_);
            dest -= localResidual().residual();
        }
        else {
426
427
428
            // we are using forward differences, i.e. we don't need to
            // calculate f(x + \epsilon) and we can recycle the
            // (already calculated) residual f(x)
429
            dest -= residual_;
430
        }
431

432
433
        // divide difference in residuals by the magnitude of the
        // deflections between the two function evaluation
434
        dest /= delta;
Andreas Lauser's avatar
Andreas Lauser committed
435

436
        // restore the orignal state of the element's volume variables
437
        curVolVars_[scvIdx] = origVolVars;
Andreas Lauser's avatar
Andreas Lauser committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454

#if HAVE_VALGRIND
        for (unsigned i = 0; i < dest.size(); ++i)
            Valgrind::CheckDefined(dest[i]);
#endif
    }

    /*!
     * \brief Returns the epsilon value which is added and removed
     *        from the current solution.
     *
     * \param scvIdx     The local index of the element's vertex for
     *                   which the local derivative ought to be calculated.
     * \param pvIdx      The index of the primary variable which gets varied
     */
    Scalar numericEpsilon_(int scvIdx,
                           int pvIdx) const
455
456
    {
        Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
Andreas Lauser's avatar
Andreas Lauser committed
457
458
459
460
        return 1e-9*(std::abs(pv) + 1);
    }

    /*!
461
     * \brief Updates the current local Jacobian matrix with the
Andreas Lauser's avatar
Andreas Lauser committed
462
     *        partial derivatives of all equations in regard to the
463
     *        primary variable 'pvIdx' at vertex 'scvIdx' .
Andreas Lauser's avatar
Andreas Lauser committed
464
465
466
     */
    void updateLocalJacobian_(int scvIdx,
                              int pvIdx,
467
                              const ElementSolutionVector &deriv)
Andreas Lauser's avatar
Andreas Lauser committed
468
    {
469
470
471
472
473
474
475
        for (int i = 0; i < fvElemGeom_.numVertices; i++)
        {
            if (jacAsm_().vertexColor(elem_(), i) == Green) {
                // Green vertices are not to be changed!
                continue;
            }

Andreas Lauser's avatar
Andreas Lauser committed
476
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
477
478
479
480
481
                // A[i][scvIdx][eqIdx][pvIdx] is the rate of change of
                // the residual of equation 'eqIdx' at vertex 'i'
                // depending on the primary variable 'pvIdx' at vertex
                // 'scvIdx'.
                this->A_[i][scvIdx][eqIdx][pvIdx] = deriv[i][eqIdx];
Andreas Lauser's avatar
Andreas Lauser committed
482
483
484
485
486
487
488
                Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
            }
        }
    }

    const Element *elemPtr_;
    FVElementGeometry fvElemGeom_;
489

Andreas Lauser's avatar
Andreas Lauser committed
490
491
492
493
    ElementBoundaryTypes bcTypes_;

    // The problem we would like to solve
    Problem *problemPtr_;
494

Andreas Lauser's avatar
Andreas Lauser committed
495
496
    // secondary variables at the previous and at the current time
    // levels
497
498
499
    ElementVolumeVariables prevVolVars_;
    ElementVolumeVariables curVolVars_;

Andreas Lauser's avatar
Andreas Lauser committed
500
    LocalResidual localResidual_;
501

Andreas Lauser's avatar
Andreas Lauser committed
502
    LocalBlockMatrix A_;
503
504
    ElementSolutionVector residual_;
};
Andreas Lauser's avatar
Andreas Lauser committed
505
506
507
}

#endif