boxlocaljacobian.hh 14.7 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// $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                          *
 *                                                                           *
 *   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, as long as this copyright notice    *
 *   is included in its original form.                                       *
 *                                                                           *
 *   This program is distributed WITHOUT ANY WARRANTY.                       *
 *****************************************************************************/
/*!
 * \file
24
 * \brief Caculates the Jacobian of the local residual for box models
Andreas Lauser's avatar
Andreas Lauser committed
25
26
27
28
29
30
 */
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH

#include <dune/istl/matrix.hh>

Andreas Lauser's avatar
Andreas Lauser committed
31
#include "boxelementboundarytypes.hh"
Andreas Lauser's avatar
Andreas Lauser committed
32
33
34
35
36

namespace Dumux
{
/*!
 * \ingroup BoxModel
37
 * \brief Caculates the Jacobian of the local residual for box models
Andreas Lauser's avatar
Andreas Lauser committed
38
39
40
41
42
43
44
45
46
47
48
 */
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;
49
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
Andreas Lauser's avatar
Andreas Lauser committed
50
51

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

54
        dim = GridView::dimension,
55
56
57
58
        dimWorld = GridView::dimensionworld,

        Red = JacobianAssembler::Red,
        Yellow = JacobianAssembler::Yellow,
59
60
61
62
        Green = JacobianAssembler::Green,

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

65

Andreas Lauser's avatar
Andreas Lauser committed
66
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
67
    typedef typename GridView::Grid::ctype CoordScalar;
Andreas Lauser's avatar
Andreas Lauser committed
68
69
70
71

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

72
73
74
    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
75

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

79
80
81
    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
82
83
84

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

87
88
    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
89
90
91
92
93
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;

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

94
95
96
    // copying a local jacobian is not a good idea
    BoxLocalJacobian(const BoxLocalJacobian &);

Andreas Lauser's avatar
Andreas Lauser committed
97
98
public:
    BoxLocalJacobian()
99
    { Valgrind::SetUndefined(problemPtr_); }
Andreas Lauser's avatar
Andreas Lauser committed
100

101
102
103
104
105
106
107
108
109
    
    /*!
     * \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
110
111
112
113
114
115
116
    void init(Problem &prob)
    {
        problemPtr_ = &prob;
        localResidual_.init(prob);
        // assume quadrilinears as elements with most vertices
        A_.setSize(2<<dim, 2<<dim);
    }
117

Andreas Lauser's avatar
Andreas Lauser committed
118
    /*!
119
120
121
122
     * \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
123
124
125
126
127
128
     */
    void assemble(const Element &element)
    {
        // set the current grid element and update the element's
        // finite volume geometry
        elemPtr_ = &element;
129
        fvElemGeom_.update(gridView_(), element);
130
        reset_();
131

Andreas Lauser's avatar
Andreas Lauser committed
132
133
134
135
136
137
138
139
140
141
142
143
144
        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
145
146
147
148
/*        prevVolVars_.resize(numVertices);
        for (int i = 0; i < numVertices; ++i)
            prevVolVars_[i].setNoHint();
*/
149
        prevVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
150
151
152
                            elem_(),
                            fvElemGeom_,
                            true /* isOldSol? */);
153
154
155
156
157

/*        curVolVars_.resize(numVertices);
        for (int i = 0; i < numVertices; ++i)
            curVolVars_[i].setHint(prevVolVars_[i]);
*/
158
        curVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
159
160
161
                           elem_(),
                           fvElemGeom_,
                           false /* isOldSol? */);
162
163
164
165
166
167
168
169
        // calculate the local residual
        localResidual().eval(elem_(),
                             fvElemGeom_,
                             prevVolVars_,
                             curVolVars_,
                             bcTypes_);
        residual_ = localResidual().residual();
        
Andreas Lauser's avatar
Andreas Lauser committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
        // 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);
            }
        }
    }

    /*!
188
189
     * \brief Returns a reference to the object which calculates the
     *        local residual.
Andreas Lauser's avatar
Andreas Lauser committed
190
191
192
193
194
     */
    const LocalResidual &localResidual() const
    { return localResidual_; }

    /*!
195
196
     * \brief Returns a reference to the object which calculates the
     *        local residual.
Andreas Lauser's avatar
Andreas Lauser committed
197
     */
