evalcflfluxcoats.hh 24.2 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
20
21
22
23
 *****************************************************************************/
#ifndef DUMUX_EVALCFLFLUX_COATS_HH
#define DUMUX_EVALCFLFLUX_COATS_HH

/**
 * @file
24
 * @brief  Cfl-flux-function to evaluate a Cfl-Condition after Coats 2003
25
26
 */

27
#include <dumux/decoupled/common/impetproperties.hh> 
28
#include "evalcflflux.hh"
29
30
31

namespace Dumux
{
Markus Wolff's avatar
Markus Wolff committed
32
/*!\ingroup IMPES
33
 * @brief  Cfl-flux-function to evaluate a Cfl-Condition after Coats 2003
34
 *
Markus Wolff's avatar
Markus Wolff committed
35
 * tparam TypeTag The problem TypeTag
36
37
 */
template<class TypeTag>
38
class EvalCflFluxCoats: public EvalCflFlux<TypeTag>
39
40
{
private:
41
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
42
43
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
44

45
46
    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
    typedef typename SpatialParams::MaterialLaw MaterialLaw;
47
48
49
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
50

51
52
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
53
    typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
54

55
    typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
56

57
    enum
58
59
60
        {
            dim = GridView::dimension, dimWorld = GridView::dimensionworld
        };
61
    enum
62
63
        {
            wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
64
            eqIdxPress = Indices::pressureEqIdx,
65
66
67
            eqIdxSat = Indices::satEqIdx,
            numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
        };
68
69

    enum
70
71
        {
            pw = Indices::pressureW,
72
            pn = Indices::pressureNw,
73
            vt = Indices::velocityTotal,
74
75
            sw = Indices::saturationW,
            sn = Indices::saturationNw
76
        };
77
78
79

    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
80
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
81
82
83
    typedef typename GridView::Intersection Intersection;

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
84
    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
85
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
86
87

public:
88
89
90
91
92
93
94
95
96
97
98
99
100
    //! \brief Initializes the cfl-flux-model
    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);
    }
101

Markus Wolff's avatar
Markus Wolff committed
102
103
104
105
    /*! \brief adds a flux to the cfl-criterion evaluation
     *
     * \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Element&,int)
     */
106
    void addFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
107
                 const Element& element, int phaseIdx = -1)
108
    {
109
        addDefaultFlux(flux, phaseIdx);
110
111
    }

Markus Wolff's avatar
Markus Wolff committed
112
113
114
115
    /*! \brief adds a flux to the cfl-criterion evaluation
     *
     * \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Intersection&,int)
     */
116
    void addFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
117
                 const Intersection& intersection, int phaseIdx = -1)
118
    {
119
        addDefaultFlux(flux, phaseIdx);
120
        addCoatsFlux(lambdaW, lambdaNw, viscosityW, viscosityNw, flux, intersection, phaseIdx);
121
    }
122

123
    /*! \brief Returns the Cfl flux-function
124
     *
125
     * \copydetails EvalCflFlux::getCflFluxFunction(const Element&)
126
     */
127
    Scalar getCflFluxFunction(const Element& element)
