fvpressuremultiphysics.hh 41.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   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_FVPRESSURE2P2C_MULTIPHYSICS_HH
#define DUMUX_FVPRESSURE2P2C_MULTIPHYSICS_HH

#include <dune/common/float_cmp.hh>

// dumux environment
#include <dumux/porousmediumflow/2p2c/sequential/fvpressure.hh>
#include <dumux/linear/vectorexchange.hh>

/**
 * @file
 * @brief  Finite Volume 2p2c Pressure Model with MultiPhysics
 */

namespace Dumux
{
//! The finite volume model for the solution of the compositional pressure equation
/*! \ingroup multiphysics
 *  Provides a Finite Volume implementation for the pressure equation of a gas-liquid
 *  system with two components. An IMPES-like method is used for the sequential
 *  solution of the problem.  Diffusion is neglected, capillarity can be regarded.
 *  Isothermal conditions and local thermodynamic
 *  equilibrium are assumed.  Gravity is included.
 *  \f[
         c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}}
         \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{alpha} \bf{v}_{\alpha}\right)
          = \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} q^{\kappa},
 *  \f]
 *  where \f$\bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \f$.
 *  \f$ c_{total} \f$ represents the total compressibility, for constant porosity this yields
 *  \f$ - \frac{\partial V_{total}}{\partial p_{\alpha}} \f$,
 *  \f$p_{\alpha} \f$ denotes the phase pressure, \f$ \bf{K} \f$ the absolute permeability,
 *  \f$ \lambda_{\alpha} \f$ the phase mobility,
 *  \f$ \rho_{\alpha} \f$ the phase density and \f$ \bf{g} \f$ the gravity constant and
 *  \f$ C^{\kappa} \f$ the total Component concentration.
 * See paper SPE 99619 or "Analysis of a Compositional Model for Fluid
 * Flow in Porous Media" by Chen, Qin and Ewing for derivation.
 *
 *  The partial derivatives of the actual fluid volume \f$ v_{total} \f$ are gained by using a secant method.
 *
 * The model domain is automatically divided
 * in a single-phase and a two-phase domain. The full 2p2c model is only evaluated within the
 * two-phase subdomain, whereas a single-phase transport model is computed in the rest of the
 * domain.
 *
 * \tparam TypeTag The Type Tag
 */
template<class TypeTag>
class FVPressure2P2CMultiPhysics : public FVPressure2P2C<TypeTag>
{
    typedef FVPressure2P2C<TypeTag> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;

    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, BoundaryTypes) BoundaryTypes;

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

    typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };
    enum
    {
        pw = Indices::pressureW
    };
    enum
    {
        wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
        wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
        contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
    };

    // typedefs to abbreviate several dune classes...
    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::Grid Grid;
    typedef typename GridView::Intersection Intersection;

    // convenience shortcuts for Vectors/Matrices
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef Dune::FieldMatrix<Scalar, dim, dim> DimMatrix;
    typedef Dune::FieldVector<Scalar, GET_PROP_VALUE(TypeTag, NumPhases)> PhaseVector;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;

    //! @copydoc FVPressure::EntryType
    typedef Dune::FieldVector<Scalar, 2> EntryType;
//! Access functions to the current problem object
    Problem& problem()
    {    return this->problem_;   }
    const Problem& problem() const
    {    return this->problem_;   }
