fvpressure2p2cmultiphysics.hh 42.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
#ifndef DUMUX_FVPRESSURE2P2C_MULTIPHYSICS_HH
#define DUMUX_FVPRESSURE2P2C_MULTIPHYSICS_HH
21

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

24
// dumux environment
25
#include <dumux/decoupled/2p2c/fvpressure2p2c.hh>
26
#include <dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh>
27
#include <dumux/linear/vectorexchange.hh>
28
29
30

/**
 * @file
31
 * @brief  Finite Volume 2p2c Pressure Model with MultiPhysics
32
33
34
35
36
37
38
39
 */

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
40
41
 *  solution of the problem.  Diffusion is neglected, capillarity can be regarded.
 *  Isothermal conditions and local thermodynamic
42
43
 *  equilibrium are assumed.  Gravity is included.
 *  \f[
44
45
         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)
46
47
48
          = \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$.
49
50
51
52
53
54
 *  \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.
55
56
57
58
 * 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.
59
 *
60
61
62
63
 * 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.
64
 *
65
66
 * \tparam TypeTag The Type Tag
 */
67
68
template<class TypeTag>
class FVPressure2P2CMultiPhysics : public FVPressure2P2C<TypeTag>
69
{
70
    typedef FVPressure2P2C<TypeTag> ParentType;
71
72
73
74
    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;
75

76
77
    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
    typedef typename SpatialParams::MaterialLaw MaterialLaw;
78

79
80
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
81

82
83
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
84

85
    typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
86
87
88
89
90
91
    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };
    enum
    {
92
        pw = Indices::pressureW
93
94
95
96
    };
    enum
    {
        wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
97
98
        wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
        contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
99
100
101
102
103
104
105
    };

    // typedefs to abbreviate several dune classes...
    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::Grid Grid;
    typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
106
    typedef typename GridView::Intersection Intersection;
107
108
109
110
    typedef typename GridView::IntersectionIterator IntersectionIterator;

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

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

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

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

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

132
    void get1pFluxOnBoundary(EntryType&,
133
134
                            const Intersection&, const CellData&);

135
136
    //initialize multi-physics-specific pressure model constituents
    void initialize(bool solveTwice = false)
137
138
139
140
141
142
143
144
145
146
147
148
    {
        // 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();
    }

149
150
151
152
153
154
155
156
157
158
    /*! \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);
159
160
        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
161
162
163
164
165
166
167
168
169
170
171
172
173
174
        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);

175
176
        int eIdxGlobal = problem().variables().index(element);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
177
178
179
180
181
        int subdomainIdx;
        instream >> subdomainIdx;
        cellData.setSubdomainAndFluidStateType(subdomainIdx);
    }

182

183
    //constitutive functions are initialized and stored in the variables object
Christoph Grueninger's avatar
Christoph Grueninger committed
184
    void updateMaterialLaws(bool postTimeStep = false);
185
    //updates singlephase secondary variables for one cell and stores in the variables object
Christoph Grueninger's avatar
Christoph Grueninger committed
186
    void update1pMaterialLawsInElement(const Element& elementI, CellData& cellData, bool postTimeStep);
187
188
189
190
191
192

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

195
        if(problem().vtkOutputLevel()>=1)
196
        {
197
198
199
200
201
202
203
204
205
            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");
206
        }
207
208
209
210

        return;
    }

211
    //! Constructs a FVPressure2P2CPC object
212
213
214
    /**
     * \param problem a problem class object
     */
215
    FVPressure2P2CMultiPhysics(Problem& problem) : FVPressure2P2C<TypeTag>(problem),
216
            gravity_(problem.gravity()), timer_(false)
217
    {}
218