128
    {
129

130
131
132
//        Scalar cflFluxDefault = getCflFluxFunctionDefault();
//
//    	return 0.99 / cflFluxDefault;
133

134
135
136
137
138
139
140
141
142
143
144
145
    	 Scalar cflFluxDefault = getCflFluxFunctionDefault();

        if (rejectForTimeStepping_)
        	return 0.99 / cflFluxDefault;

        //        std::cout<<"cflFluxFunctionCoats_ = "<<cflFluxFunctionCoats_<<", cflFluxDefault = "<<cflFluxDefault<<"\n";
        if (std::isnan(cflFluxFunctionCoatsOut_) || std::isinf(cflFluxFunctionCoatsOut_)){cflFluxFunctionCoatsOut_ = 0.0;}
        if (std::isnan(cflFluxFunctionCoatsIn_) || std::isinf(cflFluxFunctionCoatsIn_)){cflFluxFunctionCoatsIn_ = 0.0;}

        Scalar cflFluxFunctionCoats = std::max(cflFluxFunctionCoatsIn_, cflFluxFunctionCoatsOut_);
//                std::cout<<"cflFluxFunctionCoatsIn = "<<cflFluxFunctionCoatsIn_<<"cflFluxFunctionCoatsOut = "<<cflFluxFunctionCoatsOut_<<", cflFluxDefault = "<<cflFluxDefault<<"\n";
        if (cflFluxFunctionCoats <= 0)
146
147
148
        {
            return 0.99 / cflFluxDefault;
        }
149
        else if (cflFluxDefault > cflFluxFunctionCoats)
150
        {
151
        	return 0.99 / cflFluxDefault;
152
153
154
        }
        else
        {
155
            return 0.99 / cflFluxFunctionCoats;
156
157
        }
    }
158

159
    /*! \brief  Returns the Cfl time-step
160
161
162
163
164
     *
     * \copydetails EvalCflFlux::getDt(const Element&)
     */
    Scalar getDt(const Element& element)
    {
165
166
//        if (rejectForTimeStepping_)
//            return 1e100;
167
168

        Scalar porosity = std::max(problem_.spatialParams().porosity(element), porosityThreshold_);
169
        //        if (porosity > porosityThreshold_)
170
        return (getCflFluxFunction(element) * porosity * element.geometry().volume());
171
172
        //        else
        //            return 1e100;//(getCflFluxFunction(element) * element.geometry().volume());
173
    }
174

175
176
177
    //! \brief  Resets the Timestep-estimator
    void reset()
    {
178
179
        cflFluxFunctionCoatsIn_ = 0;
        cflFluxFunctionCoatsOut_ = 0;
180
        rejectForTimeStepping_ = false;
181
182
183
184
185
        fluxWettingOut_ = 0;
        fluxNonwettingOut_ = 0;
        fluxIn_ = 0;
        fluxOut_ = 0;
    }
186

187
188
189
190
191
    /*! \brief Constructs an EvalCflFluxDefault object
     *
     * \param problem A problem type object
     */
    EvalCflFluxCoats(Problem& problem) :
192
        problem_(problem), epsDerivative_(5e-3), threshold_(1e-8)
193
194
195
196
    {
        reset();
        density_[wPhaseIdx] = 0.;
        density_[nPhaseIdx] = 0.;
197
198

        porosityThreshold_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, PorosityThreshold);
199
    }
200

201
private:
202
    Scalar getCflFluxFunctionDefault()
203
204
    {
        if (std::isnan(fluxIn_) || std::isinf(fluxIn_))
205
        {
206
207
            fluxIn_ = 1e-100;
        }
208

209
        Scalar cFLFluxIn = fluxIn_;
210
211


212
        Scalar cFLFluxOut = 0;
213

214
215
216
217
218
219
        if (velocityType_ == vt)
        {
            if (std::isnan(fluxOut_) || std::isinf(fluxOut_))
            {
                fluxOut_ = 1e-100;
            }
220

221
222
223
224
225
            cFLFluxOut = fluxOut_;
        }
        else
        {
            if (std::isnan(fluxWettingOut_) || std::isinf(fluxWettingOut_))
226
            {
227
228
229
230
231
                fluxWettingOut_ = 1e-100;
            }
            if (std::isnan(fluxNonwettingOut_) || std::isinf(fluxNonwettingOut_))
            {
                fluxNonwettingOut_ = 1e-100;
232
            }
233

234
235
            cFLFluxOut = std::max(fluxWettingOut_, fluxNonwettingOut_);
        }
236
237


238
239
        //determine timestep
        Scalar cFLFluxFunction = std::max(cFLFluxIn, cFLFluxOut);
240

241
242
243
244
        return cFLFluxFunction;
    }

    void addDefaultFlux(Scalar flux,int phaseIdx);
245

246
    void addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
247
                      const Intersection& intersection, int phaseIdx);
