fvmpfal2dpressurevelocity2padaptive.hh 29.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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.                                     *
 *                                                                           *
 *   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/>.   *
 *****************************************************************************/

#ifndef DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH
#define DUMUX_MPFAL2DPRESSUREVELOCITIES2P_ADAPTIVE_HH

21
22
#include <dune/common/float_cmp.hh>

23
#include "fvmpfal2dpressure2padaptive.hh"
24
#include "fvmpfal2dvelocity2padaptive.hh"
25
26
27
28
29
30
31
32
33
34

/**
 * @file
 * @brief  Velocity Field from a finite volume solution of a pressure equation using a grid adaptive MPFA L-method.
 */

namespace Dumux
{

//! \ingroup FVPressure2p
35
/*! \brief Class for the calculation of velocities from the  pressure solution of an IMPES scheme using a grid adaptive MPFA L-method.
36
 *
37
38
 * Can be used for calculating the complete velocity field before the solution of the transport equation (more efficient),
 * or for face-wise velocity calculation directly in the transport solution (less efficient).
39
40
41
 *
 * Remark1: only for 2-D quadrilateral grids!
 *
42
 * Remark2: can use UGGrid, ALUGrid
43
44
45
46
47
 *
 * Remark3: Allowed difference in grid levels of two neighboring cells: 1
 *
 * \tparam TypeTag The problem Type Tag
 */
48
template<class TypeTag> class FvMpfaL2dPressureVelocity2pAdaptive: public FvMpfaL2dPressure2pAdaptive<TypeTag>
49
50
51
52
53
54
55
56
57
58
59
60
{
    typedef FvMpfaL2dPressure2pAdaptive<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
61
62
63
64
    typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
    typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
65

66
67
    typedef typename Dune::ReferenceElements<Scalar, dim> ReferenceElements;
    typedef typename Dune::ReferenceElement<Scalar, dim> ReferenceElement;
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
    typedef typename SpatialParams::MaterialLaw MaterialLaw;

    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;

    typedef typename GridView::Grid Grid;
    typedef typename GridView::IndexSet IndexSet;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
    typedef typename GridView::IntersectionIterator IntersectionIterator;
    typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;
83
    typedef typename GridView::Intersection Intersection;
84

85
    typedef typename Grid::template Codim<0>::Entity::Geometry Geometry;
86
    typedef typename Geometry::JacobianTransposed JacobianTransposed;
87

88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    typedef typename ParentType::InteractionVolume InteractionVolume;

    enum
    {
        pw = Indices::pressureW,
        pn = Indices::pressureNw,
        sw = Indices::saturationW,
        sn = Indices::saturationNw
    };
    enum
    {
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx,
        pressureIdx = Indices::pressureIdx,
        saturationIdx = Indices::saturationIdx,
        pressEqIdx = Indices::pressureEqIdx,
        satEqIdx = Indices::satEqIdx,
        numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
    };

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
    typedef Dune::FieldVector<Scalar, dim> DimVector;

public:
113
    //! Constructs a FvMpfaL2dPressureVelocity2pAdaptive object
114
115
116
    /*!
     * \param problem A problem class object
     */
117
    FvMpfaL2dPressureVelocity2pAdaptive(Problem& problem) :
118
            ParentType(problem), problem_(problem), velocity_(problem)
119
120
121
122
123
124
    {
        density_[wPhaseIdx] = 0.;
        density_[nPhaseIdx] = 0.;
        viscosity_[wPhaseIdx] = 0.;
        viscosity_[nPhaseIdx] = 0.;

125
        calcVelocityInTransport_ = GET_PARAM_FROM_GROUP(TypeTag, bool, MPFA, CalcVelocityInTransport);
126
127
128
129
130
    }

    //Calculates the velocities at all cell-cell interfaces.
    void calculateVelocity();

131
132
133
134
    // Calculates the velocity at a cell-cell interface.
    void calculateVelocity(const Intersection&, CellData&);
    void calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData);

135
    //! Function for updating the velocity field if iterations are necessary in the transport solution
136
137
138
139
    void updateVelocity()
    {
        this->updateMaterialLaws();

140
        this->storePressureSolution();
141

142
143
        if (!calculateVelocityInTransport())
            calculateVelocity();
144
145
146
147
    }