219
protected:
220
221
222
223
224
    #if HAVE_MPI
        typedef typename SolutionTypes::ElementMapper ElementMapper;
        typedef VectorExchange<ElementMapper, Dune::BlockVector<Dune::FieldVector<int, 1> > > DataHandle;
    #endif

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

    //! Indices of matrix and rhs entries
    /**
234
235
236
    * 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:
237
238
239
240
241
242
243
    */
    enum
    {
        rhs = 1,//!<index for the right hand side entry
        matrix = 0//!<index for the global matrix entry

    };
244
245
246
};

//! function which assembles the system of equations to be solved
Benjamin Faigle's avatar
Benjamin Faigle committed
247
/*! for first == true, this function assembles the matrix and right hand side for
248
249
250
251
252
253
254
255
256
257
 * 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)
258
259
260
261
262
263
264
{
    if(first)
    {
        ParentType::assemble(true);
        return;
    }
    // initialization: set matrix this->A_ to zero
265
266
    this->A_ = 0;
    this->f_ = 0;
267

268
269
    ElementIterator eEndIt = problem().gridView().template end<0> ();
    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eEndIt; ++eIt)
270
    {
271
        // get the global index of the cell
272
        int eIdxGlobalI = problem().variables().index(*eIt);
273

274
275
        // assemble interior element contributions
        if (eIt->partitionType() == Dune::InteriorEntity)
276
        {
277
            // get the cell data
278
            CellData& cellDataI = problem().variables().cellData(eIdxGlobalI);
279

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

282
283
284
285
286
            /*****  source term ***********/
            if(cellDataI.subdomain() != 2)
                problem().pressureModel().get1pSource(entries,*eIt, cellDataI);
            else
                problem().pressureModel().getSource(entries,*eIt, cellDataI, first);
287

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

290
291
            /*****  flux term ***********/
            // iterate over all faces of the cell
292
293
            IntersectionIterator isEndIt = problem().gridView().iend(*eIt);
            for (IntersectionIterator isIt = problem().gridView().ibegin(*eIt); isIt != isEndIt; ++isIt)
294
            {
295
296
297
                /************* handle interior face *****************/
                if (isIt->neighbor())
                {
298
                    int eIdxGlobalJ = problem().variables().index(*(isIt->outside()));
299
300

                    if (cellDataI.subdomain() != 2
301
                            or problem().variables().cellData(eIdxGlobalJ).subdomain() != 2) // cell in the 1p domain
302
303
304
305
306
                        get1pFlux(entries, *isIt, cellDataI);
                    else
                        problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);

                    //set right hand side
307
                    this->f_[eIdxGlobalI] -= entries[rhs];
308
                    // set diagonal entry
309
                    this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
310
                    // set off-diagonal entry
311
                    this->A_[eIdxGlobalI][eIdxGlobalJ] = -entries[matrix];
312
313
314
315
                }   // end neighbor


                /************* boundary face ************************/
316
                else
317
318
                {
                    if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
319
                        problem().pressureModel().get1pFluxOnBoundary(entries, *isIt, cellDataI);
320
321
322
323
                    else
                        problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);

                    //set right hand side
324
                    this->f_[eIdxGlobalI] += entries[rhs];
325
                    // set diagonal entry
326
                    this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
327
328
329
330
331
332
                }
            } //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
333
                problem().pressureModel().get1pStorage(entries, *eIt, cellDataI);
334
335
            else
                problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
336

337
            this->f_[eIdxGlobalI] += entries[rhs];
338
            // set diagonal entry
339
            this->A_[eIdxGlobalI][eIdxGlobalI] += entries[matrix];
340
341
342
343
        }
        // assemble overlap and ghost element contributions
        else 
        {
344
345
346
            this->A_[eIdxGlobalI] = 0.0;
            this->A_[eIdxGlobalI][eIdxGlobalI] = 1.0;
            this->f_[eIdxGlobalI] = this->pressure()[eIdxGlobalI];
347
        }
348
349
350
351
352
    } // 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;
}
353

354
355
356
357
358
359
360
361
362
//! 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
 */