248

249
    Problem& problem_;//problem data
250
251
    Scalar cflFluxFunctionCoatsIn_;
    Scalar cflFluxFunctionCoatsOut_;
252
253
254
255
    Scalar fluxWettingOut_;
    Scalar fluxNonwettingOut_;
    Scalar fluxOut_;
    Scalar fluxIn_;
256
    bool rejectForTimeStepping_;
257
258
259
260
261
262
    Scalar density_[numPhases];
    static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation);
    static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation);
    static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation);
    const Scalar epsDerivative_;
    const Scalar threshold_;
263
    Scalar porosityThreshold_;
264
};
265

266
267
268
269
270
271
272
273
274
template<class TypeTag>
void EvalCflFluxCoats<TypeTag>::addDefaultFlux(Scalar flux, int phaseIdx)
{
    switch (phaseIdx)
    {
    case wPhaseIdx:
        {
            //for time step criterion
            if (flux >= 0)
275
            {
276
                fluxWettingOut_ += flux;
277
            }
278
            if (flux < 0)
279
            {
280
                fluxIn_ -= flux;
281
282
            }

283
284
            break;
        }
285

286
287
288
289
        //for time step criterion if the non-wetting phase velocity is used
    case nPhaseIdx:
        {
            if (flux >= 0)
290
            {
291
292
293
294
295
                fluxNonwettingOut_ += flux;
            }
            if (flux < 0)
            {
                fluxIn_ -= flux;
296
297
            }

298
299
300
301
302
            break;
        }
    default:
        {
            if (flux >= 0)
303
            {
304
                fluxOut_ += flux;
305
            }
306
            if (flux < 0)
307
            {
308
                fluxIn_ -= flux;
309
310
            }

311
312
313
314
            break;
        }
    }
}
315

316
317
318
319
320
/*! \brief adds a flux to the cfl-criterion evaluation
 *
 * \copydetails EvalCflFlux::addFlux(Scalar&,Scalar&,Scalar&,Scalar&,Scalar,const Intersection&,int)
 */
template<class TypeTag>
321
void EvalCflFluxCoats<TypeTag>::addCoatsFlux(Scalar& lambdaW, Scalar& lambdaNw, Scalar& viscosityW, Scalar& viscosityNw, Scalar flux,
322
323
                                             const Intersection& intersection, int phaseIdx)
{
324
325
326
327
    if (rejectForTimeStepping_)
        return;
    else if (phaseIdx != wPhaseIdx)
        return;
328

329
330
331
332
    Scalar lambdaT = (lambdaW + lambdaNw);

    if (lambdaT <= threshold_)
        return;
333

334
    ElementPointer element = intersection.inside();
335

336
337
    //coordinates of cell center
    const GlobalPosition& globalPos = element->geometry().center();
338

339
340
341
342
343
    // cell index
    int globalIdxI = problem_.variables().index(*element);

    CellData& cellDataI = problem_.variables().cellData(globalIdxI);

344
    if (cellDataI.potential(wPhaseIdx) < 0.0 || cellDataI.potential(nPhaseIdx) < 0.0)
345
    {
346
347
348
349
        rejectForTimeStepping_ = true;
        cflFluxFunctionCoatsIn_ = 0;
        cflFluxFunctionCoatsOut_ = 0;
        return;
350
351
    }

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

    Scalar satI = cellDataI.saturation(wPhaseIdx);
    Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
356
    Scalar lambdaNwI = cellDataI.mobility(nPhaseIdx);
357

358
    Scalar dpc_dsI = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(*element), satI);
359
360
361
362
363
364
365
366
367
368
369

    const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();

    if (intersection.neighbor())
    {
        ElementPointer neighbor = intersection.outside();

        int globalIdxJ = problem_.variables().index(*neighbor);

        CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);

370
        if (cellDataJ.potential(wPhaseIdx) < 0.0 || cellDataJ.potential(nPhaseIdx) < 0.0 )