    /*! \brief Initializes pressure and velocity
     *
Thomas Fetzer's avatar
Thomas Fetzer committed
148
     * \copydetails FVPressure::initialize()
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
     */
    void initialize()
    {
        ElementIterator 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);

164
165
166
167
        ParentType::initialize();
        velocity_.initialize();
        if (!calculateVelocityInTransport())
            calculateVelocity();
168
169
170
171
172
173

        return;
    }

    /*! \brief Pressure and velocity update
     *
Thomas Fetzer's avatar
Thomas Fetzer committed
174
     * \copydetails FVPressure::update()
175
176
177
178
179
     *
     */
    void update()
    {
        ParentType::update();
180
181
182
        if (!calculateVelocityInTransport())
            calculateVelocity();
    }
183

184
185
186
187
188
189
190
    /*! \brief Indicates if velocity is reconstructed in the pressure step or in the transport step
     *
     * Returns true (In the standard finite volume discretization the velocity is calculated during the saturation transport.)
     */
    bool calculateVelocityInTransport()
    {
        return calcVelocityInTransport_;
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
    }

    /*! \brief Adds velocity output to the output file
     *
     * Adds the phase velocities or a total velocity (depending on the formulation) to the output.
     * If the VtkOutputLevel is equal to zero (default) only primary variables are written,
     * if it is larger than zero also secondary variables are written.
     *
     * \tparam MultiWriter Class defining the output writer
     * \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
     *
     */
    template<class MultiWriter>
    void addOutputVtkFields(MultiWriter &writer)
    {
        ParentType::addOutputVtkFields(writer);
207
        velocity_.addOutputVtkFields(writer);
208
209
210
211
    }

private:
     Problem& problem_;
212
     FvMpfaL2dVelocity2pAdaptive<TypeTag> velocity_;
213
214
215

     Scalar density_[numPhases];
     Scalar viscosity_[numPhases];
216
     bool calcVelocityInTransport_;
217

218
     //! gives kind of pressure used (\f$ 0 = p_w\f$, \f$ 1 = p_n\f$, \f$ 2 = p_{global}\f$)
219
     static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
220
     //! gives kind of saturation used (\f$ 0 = S_w\f$, \f$ 1 = S_n\f$)
221
     static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
222
223
224
};
// end of template

225
/*! \brief Calculates the velocities at a cell-cell interfaces for the entire simulation grid.
226
 *
227
 * Calculates the velocities at a cell-cell interfaces from a given pressure field for the entire simulation grid.
228
229
230
 *
 */
template<class TypeTag>
231
void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity()
232
233
{
    // run through all elements
234
235
    VertexIterator vEndIt = problem_.gridView().template end<dim>();
    for (VertexIterator vIt = problem_.gridView().template begin<dim>(); vIt != vEndIt; ++vIt)
236
    {
237
        int vIdxGlobal = problem_.variables().index(*vIt);
238

239
        InteractionVolume& interactionVolume = this->interactionVolumes_[vIdxGlobal];
240
241
242
243
244
245
246
247
248
249
250

        if (interactionVolume.isInnerVolume())
        {
            if (interactionVolume.getElementNumber() == 4)
            {
                ElementPointer & elementPointer1 = interactionVolume.getSubVolumeElement(0);
                ElementPointer & elementPointer2 = interactionVolume.getSubVolumeElement(1);
                ElementPointer & elementPointer3 = interactionVolume.getSubVolumeElement(2);
                ElementPointer & elementPointer4 = interactionVolume.getSubVolumeElement(3);

                // cell index
251
252
253
254
                int eIdxGlobal1 = problem_.variables().index(*elementPointer1);
                int eIdxGlobal2 = problem_.variables().index(*elementPointer2);
                int eIdxGlobal3 = problem_.variables().index(*elementPointer3);
                int eIdxGlobal4 = problem_.variables().index(*elementPointer4);
255
256

                //get the cell Data
257
258
259
260
                CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
                CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
                CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
                CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
261

262
263
                velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellData1, cellData2, cellData3,
                                                                  cellData4, this->innerBoundaryVolumeFaces_);
264

265
266
267
268
269
270
            }
            else if (interactionVolume.getElementNumber() == 3)
            {
                ElementPointer & elementPointer1 = interactionVolume.getSubVolumeElement(0);
                  ElementPointer & elementPointer2 = interactionVolume.getSubVolumeElement(1);
                  ElementPointer & elementPointer4 = interactionVolume.getSubVolumeElement(3);
271

272
                  // cell index
273
274
275
                  int eIdxGlobal1 = problem_.variables().index(*elementPointer1);
                  int eIdxGlobal2 = problem_.variables().index(*elementPointer2);
                  int eIdxGlobal4 = problem_.variables().index(*elementPointer4);
276

277
                  //get the cell Data
278
279
280
                  CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
                  CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
                  CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
281

282
283
                  velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume, cellData1, cellData2,
                                                                          cellData4, this->innerBoundaryVolumeFaces_);
284

285
286
287
288
289
290
            }
            else
            {
                DUNE_THROW(Dune::NotImplemented, "Unknown interactionvolume type!");
            }
        }