363
template<class TypeTag>
364
365
void FVPressure2P2CMultiPhysics<TypeTag>::get1pSource(Dune::FieldVector<Scalar, 2>& sourceEntry,
        const Element& elementI, const CellData& cellDataI)
366
367
{
    sourceEntry=0.;
368

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

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

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

380
381
    return;
}
382

383
384
385
386
387
388
389
390
391
392
393
//! 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
 */
394
395
396
397
398
399
400
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pStorage(Dune::FieldVector<Scalar, 2>& storageEntry,
                                                        const Element& elementI,
                                                        CellData& cellDataI)
{
    storageEntry = 0.;
    // cell index
401
    int eIdxGlobalI = problem().variables().index(elementI);
402
403
    int presentPhaseIdx = cellDataI.subdomain();
    Scalar volume = elementI.geometry().volume();
404

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

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

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

418
419
420
421
        Scalar sumC = (cellDataI.massConcentration(wCompIdx) + cellDataI.massConcentration(nCompIdx));
        Scalar Z1 = cellDataI.massConcentration(wCompIdx) / sumC;
        // initialize simple fluidstate object
        PseudoOnePTwoCFluidState<TypeTag> pseudoFluidState;
Benjamin Faigle's avatar
Benjamin Faigle committed
422
423
        CompositionalFlash<TypeTag> flashSolver;
        flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
424
                cellDataI.temperature(wPhaseIdx));
425

Benjamin Faigle's avatar
Benjamin Faigle committed
426
        Scalar v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
427
        cellDataI.dv_dp() = (sumC * ( v_ - (1. /cellDataI.density(presentPhaseIdx)))) /incp;
428

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

            p_ -= 2*incp;
Benjamin Faigle's avatar
Benjamin Faigle committed
435
            flashSolver.concentrationFlash1p2c(pseudoFluidState, Z1, p_, cellDataI.subdomain(),
436
                    cellDataI.temperature(wPhaseIdx));
Benjamin Faigle's avatar
Benjamin Faigle committed
437
            v_ = 1. / pseudoFluidState.density(presentPhaseIdx);
438
439
440
441
            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)
            {
442
443
                Dune::dinfo <<__FILE__<< "dv_dp still larger 0 after inverting secant. regularize"<< std::endl;
                cellDataI.dv_dp() *= -1;
444
445
            }
        }
446
447


448
        Scalar compress_term = cellDataI.dv_dp() / timestep_;
449

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

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

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

460
461
462
463
464
    // 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.
    if( problem().timeManager().episodeWillBeOver()
            || problem().timeManager().willBeFinished())
    {
465
        problem().variables().cellData(eIdxGlobalI).errorCorrection() = 0.;
466
467
        return;
    }
468

469
470
    // error reduction routine: volumetric error is damped and inserted to right hand side
    // if damping is not done, the solution method gets unstable!
471
    problem().variables().cellData(eIdxGlobalI).volumeError() /= timestep_;
472
473
474
475
476
477
478
479
480
481
482
483
484
    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)
485
            storageEntry[rhs] +=
486
                    problem().variables().cellData(eIdxGlobalI).errorCorrection() =
487
488
489
                            fac* (1-x_mi*(lofac-1)/(x_lo-x_mi) + (lofac-1)/(x_lo-x_mi)*erri/maxError)
                                * cellDataI.volumeError() * volume;
        else
490
            storageEntry[rhs] +=
491
                    problem().variables().cellData(eIdxGlobalI).errorCorrection() =
492
493
494
                            fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError)
                                * cellDataI.volumeError() * volume;
    }
Benjamin Faigle's avatar
Benjamin Faigle committed
495
    else
496
        problem().variables().cellData(eIdxGlobalI).errorCorrection()=0 ;
497

498
499
    return;
}
500

501
502
503
504
505
506
/*! \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}
507
       \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) .
508
509
510
511
512
513
   \f]
 *
 * \param entries The Matrix and RHS entries
 * \param intersection Intersection between cell I and J
 * \param cellDataI Data of cell I
 */