198
    LocalResidual &localResidual()
Andreas Lauser's avatar
Andreas Lauser committed
199
200
201
    { return localResidual_; }

    /*!
202
     * \brief Returns the Jacobian of the equations at vertex i to the
Andreas Lauser's avatar
Andreas Lauser committed
203
     *        primary variables at vertex j.
204
205
206
207
208
     *
     * \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
209
210
211
212
     */
    const MatrixBlock &mat(int i, int j) const
    { return A_[i][j]; }

213
214
215
216
217
218
219
220
221
    /*!
     * \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
222
223
224
225
226
227
228
229
230
231
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
232
    {
Andreas Lauser's avatar
Andreas Lauser committed
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
        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(); };

258
259
260
261
262
263
    /*!
     * \brief Returns a reference to the jacobian assembler.
     */
    const JacobianAssembler &jacAsm_() const
    { return model_().jacobianAssembler(); }

Andreas Lauser's avatar
Andreas Lauser committed
264
265
266
267
268
269
    /*!
     * \brief Returns a reference to the vertex mapper.
     */
    const VertexMapper &vertexMapper_() const
    { return problem_().vertexMapper(); };

270
271
272
273
274
275
276
277
278
279
280
281
282
    /*!
     * \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
283
284
285
286
287
288
289
290
291
292
293
294
    /*!
     * \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.
     */
    void evalPartialDerivative_(ElementSolutionVector &dest,
                                int scvIdx,
                                int pvIdx)
    {
        int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
295
296
        PrimaryVariables priVars(model_().curSol()[globalIdx]);
        VolumeVariables origVolVars(curVolVars_[scvIdx]);
Andreas Lauser's avatar
Andreas Lauser committed
297

298
        curVolVars_[scvIdx].setEvalPoint(&origVolVars);
Andreas Lauser's avatar
Andreas Lauser committed
299
        Scalar eps = asImp_().numericEpsilon_(scvIdx, pvIdx);
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
        Scalar delta = 0;

        if (numDiffMethod >= 0) { // not backward differences
            // 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 {
            // backward differences
325
            dest = residual_;
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
        }


        if (numDiffMethod <= 0) { // not forward differences
            // 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 {
            // forward differences
350
            dest -= residual_;
351
        }
352

353
        dest /= delta;
Andreas Lauser's avatar
Andreas Lauser committed
354
355
356

        // restore the orignal state of the element's secondary
        // variables
357
        curVolVars_[scvIdx] = origVolVars;
Andreas Lauser's avatar
Andreas Lauser committed
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374

#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
375
376
    {
        Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
Andreas Lauser's avatar
Andreas Lauser committed
377
378
379
380
381
382
383
384
385
386
        return 1e-9*(std::abs(pv) + 1);
    }

    /*!
     * \brief Updates the current local stiffness matrix with the
     *        partial derivatives of all equations in regard to the
     *        primary variable 'pvIdx' at vertex j .
     */
    void updateLocalJacobian_(int scvIdx,
                              int pvIdx,
387
                              const ElementSolutionVector &deriv)
Andreas Lauser's avatar
Andreas Lauser committed
388
    {
389
390
391
392
393
394
395
        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
396
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
397
398
399
400
401
                // 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
402
403
404
405
406
407
408
                Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
            }
        }
    }

    const Element *elemPtr_;
    FVElementGeometry fvElemGeom_;
409

Andreas Lauser's avatar
Andreas Lauser committed
410
411
412
413
    ElementBoundaryTypes bcTypes_;

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

Andreas Lauser's avatar
Andreas Lauser committed
415
416
    // secondary variables at the previous and at the current time
    // levels
417
418
419
    ElementVolumeVariables prevVolVars_;
    ElementVolumeVariables curVolVars_;

Andreas Lauser's avatar
Andreas Lauser committed
420
    LocalResidual localResidual_;
421

Andreas Lauser's avatar
Andreas Lauser committed
422
    LocalBlockMatrix A_;
423
424
    ElementSolutionVector residual_;
};
Andreas Lauser's avatar
Andreas Lauser committed
425
426
427
}

#endif