291

292
293
294
295
296
297
        // at least one face on boundary!
        else
        {
            for (int elemIdx = 0; elemIdx < 2 * dim; elemIdx++)
            {
                bool isOutside = false;
298
                for (int fIdx = 0; fIdx < dim; fIdx++)
299
                {
300
                    int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
301
302
303
304
305
                    if (interactionVolume.isOutsideFace(intVolFaceIdx))
                    {
                        isOutside = true;
                        break;
                    }
306
                }
307
                if (isOutside)
308
                {
309
                    continue;
310
                }
311
                ElementPointer & elementPointer = interactionVolume.getSubVolumeElement(elemIdx);
312

313
                // cell index
314
                int eIdxGlobal = problem_.variables().index(*elementPointer);
315
                //get the cell Data
316
                CellData& cellData = problem_.variables().cellData(eIdxGlobal);
317

318
319
                velocity_.calculateBoundaryInteractionVolumeVelocity(interactionVolume, cellData, elemIdx);
            }
320

321
        } // end boundaries
322

323
324
325
    } // end vertex iterator
    return;
} // end method calcTotalVelocity
326

327
328
329
330
331
332
333
/*! \brief Calculates the velocity at a cell-cell interface.
 *
 * Calculates the velocity at a cell-cell interface from a given pressure field.
 *
 * \param intersection Intersection of two grid cells
 * \param cellData Object containing all model relevant cell data
 */
334
335
336
337
template<class TypeTag>
void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellData)
{
    int numVertices = intersection.geometry().corners();
338

339
340
    ElementPointer elementPtrI = intersection.inside();
    ElementPointer elementPtrJ = intersection.outside();
341

342
343
    int levelI = elementPtrI->level();
    int levelJ = elementPtrJ->level();
344

345
346
    int eIdxGlobalI = problem_.variables().index(*elementPtrI);
    int eIdxGlobalJ = problem_.variables().index(*elementPtrJ);
347

348
    CellData& cellDataJ = problem_.variables().cellData(eIdxGlobalJ);
349

350
    const ReferenceElement& referenceElement = ReferenceElements::general(elementPtrI->geometry().type());
351

352
353
    int indexInInside = intersection.indexInInside();
    int indexInOutside = intersection.indexInOutside();
354

355
    int fIdx = indexInInside;
356

357
    if (levelI < levelJ)
358
        fIdx = indexInOutside;
359

360
    std::vector<CellData> cellDataTemp(0);
361

362
363
364
365
366
367
368
369
370
371
372
373
374
    if (levelI != levelJ)
    {
        DimVector vel(0);
        cellData.fluxData().setVelocity(wPhaseIdx, indexInInside, vel);
        cellData.fluxData().setVelocity(nPhaseIdx, indexInInside, vel);
        cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside, 0);
        cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside, 0);

        cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside, vel);
        cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside, vel);
        cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside, 0);
        cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside, 0);
    }
375