514
515
516
517
518
519
template<class TypeTag>
void FVPressure2P2CMultiPhysics<TypeTag>::get1pFlux(Dune::FieldVector<Scalar, 2>& entries,
        const Intersection& intersection, const CellData& cellDataI)
{
    entries = 0.;
    ElementPointer elementPointerI = intersection.inside();
520

521
522
    // get global coordinate of cell center
    const GlobalPosition& globalPos = elementPointerI->geometry().center();
523

524
    // cell index
525
//    int eIdxGlobalI = problem().variables().index(*elementPointerI);
526

527
    // get absolute permeability
528
    DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPointerI));
529

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

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

536
537
        // access neighbor
        ElementPointer neighborPointer = intersection.outside();
538
539
        int eIdxGlobalJ = problem().variables().index(*neighborPointer);
        CellData& cellDataJ = problem().variables().cellData(eIdxGlobalJ);
540
541
542
543
544
545
546
547
548
549
550
551
552

        // gemotry info of neighbor
        const GlobalPosition& globalPosNeighbor = neighborPointer->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;

553
554
        DimMatrix permeabilityJ
            = problem().spatialParams().intrinsicPermeability(*neighborPointer);
555
556

        // compute vectorized permeabilities
557
        DimMatrix meanPermeability(0);
558
559
560
561
562
563
564
565
566
567
568
        Dumux::harmonicMeanMatrix(meanPermeability, permeabilityI, permeabilityJ);

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

569
        // 1p => no pc => only 1 pressure, potential
570
        Scalar potential = (cellDataI.pressure(phaseIdx) - cellDataJ.pressure(phaseIdx)) / dist;
571
572
573
574
575

        potential += rhoMean * (unitDistVec * gravity_);

        Scalar lambda;

Benjamin Faigle's avatar
Benjamin Faigle committed
576
        if (potential > 0.)
577
        {
578
            lambda = cellDataI.mobility(phaseIdx);
Benjamin Faigle's avatar
Benjamin Faigle committed
579
580
            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
581
        }
Benjamin Faigle's avatar
Benjamin Faigle committed
582
        else if (potential < 0.)
583
        {
584
            lambda = cellDataJ.mobility(phaseIdx);
Benjamin Faigle's avatar
Benjamin Faigle committed
585
586
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, true);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, true);
Benjamin Faigle's avatar
Benjamin Faigle committed
587
588
589
590
        }
        else
        {
            lambda = harmonicMean(cellDataI.mobility(phaseIdx) , cellDataJ.mobility(phaseIdx));
Benjamin Faigle's avatar
Benjamin Faigle committed
591
592
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
            cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
593
        }
594

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

599
600
601
        return;
}

602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
/*! \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
 */
618
619
620
621
622
623
624
625
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
    ElementPointer elementPointerI = intersection.inside();
    const GlobalPosition& globalPos = elementPointerI->geometry().center();
626
//    int eIdxGlobalI = problem().variables().index(*elementPointerI);
627
628
629
630
631
632
    int phaseIdx = cellDataI.subdomain();

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

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

                // geometrical information
638
                GlobalPosition distVec(globalPosFace - globalPos);
639
                Scalar dist = distVec.two_norm();
640
                GlobalPosition unitDistVec(distVec);
641
642
643
                unitDistVec /= dist;

                //get boundary condition for boundary face center
644
645
                BoundaryTypes bcType;
                problem().boundaryTypes(bcType, intersection);
646
647
648

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

                /**********         Dirichlet Boundary        *************/
651
                if (bcType.isDirichlet(Indices::pressureEqIdx))