public:
    //function which assembles the system of equations to be solved
    void assemble(bool first);

    void get1pSource(EntryType&, const Element&, const CellData&);

    void get1pStorage(EntryType&, const Element&, CellData&);

    void get1pFlux(EntryType&, const Intersection&, const CellData&);

    void get1pFluxOnBoundary(EntryType&,
                            const Intersection&, const CellData&);

    //initialize multi-physics-specific pressure model constituents
    void initialize(bool solveTwice = false)
    {
        // assign whole domain to most complex subdomain, => 2p
        int size = this->problem().gridView().size(0);
        for (int i = 0; i < size; i++)
        {
            CellData& cellData = this->problem().variables().cellData(i);
            cellData.subdomain() = 2;
        }
        nextSubdomain.resize(size);
        ParentType::initialize();
    }

    /*! \brief  Function for serialization of the pressure field.
     *
     *  Function needed for restart option. Writes the pressure of a grid element to a restart file.
     *
     *  \param outstream Stream into the restart file.
     *  \param element Grid element
     */
    void serializeEntity(std::ostream &outstream, const Element &element)
    {
        ParentType::serializeEntity(outstream,element);
        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
        outstream <<"  "<< cellData.subdomain();
    }

    /*! \brief  Function for deserialization of the pressure field.
     *
     *  Function needed for restart option. Reads the pressure of a grid element from a restart file.
     *
     *  \param instream Stream from the restart file.
     *  \param element Grid element
     */
    void deserializeEntity(std::istream &instream, const Element &element)
    {
        ParentType::deserializeEntity(instream,element);

        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
        int subdomainIdx;
        instream >> subdomainIdx;
        cellData.setSubdomainAndFluidStateType(subdomainIdx);
    }


    //constitutive functions are initialized and stored in the variables object
    void updateMaterialLaws(bool postTimeStep = false);
    //updates singlephase secondary variables for one cell and stores in the variables object
    void update1pMaterialLawsInElement(const Element& elementI, CellData& cellData, bool postTimeStep);

    //! \brief Write data files
     /*  \param name file name */
    template<class MultiWriter>
    void addOutputVtkFields(MultiWriter &writer)
    {
        FVPressureCompositional<TypeTag>::addOutputVtkFields(writer);

        if(problem().vtkOutputLevel()>=1)
        {
            int size = problem().gridView().size(0);
            // add multiphysics stuff
            Dune::BlockVector<Dune::FieldVector<int,1> >* subdomainPtr = writer.template allocateManagedBuffer<int, 1> (size);
            for (int i = 0; i < size; i++)
            {
                CellData& cellData = problem().variables().cellData(i);
                (*subdomainPtr)[i] = cellData.subdomain();
            }
            writer.attachCellData(*subdomainPtr, "subdomain");
        }

        return;
    }

    //! Constructs a FVPressure2P2CPC object
    /**
     * \param problem a problem class object
     */
    FVPressure2P2CMultiPhysics(Problem& problem) : FVPressure2P2C<TypeTag>(problem),
            gravity_(problem.gravity()), timer_(false)
    {}

protected:
    #if HAVE_MPI
        typedef typename SolutionTypes::ElementMapper ElementMapper;
        typedef VectorExchange<ElementMapper, Dune::BlockVector<Dune::FieldVector<int, 1> > > DataHandle;
    #endif

    // subdomain map
    Dune::BlockVector<Dune::FieldVector<int,1> > nextSubdomain;  //! vector holding next subdomain
    const GlobalPosition& gravity_; //!< vector including the gravity constant
    //! gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
    static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation);
    Dune::Timer timer_; //!< A timer for the time spent on the multiphysics framework.

    //! Indices of matrix and rhs entries
    /**
    * During the assembling of the global system of equations get-functions are called (getSource(),
    * getFlux(), etc.), which return global matrix or right hand side entries in a vector.
    * These can be accessed using following indices:
    */
    enum
    {
        rhs = 1,//!<index for the right hand side entry
        matrix = 0//!<index for the global matrix entry

    };
};

//! function which assembles the system of equations to be solved
/*! for first == true, this function assembles the matrix and right hand side for
 * the solution of the pressure field in the same way as in the class FVPressure2P.
 * for first == false, the approach is changed to \f[-\frac{\partial V}{\partial p}
 * \frac{\partial p}{\partial t}+\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}\nabla\cdot
 * \left(\sum_{\alpha}C_{\alpha}^{\kappa}\mathbf{v}_{\alpha}\right)
 * =\sum_{\kappa}\frac{\partial V}{\partial m^{\kappa}}q^{\kappa} \f]. See Paper SPE 99619.
 * This is done to account for the volume effects which appear when gas and liquid are dissolved in each other.
 * \param first Flag if pressure field is unknown
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
{
    if(first)
    {
        ParentType::assemble(true);
        return;
    }
    // initialization: set matrix this->A_ to zero
    this->A_ = 0;
    this->f_ = 0;

264
    for (const auto& element : elements(problem().gridView()))
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    {
        // get the global index of the cell
        int eIdxGlobalI = problem().variables().index(element);

        // assemble interior element contributions
        if (element.partitionType() == Dune::InteriorEntity)
        {
            // get the cell data
            CellData& cellDataI = problem().variables().cellData(eIdxGlobalI);

            Dune::FieldVector<Scalar, 2> entries(0.);

            /*****  source term ***********/
            if(cellDataI.subdomain() != 2)
                problem().pressureModel().get1pSource(entries,element, cellDataI);
            else
                problem().pressureModel().getSource(entries,element, cellDataI, first);

            this->f_[eIdxGlobalI] = entries[rhs];

            /*****  flux term ***********/
            // iterate over all faces of the cell