376
377
    for (int vIdx = 0; vIdx < numVertices; vIdx++)
    {
378
        int localVertIdx = referenceElement.subEntity(fIdx, dim - 1, vIdx, dim);
379

380
        int vIdxGlobal = 0;
381
                if (levelI >= levelJ)
382
                {
383
                    vIdxGlobal = problem_.variables().index(
384
                            *((*elementPtrI).template subEntity < dim > (localVertIdx)));
385
                }
386
                else
387
                {
388
                    vIdxGlobal = problem_.variables().index(
389
                            *((*elementPtrJ).template subEntity < dim > (localVertIdx)));
390
391
                }

392
        InteractionVolume& interactionVolume = this->interactionVolumes_[vIdxGlobal];
393

394
395
396
        if (interactionVolume.isInnerVolume())
        {
            // cell index vector
397
            std::vector<int> eIdxGlobal(0);
398

399
400
401
402
403
404
            if (interactionVolume.getElementNumber() == 4)
            {
                ElementPointer & elementPointer1 = interactionVolume.getSubVolumeElement(0);
                ElementPointer & elementPointer2 = interactionVolume.getSubVolumeElement(1);
                ElementPointer & elementPointer3 = interactionVolume.getSubVolumeElement(2);
                ElementPointer & elementPointer4 = interactionVolume.getSubVolumeElement(3);
405

406
                eIdxGlobal.resize(4);
407

408
409
410
411
                eIdxGlobal[0] = problem_.variables().index(*elementPointer1);
                eIdxGlobal[1] = problem_.variables().index(*elementPointer2);
                eIdxGlobal[2] = problem_.variables().index(*elementPointer3);
                eIdxGlobal[3] = problem_.variables().index(*elementPointer4);
412

413
414
                //cell Data vector
                cellDataTemp.resize(4);
415

416
417
418
419
                cellDataTemp[0] = problem_.variables().cellData(eIdxGlobal[0]);
                cellDataTemp[1] = problem_.variables().cellData(eIdxGlobal[1]);
                cellDataTemp[2] = problem_.variables().cellData(eIdxGlobal[2]);
                cellDataTemp[3] = problem_.variables().cellData(eIdxGlobal[3]);
420

421
422
                velocity_.calculateInnerInteractionVolumeVelocity(interactionVolume, cellDataTemp[0], cellDataTemp[1],
                                                                  cellDataTemp[2], cellDataTemp[3], this->innerBoundaryVolumeFaces_);
423
424
425
426
427
428
429
            }
            else if (interactionVolume.getElementNumber() == 3)
            {
                ElementPointer & elementPointer1 = interactionVolume.getSubVolumeElement(0);
                ElementPointer & elementPointer2 = interactionVolume.getSubVolumeElement(1);
                ElementPointer & elementPointer4 = interactionVolume.getSubVolumeElement(3);

430
                eIdxGlobal.resize(3);
431

432
433
434
                eIdxGlobal[0] = problem_.variables().index(*elementPointer1);
                eIdxGlobal[1] = problem_.variables().index(*elementPointer2);
                eIdxGlobal[2] = problem_.variables().index(*elementPointer4);
435

436
437
                //cell Data vector
                cellDataTemp.resize(3);
438

439
440
441
                cellDataTemp[0] = problem_.variables().cellData(eIdxGlobal[0]);
                cellDataTemp[1] = problem_.variables().cellData(eIdxGlobal[1]);
                cellDataTemp[2] = problem_.variables().cellData(eIdxGlobal[2]);
442
443


444
445
                velocity_.calculateHangingNodeInteractionVolumeVelocity(interactionVolume, cellDataTemp[0], cellDataTemp[1],
                                                                        cellDataTemp[2], this->innerBoundaryVolumeFaces_);
446
447
448
449
450
            }
            else
            {
                DUNE_THROW(Dune::NotImplemented, "Unknown interactionvolume type!");
            }
451

452
453
454
            int size = cellDataTemp.size();
        for (int i = 0; i < size; i++)
        {
455
            if (eIdxGlobal[i] == eIdxGlobalI)
456
457
458
            {
                if (levelI >= levelJ)
                {
459
460
461
462
463
464
465
466
                 cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
                                                 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInInside));
                 cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
                                                 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInInside));
                 cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
                                                        cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInInside));
                 cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
                                                        cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInInside));
467
468
469

                 if (levelI > levelJ)
                 {
470
471
472
473
474
475
476
477
                     cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
                                                      cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInInside));
                     cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
                                                      cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInInside));
                     cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
                                                             cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInInside));
                     cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
                                                             cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInInside));
478
479
480

                 }
                }
481

482
            }
483
            else if (eIdxGlobal[i] == eIdxGlobalJ)