371
        {
372
373
374
375
            rejectForTimeStepping_ = true;
            cflFluxFunctionCoatsIn_ = 0;
            cflFluxFunctionCoatsOut_ = 0;
            return;
376
        }
377
        
378
379
        if (element.level() != neighbor.level())
        {
380
            rejectForTimeStepping_ = true;
381
382
            cflFluxFunctionCoatsIn_ = 0;
            cflFluxFunctionCoatsOut_ = 0;
383
384
385
            return;
        }

386
387
388
389
390
391
392
393
        bool takeNeighbor = (element->level() < neighbor->level());
        //get phase potentials
        bool upwindWI =
            (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(wPhaseIdx, intersection.indexInOutside()) :
            cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside);
        bool upwindNwI =
            (takeNeighbor) ? !cellDataJ.fluxData().isUpwindCell(nPhaseIdx, intersection.indexInOutside()) :
            cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside);
394

395
            const GlobalPosition& globalPosNeighbor = neighbor->geometry().center();
396

397
398
            // distance vector between barycenters
            GlobalPosition distVec = globalPosNeighbor - globalPos;
399

400
401
            // compute distance between cell centers
            Scalar dist = std::abs(distVec * unitOuterNormal);
402

403
404
405
            Scalar satJ = cellDataJ.saturation(wPhaseIdx);
            Scalar lambdaWJ = cellDataI.mobility(wPhaseIdx);
            Scalar lambdaNwJ = cellDataI.mobility(nPhaseIdx);
406

407
            Scalar dpc_dsJ = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(*neighbor), satJ);
408

409
410
411
412
413
            // compute vectorized permeabilities
            DimVector permeability(0);
            DimMatrix perm(0);
            problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(*element));
            perm.mv(unitOuterNormal, permeability);
414

415
            Scalar perm1 = permeability * unitOuterNormal;
416

417
418
419
            permeability = 0;
            problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(*neighbor));
            perm.mv(unitOuterNormal, permeability);
420

421
            Scalar perm2 = permeability * unitOuterNormal;
422

423
            Scalar meanPermeability = 2*perm1*perm2/(perm1 + perm2);
424

425
            Scalar transmissibility =  meanPermeability * intersection.geometry().volume() / dist;
426

427
428
429
430
431
432
433
434
435
            Scalar satUpw = 0;
            if (upwindWI)
            {
                satUpw = std::max(satI, 0.0);
            }
            else
            {
                satUpw = std::max(satJ, 0.0);
            }
436

437
            Scalar ds = epsDerivative_;
438

439
440
441
442
443
444
445
            Scalar satPlus = satUpw + epsDerivative_;
            Scalar satMinus = satUpw;
            if (satMinus - epsDerivative_ > 0.0)
            {
                satMinus -= epsDerivative_;
                ds += epsDerivative_;
            }
446

447
448
449
            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satPlus)) / viscosityW;
            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*neighbor), std::abs(satMinus)) / viscosityW;
            dLambdaWDs /= (ds);
450

451
452
453
454
455
456
457
458
            if (upwindNwI)
            {
                satUpw = std::max(1 - satI, 0.0);
            }
            else
            {
                satUpw = std::max(1 - satJ, 0.0);
            }
459

460
461
462
463
464
            ds = epsDerivative_;

            satPlus = satUpw + epsDerivative_;
            satMinus = satUpw;
            if (satMinus - epsDerivative_ > 0.0)