287
            for (const auto& intersection : intersections(problem().gridView(), element))
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
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
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
            {
                /************* handle interior face *****************/
                if (intersection.neighbor())
                {
                    int eIdxGlobalJ = problem().variables().index(intersection.outside());

                    if (cellDataI.subdomain() != 2
                            or problem().variables().cellData(eIdxGlobalJ).subdomain() != 2) // cell in the 1p domain
                        get1pFlux(entries, intersection, cellDataI);
                    else
                        problem().pressureModel().getFlux(entries, intersection, cellDataI, first);

                    //set right hand side
                    this->f_[eIdxGlobalI] -= entries[rhs];
                    // set diagonal entry
                    this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
                    // set off-diagonal entry
                    this->A_[eIdxGlobalI][eIdxGlobalJ] = -entries[matrix];
                }   // end neighbor


                /************* boundary face ************************/
                else
                {
                    if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
                        problem().pressureModel().get1pFluxOnBoundary(entries, intersection, cellDataI);
                    else
                        problem().pressureModel().getFluxOnBoundary(entries, intersection, cellDataI, first);

                    //set right hand side
                    this->f_[eIdxGlobalI] += entries[rhs];
                    // set diagonal entry
                    this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
                }
            } //end interfaces loop
    //        printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);

            /*****  storage term ***********/
            if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
                problem().pressureModel().get1pStorage(entries, element, cellDataI);
            else
                problem().pressureModel().getStorage(entries, element, cellDataI, first);

            this->f_[eIdxGlobalI] += entries[rhs];
            // set diagonal entry
            this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
        }
        // assemble overlap and ghost element contributions
        else
        {
            this->A_[eIdxGlobalI] = 0.0;
            this->A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
            this->f_[eIdxGlobalI] = this->pressure()[eIdxGlobalI];
        }
    } // end grid traversal
//        printmatrix(std::cout, this->A_, "global stiffness matrix after assempling", "row", 11,3);
//        printvector(std::cout, this->f_, "right hand side", "row", 10);
    return;
}

//! Assembles the source term
/** The source is translated into a volumentric source term:
 * \f[ V_i \sum_{\kappa} \frac{1}{\varrho} q^{\kappa}_i \; , \f]
 * because under singlephase conditions
 * \f[ \frac{\partial v_{t}}{\partial C^{\kappa}} \approx \frac{1}{\varrho} \f].
 * \param sourceEntry The Matrix and RHS entries
 * \param elementI The element I
 * \param cellDataI Data of cell I
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pSource(Dune::FieldVector<Scalar, 2>& sourceEntry,
        const Element& elementI, const CellData& cellDataI)
{
    sourceEntry=0.;

    // cell volume & perimeter, assume linear map here
    Scalar volume = elementI.geometry().volume();
    int subdomainIdx = cellDataI.subdomain();

    /****************implement source************************/
    PrimaryVariables source(NAN);
    problem().source(source, elementI);
    source[1+subdomainIdx] /= cellDataI.density(subdomainIdx);

    sourceEntry[1] = volume * source[1+subdomainIdx];

    return;
}

