fvmpfaovelocity2p.hh 101 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:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
 *                                                                           *
6
 *   This program is free software: you can redistribute it and/or modify    *
7
 *   it under the terms of the GNU General Public License as published by    *
8
9
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
10
 *                                                                           *
11
12
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
 *   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/>.   *
18
19
 *****************************************************************************/

Andreas Lauser's avatar
Andreas Lauser committed
20
21
#ifndef DUMUX_MPFAOVELOCITY2P_HH
#define DUMUX_MPFAOVELOCITY2P_HH
22

23
#include "fvmpfaopressure2p.hh"
24
25
26

/**
 * @file
Markus Wolff's avatar
Markus Wolff committed
27
 * @brief  Velocity Field from a finite volume MPFA O-method solution of a pressure equation.
28
29
30
31
 */

namespace Dumux
{
Andreas Lauser's avatar
Andreas Lauser committed
32
template<class TypeTag>
Andreas Lauser's avatar
Andreas Lauser committed
33
34
class FVMPFAOPressure2P;

Markus Wolff's avatar
Markus Wolff committed
35
36
37
38
39
40
//! \ingroup FVPressure2p
/*! \brief Determines the velocity from a finite volume solution of the  pressure equation of a sequential model (IMPES).
 *  Calculates the total velocity from a known pressure field applying a finite volume discretization and a MPFA O-method.
 * The global pressure has to be given as piecewise constant cell values.
 * The velocities are calculated as
 * \f[
Thomas Fetzer's avatar
Thomas Fetzer committed
41
 * \boldsymbol v_{total} = \lambda_{total} \boldsymbol K \textbf{grad}\, p_{global}.
Markus Wolff's avatar
Markus Wolff committed
42
 * \f]
43
 *
Markus Wolff's avatar
Markus Wolff committed
44
45
46
47
 * Remark1: only for 2-D quadrilateral grids!
 * Remark2: can use UGGrid or SGrid/YaspGrid!
 * Remark3: gravity is neglected!
 *
Markus Wolff's avatar
Markus Wolff committed
48
 * \tparam TypeTag The problem Type Tag
49
 */
50
template<class TypeTag> class FVMPFAOVelocity2P:public FVMPFAOPressure2P<TypeTag>
51
{
52
53
    typedef FVMPFAOPressure2P<TypeTag> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
54
55
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
56

57
58
    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
    typedef typename SpatialParams::MaterialLaw MaterialLaw;
59

60
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
61

62
63
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
64

65
66
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
67
    typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
68
    typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
69

70
71
72
73
74
75
76
    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::Grid Grid;
    typedef typename GridView::IndexSet IndexSet;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::IntersectionIterator IntersectionIterator;
    typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;

77
    typedef typename GET_PROP_TYPE(TypeTag, GridTypeIndices) GridTypeIndices;
78
79
80
81
82
83
84

    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };
    enum
    {
85
86
        sw = Indices::saturationW,
        sn = Indices::saturationNw,
Markus Wolff's avatar
Markus Wolff committed
87
88
        vt = Indices::velocityTotal,
        pGlobal = Indices::pressureGlobal
89
90
91
    };
    enum
    {
92
93
94
95
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx,
        pressureIdx = Indices::pressureIdx,
        saturationIdx = Indices::saturationIdx,
96
        pressEqIdx = Indices::pressureEqIdx,
97
98
        satEqIdx = Indices::satEqIdx,
        numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
99
100
101
    };

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
102
103
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
    typedef Dune::FieldVector<Scalar, dim> DimVector;
104
105

public:
106
107
    //! Constructs a FVMPFAOVelocity2P object
    /*!
Markus Wolff's avatar
Markus Wolff committed
108
     * \param problem A problem class object
109
     */
110
    FVMPFAOVelocity2P(Problem& problem) :
111
        ParentType(problem), problem_(problem)
112
    {
Markus Wolff's avatar
Markus Wolff committed
113
114
115
116
117
118
119
120
        if (pressureType_ != pGlobal)
        {
            DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
        }
        if (dim != 2)
        {
            DUNE_THROW(Dune::NotImplemented, "MPFA method only implemented for 2-d!");
        }
121
122
    }

Markus Wolff's avatar
Markus Wolff committed
123
    //Calculates the velocities at all cell-cell interfaces.
124
125
    void calculateVelocity();

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
    void updateVelocity()
    {
        this->updateMaterialLaws();

        //reset velocities
        int size = problem_.gridView().size(0);
        for (int i = 0; i < size; i++)
        {
            CellData& cellData = problem_.variables().cellData(i);
            cellData.fluxData().resetVelocity();
        }

        calculateVelocity();
    }

Markus Wolff's avatar
Markus Wolff committed
141
142
143
144
    /*! \brief Initializes pressure and velocity
     *
     * \copydetails ParentType::initialize()
     */
145
146
147
148
    void initialize(bool solveTwice = true)
    {
        ParentType::initialize(solveTwice);

149
150
151
152
153
154
155
156
157
158
159
160
        const Element& element = *(problem_.gridView().template begin<0> ());
        FluidState fluidState;
        fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
        fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
        fluidState.setTemperature(problem_.temperature(element));
        fluidState.setSaturation(wPhaseIdx, 1.);
        fluidState.setSaturation(nPhaseIdx, 0.);
        density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
        density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
        viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
        viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);

161
162
163
164
165
        calculateVelocity();

