boxlocaljacobian.hh 12.3 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
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// $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
 * \brief Caculates the jacobian of models based on the box scheme element-wise.
 */
#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
#define DUMUX_BOX_LOCAL_JACOBIAN_HH

#include <dune/common/exceptions.hh>

#include <dumux/common/valgrind.hh>

#include <dune/grid/common/genericreferenceelements.hh>

#include <boost/format.hpp>

#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>

40
#include "boxelementvolumevariables.hh"
Andreas Lauser's avatar
Andreas Lauser committed
41
42
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
#include "boxfvelementgeometry.hh"
#include "boxlocalresidual.hh"

#include "boxproperties.hh"



namespace Dumux
{
/*!
 * \ingroup BoxModel
 * \brief Element-wise caculation of the jacobian matrix for models
 *        based on the box scheme .
 *
 * \todo Please doc me more!
 */
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;

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

71
72
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
Andreas Lauser's avatar
Andreas Lauser committed
73
74
75
    };

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
76
    typedef typename GridView::Grid::ctype CoordScalar;
Andreas Lauser's avatar
Andreas Lauser committed
77
78
79
80

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

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

    typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
86
87
    typedef typename RefElemProp::Container ReferenceElements;
    typedef typename RefElemProp::ReferenceElement ReferenceElement;
Andreas Lauser's avatar
Andreas Lauser committed
88

89
90
91
    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
92
93
94
95

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
96
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
Andreas Lauser's avatar
Andreas Lauser committed
97

98
99
    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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;

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

public:
    BoxLocalJacobian()
    {
        Valgrind::SetUndefined(problemPtr_);
    }


    void init(Problem &prob)
    {
        problemPtr_ = &prob;
        localResidual_.init(prob);
        // assume quadrilinears as elements with most vertices
        A_.setSize(2<<dim, 2<<dim);
    }
119

Andreas Lauser's avatar
Andreas Lauser committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
    /*!
     * \brief Assemble the linear system of equations for the
     *        verts of a element, given a local solution 'localU'.
     */
    void assemble(const Element &element)
    {
        // set the current grid element and update the element's
        // finite volume geometry
        elemPtr_ = &element;
        fvElemGeom_.update(gridView_(), elem_());
        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
144
        prevVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
145
146
147
                            elem_(),
                            fvElemGeom_,
                            true /* isOldSol? */);
148
        curVolVars_.update(problem_(),
Andreas Lauser's avatar
Andreas Lauser committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
                           elem_(),
                           fvElemGeom_,
                           false /* isOldSol? */);

        // 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);
            }
        }
    }

    /*!
     * \brief Returns a reference to the local residual.
     */
    const LocalResidual &localResidual() const
    { return localResidual_; }

    /*!
     * \brief Returns a reference to the local residual.
     */
179
    LocalResidual &localResidual()
Andreas Lauser's avatar
Andreas Lauser committed
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
    { return localResidual_; }

    /*!
     * \brief Returns the jacobian of the equations at vertex i to the
     *        primary variables at vertex j.
     */
    const MatrixBlock &mat(int i, int j) const
    { return A_[i][j]; }

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
199
    {
Andreas Lauser's avatar
Andreas Lauser committed
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
        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(); };

    /*!
     * \brief Returns a reference to the vertex mapper.
     */
    const VertexMapper &vertexMapper_() const
    { return problem_().vertexMapper(); };

    /*!
     * \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);
243
244
        PrimaryVariables priVars(model_().curSol()[globalIdx]);
        VolumeVariables origVolVars(curVolVars_[scvIdx]);
Andreas Lauser's avatar
Andreas Lauser committed
245

246
        curVolVars_[scvIdx].setEvalPoint(&origVolVars);
Andreas Lauser's avatar
Andreas Lauser committed
247
        Scalar eps = asImp_().numericEpsilon_(scvIdx, pvIdx);
248

Andreas Lauser's avatar
Andreas Lauser committed
249
250
        // deflect primary variables
        priVars[pvIdx] += eps;
251

Andreas Lauser's avatar
Andreas Lauser committed
252
        // calculate the residual
253
        curVolVars_[scvIdx].update(priVars,
Andreas Lauser's avatar
Andreas Lauser committed
254
255
256
257
258
                                   problem_(),
                                   elem_(),
                                   fvElemGeom_,
                                   scvIdx,
                                   false);
259
260
261
262
        localResidual().eval(elem_(),
                             fvElemGeom_,
                             prevVolVars_,
                             curVolVars_,
Andreas Lauser's avatar
Andreas Lauser committed
263
264
265
266
267
268
269
270
271
                             bcTypes_);

        // store the residual
        dest = localResidual().residual();

        // deflect the primary variables
        priVars[pvIdx] -= 2*eps;

        // calculate residual again
272
        curVolVars_[scvIdx].update(priVars,
Andreas Lauser's avatar
Andreas Lauser committed
273
274
275
276
277
                                   problem_(),
                                   elem_(),
                                   fvElemGeom_,
                                   scvIdx,
                                   false);
278
279
280
281
        localResidual().eval(elem_(),
                             fvElemGeom_,
                             prevVolVars_,
                             curVolVars_,
Andreas Lauser's avatar
Andreas Lauser committed
282
283
284
285
286
287
288
289
                             bcTypes_);

        // central differences
        dest -= localResidual().residual();
        dest /= 2*eps;

        // restore the orignal state of the element's secondary
        // variables
290
        curVolVars_[scvIdx] = origVolVars;
Andreas Lauser's avatar
Andreas Lauser committed
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308

#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 elemSol    The current solution on the element
     * \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
309
310
    {
        Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
Andreas Lauser's avatar
Andreas Lauser committed
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
        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,
                              const ElementSolutionVector &stiffness)
    {
        for (int i = 0; i < fvElemGeom_.numVertices; i++) {
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
                // A[i][scvIdx][eqIdx][pvIdx] is the approximate 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] = stiffness[i][eqIdx];
                Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
            }
        }
    }

    const Element *elemPtr_;

    FVElementGeometry fvElemGeom_;
    ElementBoundaryTypes bcTypes_;

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

Andreas Lauser's avatar
Andreas Lauser committed
343
344
    // secondary variables at the previous and at the current time
    // levels
345
346
347
    ElementVolumeVariables prevVolVars_;
    ElementVolumeVariables curVolVars_;

Andreas Lauser's avatar
Andreas Lauser committed
348
349
350
351
352
353
    LocalResidual localResidual_;
    LocalBlockMatrix A_;
};
}

#endif