//! Assembles the storage term for a 1p cell in a multiphysics framework
/** The storage term comprises the (single-phase) compressibility (due to a change in
 * pressure from last timestep):
 *  \f[ V_i c_{i} \frac{p^t_i - p^{t-\Delta t}_i}{\Delta t} \f]
 * and the damped error introduced by the incorrect transport of the last timestep:
 *  \f[ V_i \alpha_r \frac{v_{t} - \phi}{\Delta t} \f].
 * The latter is damped according to Fritz 2011.
 * \param storageEntry The Matrix and RHS entries
 * \param elementI The element I
 * \param cellDataI Data of cell I
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar, 2>& storageEntry,
                                                        const Element& elementI,
                                                        CellData& cellDataI)
{
    storageEntry = 0.;
    // cell index
    int eIdxGlobalI = problem().variables().index(elementI);
    int presentPhaseIdx = cellDataI.subdomain();
    Scalar volume = elementI.geometry().volume();

    // determine maximum error to scale error-term
    Scalar timestep_ = problem().timeManager().timeStepSize();

    // compressibility term: 1p domain, so no dv_dp calculated
    if (true)
    {
        Scalar& incp = this->incp_;

        // numerical derivative of fluid volume with respect to pressure
        PhaseVector p_(incp);
        p_[nPhaseIdx] += cellDataI.pressure(nPhaseIdx);
        p_[wPhaseIdx] += cellDataI.pressure(wPhaseIdx);

        Scalar sumC = (cellDataI.massConcentration(wCompIdx) + cellDataI.massConcentration(nCompIdx));
        Scalar Z1 = cellDataI.massConcentration(wCompIdx) / sumC;
        // initialize simple fluidstate object
        PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
        CompositionalFlash<TypeTag> flashSolver;
        flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
                cellDataI.temperature(wPhaseIdx));

        Scalar v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
        cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;

        if (cellDataI.dv_dp()>0)
        {
            // dV_dp > 0 is unphysical: Try inverse increment for secant
            Dune::dinfo << "dv_dp larger 0 at Idx " << eIdxGlobalI << " , try and invert secant"<< std::endl;

            p_ -= 2*incp;
            flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
                    cellDataI.temperature(wPhaseIdx));
            v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
            cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;
            // dV_dp > 0 is unphysical: Try inverse increment for secant
            if (cellDataI.dv_dp()>0)
            {
                Dune::dinfo <<__FILE__<< "dv_dp still larger 0 after inverting secant. regularize"<< std::endl;
                cellDataI.dv_dp() *= -1;
            }
        }


        Scalar compress_term = cellDataI.dv_dp() / timestep_;

        storageEntry[matrix] -= compress_term*volume;
        storageEntry[rhs] -= cellDataI.pressure(pressureType) * compress_term * volume;

        if (std::isnan(compress_term) || std::isinf(compress_term))
            DUNE_THROW(Dune::MathError, "Compressibility term leads to NAN matrix entry at index " << eIdxGlobalI);

        if(!GET_PROP_VALUE(TypeTag, EnableCompressibility))
            DUNE_THROW(Dune::NotImplemented, "Compressibility is switched off???");
    }

    // Abort error damping if there will be a possibly tiny timestep compared with last one
    // This might be the case if the episode or simulation comes to an end.
456
    if( problem().timeManager().episodeWillBeFinished()
457
458
459
460
461
462
463
464
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
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
            || problem().timeManager().willBeFinished())
    {
        problem().variables().cellData(eIdxGlobalI).errorCorrection() = 0.;
        return;
    }

    // error reduction routine: volumetric error is damped and inserted to right hand side
    // if damping is not done, the solution method gets unstable!
    problem().variables().cellData(eIdxGlobalI).volumeError() /= timestep_;
    Scalar maxError = this->maxError_;
    Scalar erri = fabs(cellDataI.volumeError());
    Scalar x_lo = this->ErrorTermLowerBound_;
    Scalar x_mi = this->ErrorTermUpperBound_;
    Scalar fac  = this->ErrorTermFactor_;
    if (pressureType == pw)
        fac = 0.1*this->ErrorTermFactor_;
    Scalar lofac = 0.;
    Scalar hifac = 0.;

    if ((erri*timestep_ > 5e-5) && (erri > x_lo * maxError) && (!problem().timeManager().willBeFinished()))
    {
        if (erri <= x_mi * maxError)
            storageEntry[rhs] +=
                    problem().variables().cellData(eIdxGlobalI).errorCorrection() =
                            fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxError)
                                * cellDataI.volumeError() * volume;
        else
            storageEntry[rhs] +=
                    problem().variables().cellData(eIdxGlobalI).errorCorrection() =
                            fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError)
                                * cellDataI.volumeError() * volume;
    }
    else
        problem().variables().cellData(eIdxGlobalI).errorCorrection()=0 ;

    return;
}

/*! \brief The compositional single-phase flux in the multiphysics framework
 *
 * If only single-phase conditions are encountered, the flux expression simplifies to (written for the
 * case where the wetting phase is only present):
   \f[
          A_{\gamma} \mathbf{n}_{\gamma}^T \mathbf{K}
       \lambda_w \mathbf{d}_{ij} \left( \frac{p_{w,j}^t - p^{t}_{w,i}}{\Delta x} + \varrho_{w} \mathbf{g}^T \mathbf{d}_{ij} \right) .
   \f]
 *
 * \param entries The Matrix and RHS entries
 * \param intersection Intersection between cell I and J
 * \param cellDataI Data of cell I
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>& entries,
        const Intersection& intersection, const CellData& cellDataI)
{
    entries = 0.;
    auto elementI = intersection.inside();

    // get global coordinate of cell center
    const GlobalPosition& globalPos = elementI.geometry().center();

    // get absolute permeability
    DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(elementI));

    // get normal vector
    const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();

    // get face volume
    Scalar faceArea = intersection.geometry().volume();

        // access neighbor
        auto neighbor = intersection.outside();
        int eIdxGlobalJ = problem().variables().index(neighbor);
        CellData& cellDataJ = problem().variables().cellData(eIdxGlobalJ);

        // gemotry info of neighbor
        const GlobalPosition& globalPosNeighbor = neighbor.geometry().center();

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

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

        GlobalPosition unitDistVec(distVec);
        unitDistVec /= dist;

        DimMatrix permeabilityJ
            = problem().spatialParams().intrinsicPermeability(neighbor);

        // compute vectorized permeabilities
        DimMatrix meanPermeability(0);
549
        harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
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
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782

        Dune::FieldVector<Scalar, dim> permeability(0);
        meanPermeability.mv(unitDistVec, permeability);

        // due to "safety cell" around subdomain, both cells I and J
        // have single-phase conditions, although one is in 2p domain.
        int phaseIdx = std::min(cellDataI.subdomain(), cellDataJ.subdomain());

        Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + cellDataJ.density(phaseIdx));

        // 1p => no pc => only 1 pressure, potential
        Scalar potential = (cellDataI.pressure(phaseIdx) - cellDataJ.pressure(phaseIdx)) / dist;

        potential += rhoMean * (unitDistVec * gravity_);

        Scalar lambda;

        if (potential > 0.)
        {
            lambda = cellDataI.mobility(phaseIdx);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);  // store in cellJ since cellI is const
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);  // store in cellJ since cellI is const
        }
        else if (potential < 0.)
        {
            lambda = cellDataJ.mobility(phaseIdx);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, true);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, true);
        }
        else
        {
            lambda = harmonicMean(cellDataI.mobility(phaseIdx) , cellDataJ.mobility(phaseIdx));
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
        }

        entries[0]  = lambda * faceArea * fabs(permeability * unitOuterNormal) / (dist);
        entries[1]  = rhoMean * lambda;
        entries[1] *= faceArea * fabs(permeability * unitOuterNormal) * (unitDistVec * gravity_);

        return;
}

/*! \brief The compositional single-phase flux in the multiphysics framework
 *
 * If only single-phase conditions are encountered, the flux expression simplifies to (written for the
 * case where the wetting phase is only present):
   \f[
          A_{\gamma} \mathbf{n}_{\gamma}^T \mathbf{K}
      \varrho_w \lambda_w \mathbf{d}_{i-Boundary} \left( \frac{p_{w,Boundary}^t - p^{t}_{w,i}}{\Delta x}
      + \varrho_{w} \mathbf{g}^T \mathbf{d}_{i-Boundary} \right) .
   \f]
 *
 * If a Neumann BC is set, the given (mass-)flux is directly multiplied by the volume derivative and inserted.
 *
 * \param entries The Matrix and RHS entries
 * \param intersection Intersection between cell I and J
 * \param cellDataI Data of cell I
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entries,
        const Intersection& intersection, const CellData& cellDataI)
{
    entries = 0.;
    // get global coordinate of cell center
    auto elementI = intersection.inside();
    const GlobalPosition& globalPos = elementI.geometry().center();
//    int eIdxGlobalI = problem().variables().index(elementI);
    int phaseIdx = cellDataI.subdomain();

    // get normal vector
    const GlobalPosition& unitOuterNormal = intersection.centerUnitOuterNormal();
    // get face volume
    Scalar faceArea = intersection.geometry().volume();

                // center of face in global coordinates
                const GlobalPosition& globalPosFace = intersection.geometry().center();

                // geometrical information
                GlobalPosition distVec(globalPosFace - globalPos);
                Scalar dist = distVec.two_norm();
                GlobalPosition unitDistVec(distVec);
                unitDistVec /= dist;

                //get boundary condition for boundary face center
                BoundaryTypes bcType;
                problem().boundaryTypes(bcType, intersection);

                // prepare pressure boundary condition
                PhaseVector pressBC(0.);

                /**********         Dirichlet Boundary        *************/
                if (bcType.isDirichlet(Indices::pressureEqIdx))
                {
                    // get absolute permeability
                    DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(elementI));
                    if(this->regulateBoundaryPermeability)
                    {
                        int axis = intersection.indexInInside() / 2;
                        if(permeabilityI[axis][axis] < this->minimalBoundaryPermeability)
                            permeabilityI[axis][axis] = this->minimalBoundaryPermeability;
                    }
                    // get mobilities and fractional flow factors
                    Scalar lambdaI = cellDataI.mobility(phaseIdx);

                    //permeability vector at boundary
                    Dune::FieldVector<Scalar, dim> permeability(0);
                    permeabilityI.mv(unitDistVec, permeability);

                    // create a fluid state for the boundary
                    FluidState BCfluidState;

                    //read boundary values
                    PrimaryVariables primaryVariablesOnBoundary(NAN);
                    problem().dirichlet(primaryVariablesOnBoundary, intersection);

                    {
                        // read boundary values
                        problem().transportModel().evalBoundary(globalPosFace,
                                                                    intersection,
                                                                    BCfluidState,
                                                                    pressBC);

                        // determine fluid properties at the boundary
                        Scalar densityBound =
                            FluidSystem::density(BCfluidState, phaseIdx);
                        Scalar viscosityBound =
                            FluidSystem::viscosity(BCfluidState, phaseIdx);

                        // mobility at the boundary
                        Scalar lambdaBound = 0.;
                        switch (GET_PROP_VALUE(TypeTag, BoundaryMobility))
                        {
                        case Indices::satDependent:
                            {
                                lambdaBound = BCfluidState.saturation(phaseIdx)
                                    / viscosityBound;
                            break;
                            }
                        case Indices::permDependent:
                            {
                            if (phaseIdx == wPhaseIdx)
                                lambdaBound = MaterialLaw::krw(
                                    problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
                                    / viscosityBound;
                            else
                                lambdaBound = MaterialLaw::krn(
                                    problem().spatialParams().materialLawParams(elementI), BCfluidState.saturation(wPhaseIdx))
                                    / viscosityBound;
                            break;
                            }
                        }
                        Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + densityBound);

                        Scalar potential = 0;

                        //calculate potential gradient: pc = 0;
                        potential = (cellDataI.pressure(phaseIdx) - pressBC[phaseIdx]) / dist;

                        potential += rhoMean * (unitDistVec * gravity_);

                        Scalar lambda(0.);

                        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potential, 0.0, 1.0e-30))
                        {
                            lambda = 0.5*(lambdaI + lambdaBound);
                        }
                        else if (potential > 0.)
                        {
                            lambda = lambdaI;
                        }
                        else
                        {
                            lambda = lambdaBound;
                        }

                        //calculate current matrix entry
                        Scalar entry(0.), rightEntry(0.);
                        entry = lambda * (fabs(permeability * unitOuterNormal) / dist) * faceArea;

                        //calculate right hand side
                        rightEntry = lambda * rhoMean * fabs(permeability * unitOuterNormal)
                                * faceArea ;

                        // set diagonal entry and right hand side entry
                        entries[0] += entry;
                        entries[1] += entry * primaryVariablesOnBoundary[Indices::pressureEqIdx];
                        entries[1] -= rightEntry * (gravity_ * unitDistVec);
                    }    //end of if(first) ... else{...
                }   // end dirichlet

                /**********************************
                 * set neumann boundary condition
                 **********************************/
                else if(bcType.isNeumann(Indices::pressureEqIdx))
                {
                    PrimaryVariables J(NAN);
                    problem().neumann(J, intersection);
                    J[1+phaseIdx] /= cellDataI.density(phaseIdx);

                    entries[1] -= J[1+phaseIdx] * faceArea;
                }
                else
                    DUNE_THROW(Dune::NotImplemented, "Boundary Condition neither Dirichlet nor Neumann!");


    return;
}