484
485
486
            {
                if (levelJ >= levelI)
                {
487
488
489
490
491
492
493
494
                cellDataJ.fluxData().setVelocity(wPhaseIdx, indexInOutside,
                                                 cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInOutside));
                cellDataJ.fluxData().setVelocity(nPhaseIdx, indexInOutside,
                                                 cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInOutside));
                cellDataJ.fluxData().setUpwindPotential(wPhaseIdx, indexInOutside,
                                                        cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
                cellDataJ.fluxData().setUpwindPotential(nPhaseIdx, indexInOutside,
                                                        cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
495

496
497
                if (levelJ > levelI)
                {
498
499
500
501
502
503
504
505
                    cellData.fluxData().setVelocity(wPhaseIdx, indexInInside,
                                                    cellDataTemp[i].fluxData().velocity(wPhaseIdx, indexInOutside));
                    cellData.fluxData().setVelocity(nPhaseIdx, indexInInside,
                                                    cellDataTemp[i].fluxData().velocity(nPhaseIdx, indexInOutside));
                    cellData.fluxData().setUpwindPotential(wPhaseIdx, indexInInside,
                                                           cellDataTemp[i].fluxData().upwindPotential(wPhaseIdx, indexInOutside));
                    cellData.fluxData().setUpwindPotential(nPhaseIdx, indexInInside,
                                                           cellDataTemp[i].fluxData().upwindPotential(nPhaseIdx, indexInOutside));
506

507
508
509
510
511
512
513
514
515
516
517
518
                }
                }
            }
        }
    }
    }
    if (levelI == levelJ)
    {
        cellData.fluxData().setVelocityMarker(indexInInside);
        cellDataJ.fluxData().setVelocityMarker(indexInOutside);
    }
}
519

520
521
522
523
524
525
526
/*! \brief Calculates the velocity at a boundary.
 *
 * Calculates the velocity at a boundary from a given pressure field.
 *
 * \param intersection Boundary intersection
 * \param cellData Object containing all model relevant cell data
 */
527
528
529
530
template<class TypeTag>
void FvMpfaL2dPressureVelocity2pAdaptive<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
{
    ElementPointer element = intersection.inside();
531

532
533
    //get face index
    int isIndex = intersection.indexInInside();
534

535
536
    //get face normal
    const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
537

538
539
540
541
    BoundaryTypes bcType;
    //get boundary type
    problem_.boundaryTypes(bcType, intersection);
    PrimaryVariables boundValues(0.0);
542

543
544
545
    if (bcType.isDirichlet(pressEqIdx))
    {
        problem_.dirichlet(boundValues, intersection);
546

547
548
        // get global coordinates of cell centers
        const GlobalPosition& globalPosI = element->geometry().center();
549

550
551
        // center of face in global coordinates
        const GlobalPosition& globalPosJ = intersection.geometry().center();
552

553
554
555
        // get mobilities and fractional flow factors
        Scalar lambdaWI = cellData.mobility(wPhaseIdx);
        Scalar lambdaNwI = cellData.mobility(nPhaseIdx);
556

557
558
        // get capillary pressure
        Scalar pcI = cellData.capillaryPressure();
559

560
561
        // distance vector between barycenters
        GlobalPosition distVec = globalPosJ - globalPosI;
562

563
564
        // compute distance between cell centers
        Scalar dist = distVec.two_norm();
565

566
567
568
        //permeability vector at boundary
        // compute vectorized permeabilities
        DimMatrix meanPermeability(0);
569

570
        problem_.spatialParams().meanK(meanPermeability, problem_.spatialParams().intrinsicPermeability(*element));
571

572
573
        Dune::FieldVector<Scalar, dim> permeability(0);
        meanPermeability.mv(unitOuterNormal, permeability);
574

575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
        //determine saturation at the boundary -> if no saturation is known directly at the boundary use the cell saturation
        Scalar satW = 0;
        if (bcType.isDirichlet(satEqIdx))
        {
            switch (saturationType_)
            {
            case sw:
            {
                satW = boundValues[saturationIdx];
                break;
            }
            case sn:
            {
                satW = 1 - boundValues[saturationIdx];
                break;
            }
            }
        }
        else
        {
            satW = cellData.saturation(wPhaseIdx);
        }
597

598
599
        Scalar pressBound = boundValues[pressureIdx];
        Scalar pcBound = MaterialLaw::pc(problem_.spatialParams().materialLawParams(*element), satW);
600

601
602
603
604
605
606
607
608
609
610
611
612
613
        //determine phase pressures from primary pressure variable
        Scalar pressWBound = 0;
        Scalar pressNwBound = 0;
        if (pressureType_ == pw)
        {
            pressWBound = pressBound;
            pressNwBound = pressBound + pcBound;
        }
        else if (pressureType_ == pn)
        {
            pressWBound = pressBound - pcBound;
            pressNwBound = pressBound;
        }
614

615
616
617
618
        Scalar lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satW)
                / viscosity_[wPhaseIdx];
        Scalar lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satW)
                / viscosity_[nPhaseIdx];