465
            {
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
                satMinus -= epsDerivative_;
                ds += epsDerivative_;
            }

            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satPlus) / viscosityNw;
            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighbor), satMinus) / viscosityNw;
            dLambdaNwDs /= (ds);

            Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWJ);
            Scalar lambdaNwCap = 0.5 * (lambdaNwI + lambdaNwJ);

            Scalar potentialDiff = cellDataI.potential(wPhaseIdx) - cellDataJ.potential(wPhaseIdx);
            Scalar cflFlux = transmissibility * lambdaNw * dLambdaWDs * std::abs(potentialDiff) / lambdaT;

            potentialDiff = cellDataI.potential(nPhaseIdx) - cellDataJ.potential(nPhaseIdx);
            cflFlux -= transmissibility * lambdaW * dLambdaNwDs * std::abs(potentialDiff) / lambdaT;

            cflFlux -= transmissibility * lambdaWCap * lambdaNwCap * (dpc_dsI + dpc_dsJ) / lambdaT;

            if ((upwindWI && lambdaW > threshold_)|| (upwindNwI && lambdaW < threshold_))
            {
                cflFluxFunctionCoatsOut_ += cflFlux;
            }
            else
            {
            	cflFluxFunctionCoatsIn_ += cflFlux;
492
            }
493
494
495
    }
    else
    {
496
497
            // center of face in global coordinates
            GlobalPosition globalPosFace = intersection.geometry().center();
498

499
500
501
            //get boundary type
            BoundaryTypes bcType;
            problem_.boundaryTypes(bcType, intersection);
502

503
504
            // distance vector between barycenters
            Dune::FieldVector < Scalar, dimWorld > distVec = globalPosFace - globalPos;
505

506
507
            // compute distance between cell centers
            Scalar dist = distVec.two_norm();
508

509
510
            //permeability vector at boundary
            // compute vectorized permeabilities
511

512
513
514
515
            Dune::FieldVector<Scalar, dim> permeability(0);
            DimMatrix perm(0);
            problem_.spatialParams().meanK(perm, problem_.spatialParams().intrinsicPermeability(*element));
            perm.mv(unitOuterNormal, permeability);
516

517
            Scalar transmissibility = (unitOuterNormal * permeability) * intersection.geometry().volume() / dist;
518

519
520
            Scalar satWBound =  cellDataI.saturation(wPhaseIdx);
            if (bcType.isDirichlet(eqIdxSat))
521
            {
522
523
524
                PrimaryVariables bcValues;
                problem_.dirichlet(bcValues, intersection);
                switch (saturationType_)
525
                {
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
                case sw:
                    {
                        satWBound = bcValues[eqIdxSat];
                        break;
                    }
                case sn:
                    {
                        satWBound = 1 - bcValues[eqIdxSat];
                        break;
                    }
                default:
                    {
                        DUNE_THROW(Dune::RangeError, "saturation type not implemented");
                        break;
                    }
541
                }
542

543
            }
544

545
546
547
548
            Scalar potWBound =  cellDataI.potential(wPhaseIdx);
            Scalar potNwBound =  cellDataI.potential(nPhaseIdx);
            Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * problem_.gravity();
            if (bcType.isDirichlet(eqIdxPress))
549
            {
550
551
552
                PrimaryVariables bcValues;
                problem_.dirichlet(bcValues, intersection);
                switch (pressureType_)
553
                {
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
                case pw:
                    {
                        potWBound = bcValues[eqIdxPress] + density_[wPhaseIdx] * gdeltaZ;
                        potNwBound = bcValues[eqIdxPress] + MaterialLaw::pc(problem_.spatialParams().materialLawParams(*element), satWBound) + density_[nPhaseIdx] * gdeltaZ;
                        break;
                    }
                case pn:
                    {
                        potWBound = bcValues[eqIdxPress] - MaterialLaw::pc(problem_.spatialParams().materialLawParams(*element), satWBound) + density_[wPhaseIdx] * gdeltaZ;
                        potNwBound = bcValues[eqIdxPress] + density_[nPhaseIdx] * gdeltaZ;
                        break;
                    }
                default:
                    {
                        DUNE_THROW(Dune::RangeError, "pressure type not implemented");
                        break;
                    }
571
                }
572
573
574
575
576
577
578
            }
            else if (bcType.isNeumann(eqIdxPress))
            {
                PrimaryVariables bcValues;
                problem_.neumann(bcValues, intersection);

                if (lambdaW != 0 && bcValues[wPhaseIdx] != 0)
579
                {
580
                    potWBound += bcValues[wPhaseIdx] / (transmissibility * lambdaW);
581
                }
582
                if (lambdaNw != 0 && bcValues[nPhaseIdx] != 0)
583
                {
584
                    potNwBound += bcValues[nPhaseIdx] / (transmissibility * lambdaNw);
585
                }
586
587
            }

588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
            Scalar dpc_dsBound = MaterialLaw::dpc_dsw(problem_.spatialParams().materialLawParams(*element), satWBound);

            Scalar lambdaWBound = 0;
            Scalar lambdaNwBound = 0;

            Scalar temperature = problem_.temperature(*element);
            Scalar referencePressure = problem_.referencePressure(*element);
            FluidState fluidState;
            fluidState.setPressure(wPhaseIdx, referencePressure);
            fluidState.setPressure(nPhaseIdx, referencePressure);
            fluidState.setTemperature(temperature);

            Scalar viscosityWBound = FluidSystem::viscosity(fluidState, wPhaseIdx);
            Scalar viscosityNwBound =
                FluidSystem::viscosity(fluidState, nPhaseIdx);
            lambdaWBound = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityWBound;
            lambdaNwBound = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satWBound) / viscosityNwBound;

            Scalar satUpw = 0;
            if (cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside))