        return;
    }

Markus Wolff's avatar
Markus Wolff committed
166
167
168
169
170
    /*! \brief Pressure and velocity update
     *
     * \copydetails ParentType::update()
     *
     */
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
    void update()
    {
        ParentType::update();

        calculateVelocity();

        return;
    }

private:
    Problem& problem_;

    Scalar density_[numPhases];
    Scalar viscosity_[numPhases];

    static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
    static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
188
189
}; // end of template

Markus Wolff's avatar
Markus Wolff committed
190
191
192
193
/*! \brief Calculates the velocities at all cell-cell interfaces.
*
* Calculates the velocities at all cell-cell interfaces from a given pressure field.
*/
194
195
196
197
template<class TypeTag>
void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
{
    // introduce matrix R for vector rotation and R is initialized as zero matrix
198
    DimMatrix R(0);
199
200
201
202
203
204
205
206
207

    // evaluate matrix R
    if (dim == 2)
        for (int i = 0; i < dim; ++i)
        {
            R[0][1] = 1;
            R[1][0] = -1;
        }

208
209
    BoundaryTypes bcType;

210
    // run through all elements
211
212
    ElementIterator eItEnd = problem_.gridView().template end<0> ();
    for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
213
214
215
216
    {
        // get common geometry information for the following computation

        // cell 1 geometry type
217
        //Dune::GeometryType gt1 = eIt->geometry().type();
218
219
220
221
222

        // get global coordinate of cell 1 center
        GlobalPosition globalPos1 = eIt->geometry().center();

        // cell 1 volume
223
        Scalar volume1 = eIt->geometry().volume();
224
225

        // cell 1 index
226
227
228
        int globalIdx1 = problem_.variables().index(*eIt);

        CellData& cellData1 = problem_.variables().cellData(globalIdx1);
229
230

        // get pressure value
231
        Scalar press1 = cellData1.globalPressure();
232
233

        // get right hand side
234
        PrimaryVariables source(0.0);
235
        problem_.source(source, *eIt);
236
        Scalar q1 = source[wPhaseIdx] + source[nPhaseIdx];
237
238

        // get absolute permeability of cell 1
239
        DimMatrix K1(problem_.spatialParams().intrinsicPermeability(*eIt));
240
241

        // compute total mobility of cell 1
242
243
        Scalar lambda1 = cellData1.mobility(wPhaseIdx)
                    + cellData1.mobility(nPhaseIdx);
244
245

        // the following two variables are used to check local conservation
246
        Scalar facevol[2 * dim];
247
        GlobalPosition unitOuterNormal[2 * dim];
248

249
250
        IntersectionIterator isItBegin = problem_.gridView().ibegin(*eIt);
        IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
251
252
253
254
255
256
257
258
259
        for (IntersectionIterator isIt = isItBegin; isIt != isItEnd; ++isIt)
        {
            // intersection iterator 'nextisIt' is used to get geometry information
            IntersectionIterator tempisIt = isIt;
            IntersectionIterator tempisItBegin = isItBegin;

            IntersectionIterator nextisIt = ++tempisIt;

            //get nextisIt
260
            switch (GET_PROP_VALUE(TypeTag, GridImplementation))
261
262
263
264
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
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
            {
                // for SGrid
                case GridTypeIndices::sGrid:
                {
                    if (nextisIt == isItEnd)
                    {
                        nextisIt = isItBegin;
                    }
                    else
                    {
                        nextisIt = ++tempisIt;

                        if (nextisIt == isItEnd)
                        {
                            nextisIt = ++tempisItBegin;
                        }
                    }

                    break;
                }
                // for YaspGrid
                case GridTypeIndices::yaspGrid:
                {
                    if (nextisIt == isItEnd)
                    {
                        nextisIt = isItBegin;
                    }
                    else
                    {
                        nextisIt = ++tempisIt;

                        if (nextisIt == isItEnd)
                        {
                            nextisIt = ++tempisItBegin;
                        }
                    }

                    break;
                }
                // for UGGrid
                case GridTypeIndices::ugGrid:
                {
                    if (nextisIt == isItEnd)
                    nextisIt = isItBegin;

                    break;
                }
                default:
                {
                    DUNE_THROW(Dune::NotImplemented, "GridType can not be used with MPFAO implementation!");
311
                    break;
312
313
314
315
                }
            }

            // get geometry type of face 'isIt', i.e., the face between cell1 and cell2 (locally numbered)
316
            //Dune::GeometryType gtf12 = isIt->geometryInInside().type();
317
318
319
320
321
322
323
324

            // center of face in global coordinates, i.e., the midpoint of edge 'isIt'
            GlobalPosition globalPosFace12 = isIt->geometry().center();

            // get local number of facet 'isIt'
            int indexInInside = isIt->indexInInside();

            // get face volume
325
            Scalar face12vol = isIt->geometry().volume();
326
327
328
329
330

            // get face volume to check if local mass conservative
            facevol[indexInInside] = isIt->geometry().volume();

            // get outer normal vector scaled with half volume of face 'isIt'
331
            GlobalPosition integrationOuterNormaln1 = isIt->centerUnitOuterNormal();
332
333
334
            integrationOuterNormaln1 *= face12vol / 2.0;

            // get unit outer normal vector of face 'isIt'
335
            const GlobalPosition& unitOuterNormaln1 = isIt->centerUnitOuterNormal();
336
337
338
339
340

            // get unit outer normal vector of face 'isIt' to check if local mass conservative
            unitOuterNormal[indexInInside] = unitOuterNormaln1;

            // get geometry type of 'nextisIt', i.e., face between cell1 and cell3 (locally numbered)
341
            //Dune::GeometryType gtf13 = nextisIt->geometryInInside().type();
342
343
344
345
346
347
348
349

            // center of face in global coordinates, i.e., the midpoint of edge 'nextisIt'
            GlobalPosition globalPosFace13 = nextisIt->geometry().center();

            // get local number of facet 'nextisIt'
            int nextindexInInside = nextisIt->indexInInside();

            // get face volume
350
            Scalar face13vol = nextisIt->geometry().volume();
351
352

            // get outer normal vector scaled with half volume of face 'nextisIt'
353
            GlobalPosition integrationOuterNormaln3 = nextisIt->centerUnitOuterNormal();
354
355
356
            integrationOuterNormaln3 *= face13vol / 2.0;

            // get unit outer normal vector of face 'nextisIt'
357
            const GlobalPosition& unitOuterNormaln3 = nextisIt->centerUnitOuterNormal();
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384

            // get the intersection node /bar^{x_3} between 'isIt' and 'nextisIt', denoted as 'corner1234'
            // initialization of corner1234
            GlobalPosition corner1234(0);

            // get the global coordinate of corner1234
            for (int i = 0; i < isIt->geometry().corners(); ++i)
            {
                GlobalPosition isItcorner = isIt->geometry().corner(i);

                for (int j = 0; j < nextisIt->geometry().corners(); ++j)
                {
                    GlobalPosition nextisItcorner = nextisIt->geometry().corner(j);

                    if (nextisItcorner == isItcorner)
                    {
                        corner1234 = isItcorner;
                        continue;
                    }
                }
            }

            // handle interior face
            if (isIt->neighbor())
            {
                // access neighbor cell 2 of 'isIt'
                ElementPointer outside = isIt->outside();
385
386
387
                int globalIdx2 = problem_.variables().index(*outside);

                CellData& cellData2 = problem_.variables().cellData(globalIdx2);
388
389

                // get pressure value
390
                Scalar press2 = cellData2.globalPressure();
391
392

                // neighbor cell 2 geometry type
393
                //Dune::GeometryType gt2 = outside->geometry().type();
394
395
396
397
398

                // get global coordinate of neighbor cell 2 center
                GlobalPosition globalPos2 = outside->geometry().center();

                // get absolute permeability of neighbor cell 2
399
                DimMatrix K2(problem_.spatialParams().intrinsicPermeability(*outside));
400
401

                // get total mobility of neighbor cell 2
402
403
                Scalar lambda2 = cellData2.mobility(wPhaseIdx)
                            + cellData2.mobility(nPhaseIdx);
404
405
406
407
408
409
410
411

                // 'nextisIt' is an interior face
                if (nextisIt->neighbor())
                {
                    // get basic information of cell 1,2's neighbor cell 3,4
                    // neighbor cell 3
                    // access neighbor cell 3
                    ElementPointer nextisItoutside = nextisIt->outside();
412
413
414
                    int globalIdx3 = problem_.variables().index(*nextisItoutside);

                    CellData& cellData3 = problem_.variables().cellData(globalIdx3);
415
416

                    // get pressure value
417
                    Scalar press3 = cellData3.globalPressure();
418
419

                    // neighbor cell 3 geometry type
420
                    //Dune::GeometryType gt3 = nextisItoutside->geometry().type();
421
422
423
424
425

                    // get global coordinate of neighbor cell 3 center
                    GlobalPosition globalPos3 = nextisItoutside->geometry().center();

                    // get absolute permeability of neighbor cell 3
426
                    DimMatrix K3(problem_.spatialParams().intrinsicPermeability(*nextisItoutside));
427
428

                    // get total mobility of neighbor cell 3
429
430
                    Scalar lambda3 = cellData3.mobility(wPhaseIdx)
                                + cellData3.mobility(nPhaseIdx);
431
432
433

                    // neighbor cell 4
                    GlobalPosition globalPos4(0);
434
                    DimMatrix K4(0);
435
                    Scalar lambda4 = 0;
436
                    int globalIdx4 = 0;
437
                    Scalar press4 = 0;
438

439
440
441
                    IntersectionIterator innerisItEnd = problem_.gridView().iend(*outside);
                    IntersectionIterator innernextisItEnd = problem_.gridView().iend(*nextisItoutside);
                    for (IntersectionIterator innerisIt = problem_.gridView().ibegin(*outside); innerisIt
442
                            != innerisItEnd; ++innerisIt)
443
                        for (IntersectionIterator innernextisIt = problem_.gridView().ibegin(
444
445
446
447
448
449
450
451
452
453
454
                                *nextisItoutside); innernextisIt != innernextisItEnd; ++innernextisIt)
                        {
                            if (innerisIt->neighbor() && innernextisIt->neighbor())
                            {
                                ElementPointer innerisItoutside = innerisIt->outside();
                                ElementPointer innernextisItoutside = innernextisIt->outside();

                                // find the common neighbor cell between cell 2 and cell 3, except cell 1
                                if (innerisItoutside == innernextisItoutside && innerisItoutside != isIt->inside())
                                {
                                    // access neighbor cell 4
455
456
457
                                    globalIdx4 = problem_.variables().index(*innerisItoutside);

                                    CellData& cellData4 = problem_.variables().cellData(globalIdx4);
458
459

                                    // neighbor cell 4 geometry type
460
                                    //Dune::GeometryType gt4 = innerisItoutside->geometry().type();
461
462
463
464
465

                                    // get global coordinate of neighbor cell 4 center
                                    globalPos4 = innerisItoutside->geometry().center();

                                    // get absolute permeability of neighbor cell 4
466
                                    K4 += problem_.spatialParams().intrinsicPermeability(*innerisItoutside);
467
468
469

                                    lambda4 = cellData4.mobility(wPhaseIdx)
                                                + cellData4.mobility(nPhaseIdx);
470

471
472
                                    // get pressure value
                                    press4 = cellData4.globalPressure();
473
474
475
476
                                }
                            }
                        }

477

478
479
480
481
482

                    // computation of flux through the first half edge of 'isIt' and the flux
                    // through the second half edge of 'nextisIt'

                    // get the information of the face 'isIt24' between cell2 and cell4 (locally numbered)
483
                    IntersectionIterator isIt24 = problem_.gridView().ibegin(*outside);
484

485
                    for (IntersectionIterator innerisIt = problem_.gridView().ibegin(*outside); innerisIt
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
                            != innerisItEnd; ++innerisIt)
                    {
                        if (innerisIt->neighbor())
                        {
                            if (innerisIt->outside() != isIt->inside())
                            {
                                for (int i = 0; i < innerisIt->geometry().corners(); ++i)
                                {
                                    GlobalPosition innerisItcorner = innerisIt->geometry().corner(i);

                                    if (innerisItcorner == corner1234)
                                    {
                                        isIt24 = innerisIt;
                                        continue;
                                    }
                                }
                            }
                        }
                    }

                    // get geometry type of face 'isIt24'
507
                    //Dune::GeometryType gtf24 = isIt24->geometryInInside().type();
508
509
510
511
512

                    // center of face in global coordinates, i.e., the midpoint of edge 'isIt24'
                    GlobalPosition globalPosFace24 = isIt24->geometry().center();

                    // get face volume
513
                    Scalar face24vol = isIt24->geometry().volume();
514
515

                    // get outer normal vector scaled with half volume of face 'isIt24'
516
                    GlobalPosition integrationOuterNormaln4 = isIt24->centerUnitOuterNormal();
517
518
519
                    integrationOuterNormaln4 *= face24vol / 2.0;

                    // get the information of the face 'isIt34' between cell3 and cell4 (locally numbered)
520
                    IntersectionIterator isIt34 = problem_.gridView().ibegin(*nextisItoutside);
521

522
                    for (IntersectionIterator innerisIt = problem_.gridView().ibegin(*nextisItoutside); innerisIt
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
                            != innernextisItEnd; ++innerisIt)
                    {
                        if (innerisIt->neighbor())
                        {
                            if (innerisIt->outside() != isIt->inside())
                            {
                                for (int i = 0; i < innerisIt->geometry().corners(); ++i)
                                {
                                    GlobalPosition innerisItcorner = innerisIt->geometry().corner(i);

                                    if (innerisItcorner == corner1234)
                                    {
                                        isIt34 = innerisIt;
                                        continue;
                                    }
                                }
                            }
                        }
                    }

                    // get geometry type of face 'isIt34'
544
                    //Dune::GeometryType gtf34 = isIt34->geometryInInside().type();
545
546
547
548
549

                    // center of face in global coordinates, i.e., the midpoint of edge 'isIt34'
                    GlobalPosition globalPosFace34 = isIt34->geometry().center();

                    // get face volume
550
                    Scalar face34vol = isIt34->geometry().volume();
551
552

                    // get outer normal vector scaled with half volume of face 'isIt34'
553
                    GlobalPosition integrationOuterNormaln2 = isIt34->centerUnitOuterNormal();
554
555
556
                    integrationOuterNormaln2 *= face34vol / 2.0;

                    // compute normal vectors nu11,nu21; nu12, nu22; nu13, nu23; nu14, nu24;
557
                    DimVector nu11(0);
558
559
                    R.umv(globalPosFace13 - globalPos1, nu11);

560
                    DimVector nu21(0);
561
562
                    R.umv(globalPos1 - globalPosFace12, nu21);

563
                    DimVector nu12(0);
564
565
                    R.umv(globalPosFace24 - globalPos2, nu12);

566
                    DimVector nu22(0);
567
568
                    R.umv(globalPosFace12 - globalPos2, nu22);

569
                    DimVector nu13(0);
570
571
                    R.umv(globalPos3 - globalPosFace13, nu13);

572
                    DimVector nu23(0);
573
574
                    R.umv(globalPos3 - globalPosFace34, nu23);

575
                    DimVector nu14(0);
576
577
                    R.umv(globalPos4 - globalPosFace24, nu14);

578
                    DimVector nu24(0);
579
580
581
                    R.umv(globalPosFace34 - globalPos4, nu24);

                    // compute dF1, dF2, dF3, dF4 i.e., the area of quadrilateral made by normal vectors 'nu'
582
                    DimVector Rnu21(0);
583
                    R.umv(nu21, Rnu21);
584
                    Scalar dF1 = fabs(nu11 * Rnu21);
585

586
                    DimVector Rnu22(0);
587
                    R.umv(nu22, Rnu22);
588
                    Scalar dF2 = fabs(nu12 * Rnu22);
589

590
                    DimVector Rnu23(0);
591
                    R.umv(nu23, Rnu23);
592
                    Scalar dF3 = fabs(nu13 * Rnu23);
593

594
                    DimVector Rnu24(0);
595
                    R.umv(nu24, Rnu24);
596
                    Scalar dF4 = fabs(nu14 * Rnu24);
597
598

                    // compute components needed for flux calculation, denoted as 'g'
599
                    DimVector K1nu11(0);
600
                    K1.umv(nu11, K1nu11);
601
                    DimVector K1nu21(0);
602
                    K1.umv(nu21, K1nu21);
603
                    DimVector K2nu12(0);
604
                    K2.umv(nu12, K2nu12);
605
                    DimVector K2nu22(0);
606
                    K2.umv(nu22, K2nu22);
607
                    DimVector K3nu13(0);
608
                    K3.umv(nu13, K3nu13);
609
                    DimVector K3nu23(0);
610
                    K3.umv(nu23, K3nu23);
611
                    DimVector K4nu14(0);
612
                    K4.umv(nu14, K4nu14);
613
                    DimVector K4nu24(0);
614
                    K4.umv(nu24, K4nu24);
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
                    Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1;
                    Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1;
                    Scalar g211 = lambda1 * (integrationOuterNormaln3 * K1nu11) / dF1;
                    Scalar g221 = lambda1 * (integrationOuterNormaln3 * K1nu21) / dF1;
                    Scalar g112 = lambda2 * (integrationOuterNormaln1 * K2nu12) / dF2;
                    Scalar g122 = lambda2 * (integrationOuterNormaln1 * K2nu22) / dF2;
                    Scalar g212 = lambda2 * (integrationOuterNormaln4 * K2nu12) / dF2;
                    Scalar g222 = lambda2 * (integrationOuterNormaln4 * K2nu22) / dF2;
                    Scalar g113 = lambda3 * (integrationOuterNormaln2 * K3nu13) / dF3;
                    Scalar g123 = lambda3 * (integrationOuterNormaln2 * K3nu23) / dF3;
                    Scalar g213 = lambda3 * (integrationOuterNormaln3 * K3nu13) / dF3;
                    Scalar g223 = lambda3 * (integrationOuterNormaln3 * K3nu23) / dF3;
                    Scalar g114 = lambda4 * (integrationOuterNormaln2 * K4nu14) / dF4;
                    Scalar g124 = lambda4 * (integrationOuterNormaln2 * K4nu24) / dF4;
                    Scalar g214 = lambda4 * (integrationOuterNormaln4 * K4nu14) / dF4;
                    Scalar g224 = lambda4 * (integrationOuterNormaln4 * K4nu24) / dF4;
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687

                    // compute transmissibility matrix T = CA^{-1}B+F
                    Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> C(0), F(0), A(0), B(0);

                    // evaluate matrix C, F, A, B
                    C[0][0] = -g111;
                    C[0][2] = -g121;
                    C[1][1] = g114;
                    C[1][3] = g124;
                    C[2][1] = -g213;
                    C[2][2] = g223;
                    C[3][0] = g212;
                    C[3][3] = -g222;

                    F[0][0] = g111 + g121;
                    F[1][3] = -g114 - g124;
                    F[2][2] = g213 - g223;
                    F[3][1] = -g212 + g222;

                    A[0][0] = g111 + g112;
                    A[0][2] = g121;
                    A[0][3] = -g122;
                    A[1][1] = g114 + g113;
                    A[1][2] = -g123;
                    A[1][3] = g124;
                    A[2][0] = g211;
                    A[2][1] = -g213;
                    A[2][2] = g223 + g221;
                    A[3][0] = -g212;
                    A[3][1] = g214;
                    A[3][3] = g222 + g224;

                    B[0][0] = g111 + g121;
                    B[0][1] = g112 - g122;
                    B[1][2] = g113 - g123;
                    B[1][3] = g114 + g124;
                    B[2][0] = g211 + g221;
                    B[2][2] = -g213 + g223;
                    B[3][1] = -g212 + g222;
                    B[3][3] = g214 + g224;

                    // compute T
                    A.invert();
                    F += B.leftmultiply(C.rightmultiply(A));
                    Dune::FieldMatrix<Scalar, 2 * dim, 2 * dim> T(F);

                    // use the pressure values to compute the fluxes
                    Dune::FieldVector<Scalar, 2 * dim> Tu(0);
                    Dune::FieldVector<Scalar, 2 * dim> u(0);
                    u[0] = press1;
                    u[1] = press2;
                    u[2] = press3;
                    u[3] = press4;

                    T.umv(u, Tu);

                    // evaluate velocity of facet 'isIt'
688
                    DimVector vector1 = unitOuterNormaln1;
689
                    vector1 *= Tu[0] / face12vol;
690
691
                    vector1 += cellData1.fluxData().velocityTotal(indexInInside);
                    cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1);
692
693

                    // evaluate velocity of facet 'nextisIt'
694
                    DimVector vector3 = unitOuterNormaln3;
695
                    vector3 *= Tu[2] / face13vol;
696
697
                    vector3 += cellData1.fluxData().velocityTotal(nextindexInInside);
                    cellData1.fluxData().setVelocity(wPhaseIdx, nextindexInInside, vector3);
698
699
700
701
702
703
704
705
706
707

                }
                // 'nextisIt' is on the boundary
                else
                {
                    // computation of flux through the first half edge of 'isIt' and the flux
                    // through the second half edge of 'nextisIt'

                    // get common geometry information for the following computation
                    // get the information of the face 'isIt24' between cell2 and cell4 (locally numbered)
708
709
710
                    IntersectionIterator isIt24 = problem_.gridView().ibegin(*outside);
                    IntersectionIterator innerisItEnd = problem_.gridView().iend(*outside);
                    for (IntersectionIterator innerisIt = problem_.gridView().ibegin(*outside); innerisIt
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
                            != innerisItEnd; ++innerisIt)
                    {
                        if (innerisIt->boundary())
                        {
                            for (int i = 0; i < innerisIt->geometry().corners(); ++i)
                            {
                                GlobalPosition innerisItcorner = innerisIt->geometry().corner(i);

                                if (innerisItcorner == corner1234)
                                {
                                    isIt24 = innerisIt;
                                    continue;
                                }
                            }
                        }
                    }

                    // get geometry type of face 'isIt24'
729
                    //Dune::GeometryType gtf24 = isIt24->geometryInInside().type();
730
731
732
733
734

                    // center of face in global coordinates, i.e., the midpoint of edge 'isIt24'
                    GlobalPosition globalPosFace24 = isIt24->geometry().center();

                    // get face volume
735
                    Scalar face24vol = isIt24->geometry().volume();
736
737

                    // get outer normal vector scaled with half volume of face 'isIt24'
738
                    GlobalPosition integrationOuterNormaln4 = isIt24->centerUnitOuterNormal();
739
740
                    integrationOuterNormaln4 *= face24vol / 2.0;

741
                    BoundaryTypes nextIsItBcType;
742
                    problem_.boundaryTypes(nextIsItBcType, *nextisIt);
743
744
745

                    // get boundary condition for boundary face (isIt24) center
                    BoundaryTypes isIt24BcType;
746
                    problem_.boundaryTypes(isIt24BcType, *isIt24);
747
748

                    PrimaryVariables boundValues(0.0);
749
750

                    // 'nextisIt': Neumann boundary
751
                    if (nextIsItBcType.isNeumann(pressEqIdx))
752
753
                    {
                        // get Neumann boundary value of 'nextisIt'
754
755
                        problem_.neumann(boundValues, *nextisIt);
                        Scalar J3 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]);
756
757

                        // 'isIt24': Neumann boundary
758
                        if (isIt24BcType.isNeumann(pressEqIdx))
759
760
                        {
                            // get neumann boundary value of 'isIt24'
761
762
                            problem_.neumann(boundValues, *isIt24);
                            Scalar J4 = (boundValues[wPhaseIdx]/density_[wPhaseIdx]+boundValues[nPhaseIdx]/density_[nPhaseIdx]);
763
764

                            // compute normal vectors nu11,nu21; nu12, nu22;
765
                            DimVector nu11(0);
766
767
                            R.umv(globalPosFace13 - globalPos1, nu11);

768
                            DimVector nu21(0);
769
770
                            R.umv(globalPos1 - globalPosFace12, nu21);

771
                            DimVector nu12(0);
772
773
                            R.umv(globalPosFace24 - globalPos2, nu12);

774
                            DimVector nu22(0);
775
776
777
                            R.umv(globalPosFace12 - globalPos2, nu22);

                            // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu'
778
                            DimVector Rnu21(0);
779
                            R.umv(nu21, Rnu21);
780
                            Scalar dF1 = fabs(nu11 * Rnu21);
781

782
                            DimVector Rnu22(0);
783
                            R.umv(nu22, Rnu22);
784
                            Scalar dF2 = fabs(nu12 * Rnu22);
785
786

                            // compute components needed for flux calculation, denoted as 'g'
787
                            DimVector K1nu11(0);
788
                            K1.umv(nu11, K1nu11);
789
                            DimVector K1nu21(0);
790
                            K1.umv(nu21, K1nu21);
791
                            DimVector K2nu12(0);
792
                            K2.umv(nu12, K2nu12);
793
                            DimVector K2nu22(0);
794
                            K2.umv(nu22, K2nu22);
795
796
797
798
799
800
801
802
                            Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1;
                            Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1;
                            Scalar g211 = lambda1 * (integrationOuterNormaln3 * K1nu11) / dF1;
                            Scalar g221 = lambda1 * (integrationOuterNormaln3 * K1nu21) / dF1;
                            Scalar g112 = lambda2 * (integrationOuterNormaln1 * K2nu12) / dF2;
                            Scalar g122 = lambda2 * (integrationOuterNormaln1 * K2nu22) / dF2;
                            Scalar g212 = lambda2 * (integrationOuterNormaln4 * K2nu12) / dF2;
                            Scalar g222 = lambda2 * (integrationOuterNormaln4 * K2nu22) / dF2;
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833

                            // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r
                            Dune::FieldMatrix<Scalar, 2 * dim - 1, 2 * dim - 1> A(0);
                            Dune::FieldMatrix<Scalar, 2 * dim - 1, dim> B(0);
                            Dune::FieldVector<Scalar, 2 * dim - 1> r1(0), r(0);

                            // evaluate matrix A, B
                            A[0][0] = g111 + g112;
                            A[0][1] = g121;
                            A[0][2] = -g122;
                            A[1][0] = g211;
                            A[1][1] = g221;
                            A[2][0] = -g212;
                            A[2][2] = g222;

                            B[0][0] = g111 + g121;
                            B[0][1] = g112 - g122;
                            B[1][0] = g211 + g221;
                            B[2][1] = g222 - g212;

                            // evaluate vector r1
                            r1[1] = -J3 * nextisIt->geometry().volume() / 2.0;
                            r1[2] = -J4 * isIt24->geometry().volume() / 2.0;

                            // compute T and r
                            A.invert();
                            B.leftmultiply(A);
                            Dune::FieldMatrix<Scalar, 2 * dim - 1, dim> T(B);
                            A.umv(r1, r);

                            // use the pressure values to compute the fluxes
834
                            Scalar f1 = (g111 + g121 - g111 * T[0][0] - g121 * T[1][0]) * press1 - (g111 * T[0][1]
835
836
837
                                    + g121 * T[1][1]) * press2 - (g111 * r[0] + g121 * r[1]);

                            // evaluate velocity of facet 'isIt'
838
                            DimVector vector1 = unitOuterNormaln1;
839
                            vector1 *= f1 / face12vol;
840
841
                            vector1 += cellData1.fluxData().velocityTotal(indexInInside);
                            cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1);
842
843
844

                        }
                        // 'isIt24': Dirichlet boundary
845
                        else if (isIt24BcType.isDirichlet(pressEqIdx))
846
847
                        {
                            // get Dirichlet boundary value on 'isIt24'
848
                            problem_.dirichlet(boundValues, *isIt24);
849
                             Scalar g4 = boundValues[pressureIdx];
850
851
852

                            // compute total mobility for Dirichlet boundary 'isIt24'
                            //determine lambda at the boundary -> if no saturation is known directly at the boundary use the cell saturation
853
                            Scalar alambda2 = 0;
854
                            if (isIt24BcType.isDirichlet(satEqIdx))
855
                            {
856
                                Scalar satBound = boundValues[saturationIdx];
857
858
859

                                //determine phase saturations from primary saturation variable
                                Scalar satW = 0;
860
                                switch (saturationType_)
861
                                {
862
                                case sw:
863
864
865
866
                                {
                                    satW = satBound;
                                    break;
                                }
867
                                case sn:
868
869
                                {
                                    satW = 1 - satBound;
870
                                    break;
871
872
873
874
875
876
877
                                }
                                }

                                Scalar lambdaWBound = 0;
                                Scalar lambdaNWBound = 0;

                                lambdaWBound = MaterialLaw::krw(
878
                                        problem_.spatialParams().materialLawParams(*eIt), satW)
879
                                        / viscosity_[wPhaseIdx];
880
                                lambdaNWBound = MaterialLaw::krn(
881
                                        problem_.spatialParams().materialLawParams(*eIt), satW)
882
                                        / viscosity_[nPhaseIdx];
883
884
885
886
887
888
889
890
                                alambda2 = lambdaWBound + lambdaNWBound;
                            }
                            else
                            {
                                alambda2 = lambda2;
                            }

                            // compute normal vectors nu11,nu21; nu12, nu22;
891
                            DimVector nu11(0);
892
893
                            R.umv(globalPosFace13 - globalPos1, nu11);

894
                            DimVector nu21(0);
895
896
                            R.umv(globalPos1 - globalPosFace12, nu21);

897
                            DimVector nu12(0);
898
899
                            R.umv(globalPosFace24 - globalPos2, nu12);

900
                            DimVector nu22(0);
901
902
903
                            R.umv(globalPosFace12 - globalPos2, nu22);

                            // compute dF1, dF2 i.e., the area of quadrilateral made by normal vectors 'nu'
904
                            DimVector Rnu21(0);
905
                            R.umv(nu21, Rnu21);
906
                            Scalar dF1 = fabs(nu11 * Rnu21);
907

908
                            DimVector Rnu22(0);
909
                            R.umv(nu22, Rnu22);
910
                            Scalar dF2 = fabs(nu12 * Rnu22);
911
912

                            // compute components needed for flux calculation, denoted as 'g'
913
                            DimVector K1nu11(0);
914
                            K1.umv(nu11, K1nu11);
915
                            DimVector K1nu21(0);
916
                            K1.umv(nu21, K1nu21);
917
                            DimVector K2nu12(0);
918
                            K2.umv(nu12, K2nu12);
919
                            DimVector K2nu22(0);
920
                            K2.umv(nu22, K2nu22);
921
922
923
924
925
926
                            Scalar g111 = lambda1 * (integrationOuterNormaln1 * K1nu11) / dF1;
                            Scalar g121 = lambda1 * (integrationOuterNormaln1 * K1nu21) / dF1;
                            Scalar g211 = lambda1 * (integrationOuterNormaln3 * K1nu11) / dF1;
                            Scalar g221 = lambda1 * (integrationOuterNormaln3 * K1nu21) / dF1;
                            Scalar g112 = alambda2 * (integrationOuterNormaln1 * K2nu12) / dF2;
                            Scalar g122 = alambda2 * (integrationOuterNormaln1 * K2nu22) / dF2;
927
928

                            // compute the matrix T & vector r in v = A^{-1}(Bu + r1) = Tu + r
929
930
                            DimMatrix A(0), B(0);
                            DimVector r1(0), r(0);
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948

                            // evaluate matrix A, B
                            A[0][0] = g111 + g112;
                            A[0][1] = g121;
                            A[1][0] = g211;
                            A[1][1] = g221;

                            B[0][0] = g111 + g121;
                            B[0][1] = g112 - g122;
                            B[1][0] = g211 + g221;

                            // evaluate vector r1
                            r1[0] = g122 * g4;
                            r1[1] = -J3 * nextisIt->geometry().volume() / 2.0;

                            // compute T and r
                            A.invert();
                            B.leftmultiply(A);
949
                            DimMatrix T(B);
950
951
952
                            A.umv(r1, r);

                            // use the pressure values to compute the fluxes
953
                            Scalar f1 = (g111 + g121 - g111 * T[0][0] - g121 * T[1][0]) * press1 - (g111 * T[0][1]
954
955
956
                                    + g121 * T[1][1]) * press2 - (g111 * r[0] + g121 * r[1]);

                            // evaluate velocity of facet 'isIt'
957
                            DimVector vector1 = unitOuterNormaln1;
958
                            vector1 *= f1 / face12vol;
959
960
                            vector1 += cellData1.fluxData().velocityTotal(indexInInside);
                            cellData1.fluxData().setVelocity(wPhaseIdx, indexInInside, vector1);
961
962
963
964

                        }
                    }
                    // 'nextisIt': Dirichlet boundary
965
                    else if (nextIsItBcType.isDirichlet(pressEqIdx))
966
967
                    {
                        // get Dirichlet boundary value of 'nextisIt'
968
                        problem_.dirichlet(boundValues, *nextisIt);
969
                        Scalar g3 = boundValues[pressureIdx];
970
971
972

                        // compute total mobility for Dirichlet boundary 'nextisIt'
                        //determine lambda at the boundary -> if no saturation is known directly at the boundary use the cell saturation
973
                        Scalar alambda1 = 0;
974
                        if (nextIsItBcType.isDirichlet(satEqIdx))
975
                        {
976
                            Scalar satBound = boundValues[saturationIdx];
977
978
979

                            //determine phase saturations from primary saturation variable
                            Scalar satW = 0;
980
                            switch (saturationType_)
981
                            {
982
                            case sw:
983
984
985
986
                            {
                                satW = satBound;
                                break;
                            }
987
                            case sn:
988
989
                            {
                                satW = 1 - satBound;
990
                                break;
991
992
993
994
995
996
997
                            }
                            }

                            Scalar lambdaWBound = 0;
                            Scalar lambdaNWBound = 0;

                            lambdaWBound = MaterialLaw::krw(
998
                                    problem_.spatialParams().materialLawParams(*eIt), satW)
999
                                    / viscosity_[wPhaseIdx];
1000
                            lambdaNWBound = MaterialLaw::krn(