//! constitutive functions are updated once if new concentrations are calculated and stored in the variables container
/*!
 * In contrast to the standard sequential 2p2c model ( FVPressure2P2C<TypeTag>::updateMaterialLaws() ),
 * this method also holds routines to adapt the subdomain. The subdomain indicates weather we are in 1p domain (value = 1)
 * or in the two phase subdomain (value = 2).
 * Note that the type of flash, i.e. the type of FluidState (FS), present in each cell does not have to
 * coincide with the subdomain. If a cell will be simple and was complex, a complex FS is available, so next time step
 * will use this complex FS, but updateMaterialLaw afterwards will finally transform that to simple FS.
 *  \param postTimeStep Flag indicating method is called from Problem::postTimeStep()
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
{
    //get timestep for error term
    Scalar maxError = 0.;

    // next subdomain map
    if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().time(), 0.0, 1.0e-30))
        nextSubdomain = 2;  // start with complicated sub in initialization
    else
        nextSubdomain = -1;  // reduce complexity after first TS

    // Loop A) through leaf grid
783
    for (const auto& element : elements(problem().gridView()))
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
    {
        // get global coordinate of cell center
        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);

        if(cellData.subdomain() == 2)    // complex
        {
            this->updateMaterialLawsInElement(element, postTimeStep);

            // check subdomain consistency
            timer_.start();
            // enshure we are not at source
            // get sources from problem
            PrimaryVariables source(NAN);
            problem().source(source, element);

            if ((cellData.saturation(wPhaseIdx) > 0.0 && cellData.saturation(wPhaseIdx) < 1.0)
                || Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(source.one_norm(), 0.0, 1.0e-30)) // cell still 2p
            {
                // mark this element
                nextSubdomain[eIdxGlobal] = 2;

                // mark neighbors
807
                for (const auto& intersection : intersections(problem().gridView(), element))
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
834
835
836
837
838
839
840
841
842
843
                {
                    if (intersection.neighbor())
                    {
                        int eIdxGlobalJ = problem().variables().index(intersection.outside());
                        // mark neighbor Element
                        nextSubdomain[eIdxGlobalJ] = 2;
                    }
                }
            }
            else if(nextSubdomain[eIdxGlobal] != 2)// update next subdomain if possible
            {
                if(Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(wPhaseIdx), 0.0, 1.0e-30))
                    nextSubdomain[eIdxGlobal] = wPhaseIdx;
                else if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(nPhaseIdx), 0.0, 1.0e-30))
                    nextSubdomain[eIdxGlobal] = nPhaseIdx;
            }
            timer_.stop();
            // end subdomain check
        }// end complex domain
        else if (nextSubdomain[eIdxGlobal] != 2) //check if cell remains in simple subdomain
            nextSubdomain[eIdxGlobal] = cellData.subdomain();

    } //end define complex area of next subdomain

    timer_.start();
    //communicate next subdomain if parallel
    #if HAVE_MPI
    // communicate updated values
    DataHandle dataHandle(problem().variables().elementMapper(), nextSubdomain);
    problem().gridView().template communicate<DataHandle>(dataHandle,
                                                        Dune::InteriorBorder_All_Interface,
                                                        Dune::ForwardCommunication);
    #endif

    // Loop B) thorugh leaf grid
    // investigate cells that were "simple" in current TS
844
    for (const auto& element : elements(problem().gridView()))
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
    {
        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);

        // store old subdomain information and assign new info
        int oldSubdomainI = cellData.subdomain();
        cellData.subdomain() = nextSubdomain[eIdxGlobal];

        //first check if simple will become complicated
        if(oldSubdomainI != 2
                    && nextSubdomain[eIdxGlobal] == 2)
        {
            // use complex update of the fluidstate
            timer_.stop();
            this->updateMaterialLawsInElement(element, postTimeStep);
            timer_.start();
        }
        else if(oldSubdomainI != 2
                    && nextSubdomain[eIdxGlobal] != 2)    // will be simple and was simple
        {
            // perform simple update
            this->update1pMaterialLawsInElement(element, cellData, postTimeStep);
        }
        //else
        // a) will remain complex -> everything already done in loop A
        // b) will be simple and was complex: complex FS available, so next TS
        //         will use comlex FS, next updateMaterialLaw will transform to simple

        maxError = std::max(maxError, fabs(cellData.volumeError()));
    }// end grid traversal
    this->maxError_ = maxError/problem().timeManager().timeStepSize();

    timer_.stop();

879
    if(problem().timeManager().willBeFinished() or problem().timeManager().episodeWillBeFinished())
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
        Dune::dinfo << "Subdomain routines took " << timer_.elapsed() << " seconds" << std::endl;

    return;
}

//! updates secondary variables of one single phase cell
/*! For each element, the secondary variables are updated according to the
 * primary variables. Only a simple flash calulation has to be carried out,
 * as phase distribution is already known: single-phase.
 * \param elementI The element
 * \param cellData The cell data of the current element
 * \param postTimeStep Flag indicating if we have just completed a time step
 */
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::update1pMaterialLawsInElement(const Element& elementI, CellData& cellData, bool postTimeStep)
{
    // get global coordinate of cell center
    GlobalPosition globalPos = elementI.geometry().center();
    int eIdxGlobal = problem().variables().index(elementI);

    // determine which phase should be present
    int presentPhaseIdx = cellData.subdomain(); // this is already =nextSubomainIdx

    // reset to calculate new timeStep if we are at the end of a time step
    if(postTimeStep)
        cellData.reset();

    // acess the simple fluid state and prepare for manipulation
    PseudoOnePTwoCFluidState<TypeTag>& pseudoFluidState
        = cellData.manipulateSimpleFluidState();

    // prepare phase pressure for fluid state
    // both phase pressures are necessary for the case 1p domain is assigned for
    // the next 2p subdomain
    PhaseVector pressure(0.);
    Scalar pc = 0;
    if(GET_PROP_VALUE(TypeTag, EnableCapillarity))
        pc = MaterialLaw::pc(problem().spatialParams().materialLawParams(elementI),
            ((presentPhaseIdx == wPhaseIdx) ? 1. : 0.)); // assign sw = 1 if wPhase present, else 0
    if(pressureType == wPhaseIdx)
    {
        pressure[wPhaseIdx] = this->pressure(eIdxGlobal);
        pressure[nPhaseIdx] = this->pressure(eIdxGlobal)+pc;
    }
    else
    {
        pressure[wPhaseIdx] = this->pressure(eIdxGlobal)-pc;
        pressure[nPhaseIdx] = this->pressure(eIdxGlobal);
    }

    // get the overall mass of component 1:  Z1 = C^k / (C^1+C^2) [-]
    Scalar sumConc = cellData.massConcentration(wCompIdx)
            + cellData.massConcentration(nCompIdx);
    Scalar Z1 = cellData.massConcentration(wCompIdx)/ sumConc;

    CompositionalFlash<TypeTag> flashSolver;
    flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, pressure, presentPhaseIdx, problem().temperatureAtPos(globalPos));

    // write stuff in fluidstate
    assert(presentPhaseIdx == pseudoFluidState.presentPhaseIdx());

//    cellData.setSimpleFluidState(pseudoFluidState);

    // initialize viscosities
    cellData.setViscosity(presentPhaseIdx, FluidSystem::viscosity(pseudoFluidState, presentPhaseIdx));

    // initialize mobilities
    if(presentPhaseIdx == wPhaseIdx)
    {
        cellData.setMobility(wPhaseIdx,
            MaterialLaw::krw(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
                / cellData.viscosity(wPhaseIdx));
        cellData.setMobility(nPhaseIdx, 0.);
    }
    else
    {
        cellData.setMobility(nPhaseIdx,
            MaterialLaw::krn(problem().spatialParams().materialLawParams(elementI), pseudoFluidState.saturation(wPhaseIdx))
                / cellData.viscosity(nPhaseIdx));
        cellData.setMobility(wPhaseIdx, 0.);
    }

    // error term handling
    Scalar vol(0.);
    vol = sumConc / pseudoFluidState.density(presentPhaseIdx);

    if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(problem().timeManager().timeStepSize(), 0.0, 1.0e-30))
        cellData.volumeError() = (vol - problem().spatialParams().porosity(elementI));
    return;
}

}//end namespace Dumux
#endif