652
                {
653
                    // get absolute permeability
654
                    DimMatrix permeabilityI(problem().spatialParams().intrinsicPermeability(*elementPointerI));
655
656
657
658
659
660
                    if(this->regulateBoundaryPermeability)
                    {
                        int axis = intersection.indexInInside() / 2;
                        if(permeabilityI[axis][axis] < this->minimalBoundaryPermeability)
                            permeabilityI[axis][axis] = this->minimalBoundaryPermeability;
                    }
661
662
663
                    // get mobilities and fractional flow factors
                    Scalar lambdaI = cellDataI.mobility(phaseIdx);

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

668
669
                    // create a fluid state for the boundary
                    FluidState BCfluidState;
670

671
672
                    //read boundary values
                    PrimaryVariables primaryVariablesOnBoundary(NAN);
673
                    problem().dirichlet(primaryVariablesOnBoundary, intersection);
674
675

                    {
676
677
                        // read boundary values
                        problem().transportModel().evalBoundary(globalPosFace,
678
                                                                    intersection,
679
680
                                                                    BCfluidState,
                                                                    pressBC);
681

682
                        // determine fluid properties at the boundary
683
684
685
686
                        Scalar densityBound =
                            FluidSystem::density(BCfluidState, phaseIdx);
                        Scalar viscosityBound =
                            FluidSystem::viscosity(BCfluidState, phaseIdx);
687
688

                        // mobility at the boundary
689
                        Scalar lambdaBound = 0.;
690
                        switch (GET_PROP_VALUE(TypeTag, BoundaryMobility))
691
692
693
                        {
                        case Indices::satDependent:
                            {
694
695
                                lambdaBound = BCfluidState.saturation(phaseIdx)
                                    / viscosityBound;
696
697
698
699
                            break;
                            }
                        case Indices::permDependent:
                            {
700
701
                            if (phaseIdx == wPhaseIdx)
                                lambdaBound = MaterialLaw::krw(
702
                                    problem().spatialParams().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
703
704
705
                                    / viscosityBound;
                            else
                                lambdaBound = MaterialLaw::krn(
706
                                    problem().spatialParams().materialLawParams(*elementPointerI), BCfluidState.saturation(wPhaseIdx))
707
                                    / viscosityBound;
708
                            break;
709
710
                            }
                        }
711
                        Scalar rhoMean = 0.5 * (cellDataI.density(phaseIdx) + densityBound);
712

713
                        Scalar potential = 0;
714

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

718
                        potential += rhoMean * (unitDistVec * gravity_);
719

720
                        Scalar lambda(0.);
721

Bernd Flemisch's avatar
Bernd Flemisch committed
722
                        if (Dune::FloatCmp::eq<Scalar, Dune::FloatCmp::absolute>(potential, 0.0, 1.0e-30))
723
                        {
Bernd Flemisch's avatar
Bernd Flemisch committed
724
725
726
727
728
                            lambda = 0.5*(lambdaI + lambdaBound);
                        }
                        else if (potential > 0.)
                        {
                            lambda = lambdaI;
729
                        }
730
                        else
731
                        {
732
733
                            lambda = lambdaBound;
                        }
734

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

739
                        //calculate right hand side
740
                        rightEntry = lambda * rhoMean * fabs(permeability * unitOuterNormal)
741
                                * faceArea ;
742
743

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

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

759
                    entries[1] -= J[1+phaseIdx] * faceArea;
760
                }
761
762
                else
                    DUNE_THROW(Dune::NotImplemented, "Boundary Condition neither Dirichlet nor Neumann!");
763
764
765
766
767


    return;
}

768

769
770
//! constitutive functions are updated once if new concentrations are calculated and stored in the variables container
/*!
771
 * In contrast to the standard sequential 2p2c model ( FVPressure2P2C<TypeTag>::updateMaterialLaws() ),
772
 * this method also holds routines to adapt the subdomain. The subdomain indicates weather we are in 1p domain (value = 1)
773
 * or in the two phase subdomain (value = 2).
774
775
776
777
 * 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()
778
779
 */
template<class TypeTag>
Christoph Grueninger's avatar
Christoph Grueninger committed
780
void FVPressure2P2CMultiPhysics<TypeTag>::updateMaterialLaws(bool postTimeStep)
781
782
{
    //get timestep for error term
783
    Scalar maxError = 0.;
784
785

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

791
    // Loop A) through leaf grid
792
793
    ElementIterator eEndIt = problem().gridView().template end<0> ();
    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eEndIt; ++eIt)
794
795
    {
        // get global coordinate of cell center
796
797
        int eIdxGlobal = problem().variables().index(*eIt);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
798

799
        if(cellData.subdomain() == 2)    // complex
800
        {
801
            this->updateMaterialLawsInElement(*eIt, postTimeStep);
802
803

            // check subdomain consistency
804
            timer_.start();
805
806
807
808
809
            // enshure we are not at source
            // get sources from problem
            PrimaryVariables source(NAN);
            problem().source(source, *eIt);

810
811
            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
812
813
            {
                // mark this element
814
                nextSubdomain[eIdxGlobal] = 2;
815
816

                // mark neighbors
817
818
                IntersectionIterator isEndIt = problem().gridView().iend(*eIt);
                for (IntersectionIterator isIt = problem().gridView().ibegin(*eIt); isIt!=isEndIt; ++isIt)
819
820
821
                {
                    if (isIt->neighbor())
                    {
822
                        int eIdxGlobalJ = problem().variables().index(*(isIt->outside()));
823
                        // mark neighbor Element
824
                        nextSubdomain[eIdxGlobalJ] = 2;
825
826
                    }
                }
827
            }
828
            else if(nextSubdomain[eIdxGlobal] != 2)// update next subdomain if possible
829
            {
830
                if(Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(wPhaseIdx), 0.0, 1.0e-30))
831
                    nextSubdomain[eIdxGlobal] = wPhaseIdx;
832
                else if (Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(cellData.saturation(nPhaseIdx), 0.0, 1.0e-30))
833
                    nextSubdomain[eIdxGlobal] = nPhaseIdx;
834
            }
835
            timer_.stop();
836
            // end subdomain check
837
        }// end complex domain
838
839
        else if (nextSubdomain[eIdxGlobal] != 2) //check if cell remains in simple subdomain
            nextSubdomain[eIdxGlobal] = cellData.subdomain();
840

841
    } //end define complex area of next subdomain