608
            {
609
                satUpw = std::max(satI, 0.0);
610
            }
611
            else
612
            {
613
                satUpw = std::max(satWBound, 0.0);
614
615
            }

616
            Scalar ds = epsDerivative_;
617

618
619
620
621
622
623
624
            Scalar satPlus = satUpw + epsDerivative_;
            Scalar satMinus = satUpw;
            if (satMinus - epsDerivative_ > 0.0)
            {
                satMinus -= epsDerivative_;
                ds += epsDerivative_;
            }
625

626
627
628
            Scalar dLambdaWDs = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityW;
            dLambdaWDs -= MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityW;
            dLambdaWDs /= (ds);
629

630
631
632
633
634
635
636
637
            if (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside))
            {
                satUpw = std::max(1 - satI, 0.0);
            }
            else
            {
                satUpw = std::max(1 - satWBound, 0.0);
            }
638

639
            ds = epsDerivative_;
640

641
642
643
644
645
646
647
            satPlus = satUpw + epsDerivative_;
            satMinus = satUpw;
            if (satMinus - epsDerivative_ > 0.0)
            {
                satMinus -= epsDerivative_;
                ds += epsDerivative_;
            }
648

649
650
651
            Scalar dLambdaNwDs = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satPlus) / viscosityNw;
            dLambdaNwDs -= MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satMinus) / viscosityNw;
            dLambdaNwDs /= (ds);
Markus Wolff's avatar
Markus Wolff committed
652

653
654
            Scalar lambdaWCap = 0.5 * (lambdaWI + lambdaWBound);
            Scalar lambdaNwCap = 0.5 * (lambdaNwI + lambdaNwBound);
655

656
657
            Scalar potDiff = cellDataI.potential(wPhaseIdx) - potWBound;
            Scalar cflFlux = transmissibility * lambdaNw * dLambdaWDs * std::abs(potDiff) / lambdaT;
658

659
            cflFlux -= transmissibility * lambdaWCap * lambdaNwCap * (dpc_dsI + dpc_dsBound) / lambdaT;
660

661
662
            potDiff = cellDataI.potential(nPhaseIdx) - potNwBound;
            cflFlux -= transmissibility * lambdaW * dLambdaNwDs * std::abs(potDiff) / lambdaT;
663

664
            if ((cellDataI.fluxData().isUpwindCell(wPhaseIdx, indexInInside) && lambdaW > threshold_) || (cellDataI.fluxData().isUpwindCell(nPhaseIdx, indexInInside) && lambdaW < threshold_))
665
            {
666
667
668
669
670
                cflFluxFunctionCoatsOut_ += cflFlux;
            }
            else
            {
            	cflFluxFunctionCoatsIn_ += cflFlux;
671
672
673
            }
    }
}
674
675
676
677

}

#endif