619

620
621
        Scalar potentialDiffW = cellData.fluxData().upwindPotential(wPhaseIdx, isIndex);
        Scalar potentialDiffNw = cellData.fluxData().upwindPotential(nPhaseIdx, isIndex);
622

623
624
625
        //calculate potential gradient
        potentialDiffW = (cellData.pressure(wPhaseIdx) - pressWBound);
        potentialDiffNw = (cellData.pressure(nPhaseIdx) - pressNwBound);
626

627
628
        potentialDiffW += density_[wPhaseIdx] * (distVec * problem_.gravity());
        potentialDiffNw += density_[nPhaseIdx] * (distVec * problem_.gravity());
629

630
631
632
        //store potential gradients for further calculations
        cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, potentialDiffW);
        cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, potentialDiffNw);
633

634
635
        //do the upwinding of the mobility depending on the phase potentials
        Scalar lambdaW = (potentialDiffW > 0.) ? lambdaWI : lambdaWBound;
636
        lambdaW = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffW, 0.0, 1.0e-30)) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
637
        Scalar lambdaNw = (potentialDiffNw > 0.) ? lambdaNwI : lambdaNwBound;
638
        lambdaNw = (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potentialDiffNw, 0.0, 1.0e-30)) ? 0.5 * (lambdaNwI + lambdaNwBound) : lambdaNw;
639
640


641
        Scalar scalarPerm = permeability.two_norm();
642

643
644
645
        //calculate the gravity term
        Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
        Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
646

647
648
649
        //calculate unit distVec
        distVec /= dist;
        Scalar areaScaling = (unitOuterNormal * distVec);
650
651
        //this treatment of g allows to account for gravity flux through faces where the face normal
        //has no z component (e.g. parallelepiped grids)
652
653
        Scalar gravityTermW = (problem_.gravity() * distVec) * density_[wPhaseIdx] * areaScaling;
        Scalar gravityTermNw = (problem_.gravity() * distVec) * density_[nPhaseIdx] * areaScaling;
654

655
656
657
658
659
660
661
662
663
        //calculate velocity depending on the pressure used -> use pc = pn - pw
        switch (pressureType_)
        {
        case pw:
        {
            velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
            velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNw)
                    + 0.5 * (lambdaNwI + lambdaNwBound) * scalarPerm * (pcI - pcBound) / dist;
            break;
664
        }
665
        case pn:
666
        {
667
668
669
670
671
672
            velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
                    - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
            velocityNw *= lambdaNw * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNw);
            break;
        }
        }
673

674
675
676
677
        //store velocities
        cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
        cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
        cellData.fluxData().setVelocityMarker(isIndex);
678

679
    } //end dirichlet boundary
680

681
682
683
    else if (bcType.isNeumann(pressEqIdx))
    {
        problem_.neumann(boundValues, intersection);
684

685
686
        Dune::FieldVector<Scalar, dimWorld> velocityW(unitOuterNormal);
        Dune::FieldVector<Scalar, dimWorld> velocityNw(unitOuterNormal);
687

688
689
        velocityW *= boundValues[wPhaseIdx];
        velocityNw *= boundValues[nPhaseIdx];
690

691
692
            velocityW /= density_[wPhaseIdx];
            velocityNw /= density_[nPhaseIdx];
693

694
695
696
        //store potential gradients for further calculations
        cellData.fluxData().setUpwindPotential(wPhaseIdx, isIndex, boundValues[wPhaseIdx]);
        cellData.fluxData().setUpwindPotential(nPhaseIdx, isIndex, boundValues[nPhaseIdx]);
697

698
699
700
701
702
703
704
705
706
        cellData.fluxData().setVelocity(wPhaseIdx, isIndex, velocityW);
        cellData.fluxData().setVelocity(nPhaseIdx, isIndex, velocityNw);
        cellData.fluxData().setVelocityMarker(isIndex);
    } //end neumann boundary
    else
    {
        DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
    }
}
707
708
709
710

}
// end of Dune namespace
#endif