842

843
    timer_.start();
844
845
846
847
848
849
850
851
852
    //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

853
854
    // Loop B) thorugh leaf grid
    // investigate cells that were "simple" in current TS
855
    for (ElementIterator eIt = problem().gridView().template begin<0> (); eIt != eEndIt; ++eIt)
856
    {
857
858
        int eIdxGlobal = problem().variables().index(*eIt);
        CellData& cellData = problem().variables().cellData(eIdxGlobal);
859

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

864
865
        //first check if simple will become complicated
        if(oldSubdomainI != 2
866
                    && nextSubdomain[eIdxGlobal] == 2)
867
        {
Benjamin Faigle's avatar
Benjamin Faigle committed
868
869
            // use complex update of the fluidstate
            timer_.stop();
870
            this->updateMaterialLawsInElement(*eIt, postTimeStep);
871
872
873
            timer_.start();
        }
        else if(oldSubdomainI != 2
874
                    && nextSubdomain[eIdxGlobal] != 2)    // will be simple and was simple
875
        {
Benjamin Faigle's avatar
Benjamin Faigle committed
876
			// perform simple update
877
            this->update1pMaterialLawsInElement(*eIt, cellData, postTimeStep);
878
        }
879
880
881
        //else
        // a) will remain complex -> everything already done in loop A
        // b) will be simple and was complex: complex FS available, so next TS
882
        //         will use comlex FS, next updateMaterialLaw will transform to simple
883

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

888
    timer_.stop();
Benjamin Faigle's avatar
Benjamin Faigle committed
889
890
891
    
    if(problem().timeManager().willBeFinished() or problem().timeManager().episodeWillBeOver())
    	Dune::dinfo << "Subdomain routines took "<