2p2cmodel.hh 29.5 KB
Newer Older
1
// $Id$
2
/*****************************************************************************
3
 *   Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, Bernd Flemisch     *
4
5
6
7
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
8
 *   This program is free software: you can redistribute it and/or modify    *
9
 *   it under the terms of the GNU General Public License as published by    *
10
11
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
12
 *                                                                           *
13
14
15
16
17
18
19
 *   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/>.   *
20
 *****************************************************************************/
Melanie Darcis's avatar
doc    
Melanie Darcis committed
21
22
23
24
25
26
27

/*!
* \file
*
* \brief Adaption of the BOX scheme to the two-phase two-component flow model.
*/

28
29
#ifndef DUMUX_2P2C_MODEL_HH
#define DUMUX_2P2C_MODEL_HH
30

31
#include "2p2cproperties.hh"
32
33
#include "2p2clocalresidual.hh"
#include "2p2cproblem.hh"
34
35
36

namespace Dumux
{
37
38
39
40
41
42
43
44
45
46
47
48
49
50
/*!
 * \ingroup BoxModels
 * \defgroup TwoPTwoCModel Two-phase two-component box model
 */

/*!
 * \ingroup TwoPTwoCModel
 * \brief Adaption of the BOX scheme to the two-phase two-component flow model.
 *
 * This model implements two-phase two-component flow of two compressible and
 * partially miscible fluids \f$\alpha \in \{ w, n \}\f$ composed of the two components
 * \f$\kappa \in \{ w, a \}\f$. The standard multiphase Darcy
 * approach is used as the equation for the conservation of momentum:
 * \f[
51
52
 v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
 \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
53
 * \f]
54
 *
55
56
 * By inserting this into the equations for the conservation of the
 * components, one gets one transport equation for each component
57
 * \f{eqnarray}
58
 && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}
59
 {\partial t}
60
 - \sum_\alpha  \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
61
 \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
62
 (\text{grad}\, p_\alpha - \varrho_{\alpha}  \mbox{\bf g}) \right\}
63
 \nonumber \\ \nonumber \\
64
65
    &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\}
 - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, ,
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
 \alpha \in \{w, g\}
 \f}
 *
 * This is discretized using a fully-coupled vertex
 * centered finite volume (box) scheme as spatial and
 * the implicit Euler method as temporal discretization.
 *
 * By using constitutive relations for the capillary pressure \f$p_c =
 * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
 * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
 * unknowns can be reduced to two.
 * The used primary variables are, like in the two-phase model, either \f$p_w\f$ and \f$S_n\f$
 * or \f$p_n\f$ and \f$S_w\f$. The formulation which ought to be used can be
 * specified by setting the <tt>Formulation</tt> property to either
 * TwoPTwoCIndices::pWsN or TwoPTwoCIndices::pNsW. By
 * default, the model uses \f$p_w\f$ and \f$S_n\f$.
 * Moreover, the second primary variable depends on the phase state, since a
 * primary variable switch is included. The phase state is stored for all nodes
 * of the system. Following cases can be distinguished:
 * <ul>
 *  <li> Both phases are present: The saturation is used (either \f$S_n\f$ or \f$S_w\f$, dependent on the chosen <tt>Formulation</tt>),
 *      as long as \f$ 0 < S_\alpha < 1\f$</li>.
 *  <li> Only wetting phase is present: The mass fraction of, e.g., air in the wetting phase \f$X^a_w\f$ is used,
 *      as long as the maximum mass fraction is not exceeded (\f$X^a_w<X^a_{w,max}\f$)</li>
 *  <li> Only non-wetting phase is present: The mass fraction of, e.g., water in the non-wetting phase, \f$X^w_n\f$, is used,
 *      as long as the maximum mass fraction is not exceeded (\f$X^w_n<X^w_{n,max}\f$)</li>
 * </ul>
93
 */
94

95
template<class TypeTag>
96
class TwoPTwoCModel: public BoxModel<TypeTag>
97
{
98
99
    typedef TwoPTwoCModel<TypeTag> ThisType;
    typedef BoxModel<TypeTag> ParentType;
100
101

    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
102
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
103
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
104
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
105

106
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
Andreas Lauser's avatar
Andreas Lauser committed
107
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
108
109
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
110
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
111
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
112
113
114
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
115
116
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;

117
    enum {
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
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld,

        numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
        numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
        numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),

        pressureIdx = Indices::pressureIdx,
        switchIdx = Indices::switchIdx,

        contiLEqIdx = Indices::contiLEqIdx,
        contiGEqIdx = Indices::contiGEqIdx,

        lPhaseIdx = Indices::lPhaseIdx,
        gPhaseIdx = Indices::gPhaseIdx,

        lCompIdx = Indices::lCompIdx,
        gCompIdx = Indices::gCompIdx,

        lPhaseOnly = Indices::lPhaseOnly,
        gPhaseOnly = Indices::gPhaseOnly,
        bothPhases = Indices::bothPhases,

        plSg = TwoPTwoCFormulation::plSg,
        pgSl = TwoPTwoCFormulation::pgSl,
        formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
    };

146
147
    typedef TwoPTwoCFluidState<TypeTag> FluidState;

148
    typedef typename GridView::template Codim<dim>::Entity Vertex;
149
150
151
    typedef typename GridView::template Codim<0>::Entity Element;
    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
152

153
154
155
    typedef Dune::FieldVector<Scalar, dim> LocalPosition;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;

156
157
158
    static const Scalar mobilityUpwindAlpha =
            GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));

159
160
161
public:
    /*!
     * \brief Initialize the static data with the initial solution.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
162
163
     *
     * \param problem The problem to be solved
164
     */
165
    void init(Problem &problem)
166
    {
167
        ParentType::init(problem);
168

169
        staticVertexDat_.resize(this->gridView_().size(dim));
170

171
172
173
174
        setSwitched_(false);

        VertexIterator it = this->gridView_().template begin<dim> ();
        const VertexIterator &endit = this->gridView_().template end<dim> ();
175
176
        for (; it != endit; ++it)
        {
177
            int globalIdx = this->dofMapper().map(*it);
178
179
180
181
            const GlobalPosition &globalPos = it->geometry().corner(0);

            // initialize phase presence
            staticVertexDat_[globalIdx].phasePresence
182
                = this->problem_().initialPhasePresence(*it, globalIdx,
183
184
185
186
187
188
189
190
                            globalPos);
            staticVertexDat_[globalIdx].wasSwitched = false;

            staticVertexDat_[globalIdx].oldPhasePresence
                    = staticVertexDat_[globalIdx].phasePresence;
        }
    }

Andreas Lauser's avatar
Andreas Lauser committed
191
192
193
    /*!
     * \brief Compute the total storage inside one phase of all
     *        conservation quantities.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
194
195
196
     *
     * \param dest Contains the storage of each component for one phase
     * \param phaseIdx The phase index
Andreas Lauser's avatar
Andreas Lauser committed
197
198
199
200
201
202
203
204
205
206
207
208
209
210
     */
    void globalPhaseStorage(PrimaryVariables &dest, int phaseIdx)
    {
        dest = 0;

        ElementIterator elemIt = this->gridView_().template begin<0>();
        const ElementIterator elemEndIt = this->gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt) {
            this->localResidual().evalPhaseStorage(*elemIt, phaseIdx);

            for (int i = 0; i < elemIt->template count<dim>(); ++i)
                dest += this->localResidual().residual(i);
        };

Bernd Flemisch's avatar
Bernd Flemisch committed
211
212
        if (this->gridView_().comm().size() > 1)
        	dest = this->gridView_().comm().sum(dest);
Andreas Lauser's avatar
Andreas Lauser committed
213
    }
214
215

    /*!
216
217
     * \brief Called by the update() method if applying the newton
     *         method was unsuccessful.
218
     */
219
    void updateFailed()
220
    {
221
        ParentType::updateFailed();
222

223
224
        setSwitched_(false);
        resetPhasePresence_();
225
226
227
        /*this->localJacobian().updateStaticData(this->curSolFunction(),
         this->prevSolFunction());
         */
Andreas Lauser's avatar
Andreas Lauser committed
228
    };
229

230
231
232
    /*!
     * \brief Returns the relative weight of a primary variable for
     *        calculating relative errors.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
233
234
235
     *
     * \param globalVertexIdx The global vertex index
     * \param pvIdx The primary variable index
236
     */
237
    Scalar primaryVarWeight(int globalVertexIdx, int pvIdx) const
238
239
    {
        if (Indices::pressureIdx == pvIdx)
240
            return std::min(1.0/this->prevSol()[globalVertexIdx][pvIdx], 1.0);
241
242
243
        return 1;
    }

244
    /*!
Melanie Darcis's avatar
doc    
Melanie Darcis committed
245
     * \brief Called by the problem if a time integration was
246
247
     *        successful, post processing of the solution is done and the
     *        result has been written to disk.
248
     *
Melanie Darcis's avatar
doc    
Melanie Darcis committed
249
     * This should prepare the model for the next time integration.
250
     */
251
    void advanceTimeLevel()
252
    {
253
        ParentType::advanceTimeLevel();
254

255
        // update the phase state
256
257
        updateOldPhasePresence_();
        setSwitched_(false);
258
259
260
261
262
263
264
265
266
267
268
269
    }

    /*!
     * \brief Return true if the primary variables were switched for
     *        at least one vertex after the last timestep.
     */
    bool switched() const
    {
        return switchFlag_;
    }

    /*!
270
     * \brief Returns the phase presence of the current or the old solution of a vertex.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
271
272
273
     *
     * \param globalVertexIdx The global vertex index
     * \param oldSol Evaluate function with solution of current or previous time step
274
     */
275
    int phasePresence(int globalVertexIdx, bool oldSol) const
276
    {
277
278
        return oldSol ? staticVertexDat_[globalVertexIdx].oldPhasePresence
                : staticVertexDat_[globalVertexIdx].phasePresence;
279
280
281
    }

    /*!
282
283
284
     * \brief Append all quantities of interest which can be derived
     *        from the solution of the current time step to the VTK
     *        writer.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
285
286
287
     *
     * \param sol The solution vector
     * \param writer The writer for multi-file VTK datasets
288
289
     */
    template<class MultiWriter>
290
    void addOutputVtkFields(const SolutionVector &sol,
291
                            MultiWriter &writer)
292
293
294
295
    {
        typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;

        // create the required scalar fields
296
        unsigned numVertices = this->problem_().gridView().size(dim);
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315

        ScalarField *Sg = writer.template createField<Scalar, 1> (numVertices);
        ScalarField *Sl = writer.template createField<Scalar, 1> (numVertices);
        ScalarField *pg = writer.template createField<Scalar, 1> (numVertices);
        ScalarField *pl = writer.template createField<Scalar, 1> (numVertices);
        ScalarField *pc = writer.template createField<Scalar, 1> (numVertices);
        ScalarField *rhoL =
                writer.template createField<Scalar, 1> (numVertices);
        ScalarField *rhoG =
                writer.template createField<Scalar, 1> (numVertices);
        ScalarField *mobL =
                writer.template createField<Scalar, 1> (numVertices);
        ScalarField *mobG =
                writer.template createField<Scalar, 1> (numVertices);
        ScalarField *phasePresence = writer.template createField<Scalar, 1> (
                numVertices);
        ScalarField *massFrac[numPhases][numComponents];
        for (int i = 0; i < numPhases; ++i)
            for (int j = 0; j < numComponents; ++j)
316
317
318
                massFrac[i][j] = writer.template createField<Scalar, 1>(numVertices);
        ScalarField *temperature = writer.template createField<Scalar, 1>(numVertices);
        ScalarField *poro = writer.template createField<Scalar, 1>(numVertices);
319

320
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
321
322
323
        ScalarField *velocityX = writer.template createField<Scalar, 1>(numVertices);
        ScalarField *velocityY = writer.template createField<Scalar, 1>(numVertices);
        ScalarField *velocityZ = writer.template createField<Scalar, 1>(numVertices);
324
//        Scalar maxV=0.; // variable to store the maximum face velocity
325
326
327
328
329
330
331
332
333
334
335
336
337
338

        // initialize velocity fields
        Scalar boxSurface[numVertices];
        for (int i = 0; i < numVertices; ++i)
        {
            (*velocityX)[i] = 0;
            if (dim > 1)
            (*velocityY)[i] = 0;
            if (dim > 2)
            (*velocityZ)[i] = 0;
            boxSurface[i] = 0.0; // initialize the boundary surface of the fv-boxes
        }
#endif

339
        unsigned numElements = this->gridView_().size(0);
340
341
342
        ScalarField *rank =
                writer.template createField<Scalar, 1> (numElements);

343
        FVElementGeometry fvElemGeom;
344
        VolumeVariables volVars;
345

346
347
348
        ElementIterator elemIt = this->gridView_().template begin<0>();
        ElementIterator elemEndIt = this->gridView_().template end<0>();
        for (; elemIt != elemEndIt; ++elemIt)
349
        {
350
351
352
            int idx = this->problem_().elementMapper().map(*elemIt);
            (*rank)[idx] = this->gridView_().comm().rank();
            fvElemGeom.update(this->gridView_(), *elemIt);
353

354
355
            int numVerts = elemIt->template count<dim> ();
            for (int i = 0; i < numVerts; ++i)
356
            {
357
                int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
358
                volVars.update(sol[globalIdx],
359
360
                               this->problem_(),
                               *elemIt,
361
                               fvElemGeom,
362
363
                               i,
                               false);
364
365
366
367
368
369
370
371
372
                (*Sg)[globalIdx] = volVars.saturation(gPhaseIdx);
                (*Sl)[globalIdx] = volVars.saturation(lPhaseIdx);
                (*pg)[globalIdx] = volVars.pressure(gPhaseIdx);
                (*pl)[globalIdx] = volVars.pressure(lPhaseIdx);
                (*pc)[globalIdx] = volVars.capillaryPressure();
                (*rhoL)[globalIdx] = volVars.fluidState().density(lPhaseIdx);
                (*rhoG)[globalIdx] = volVars.fluidState().density(gPhaseIdx);
                (*mobL)[globalIdx] = volVars.mobility(lPhaseIdx);
                (*mobG)[globalIdx] = volVars.mobility(gPhaseIdx);
373
374
375
376
                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                    {
                        (*massFrac[phaseIdx][compIdx])[globalIdx]
377
                                = volVars.fluidState().massFrac(phaseIdx,
378
                                        compIdx);
379
380

                        Valgrind::CheckDefined(
381
                            (*massFrac[phaseIdx][compIdx])[globalIdx][0]);
382
                    }
383
384
                (*poro)[globalIdx] = volVars.porosity();
                (*temperature)[globalIdx] = volVars.temperature();
385
386
387
388
                (*phasePresence)[globalIdx]
                        = staticVertexDat_[globalIdx].phasePresence;
            };

389
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
390
391
392
            // In the box method, the velocity is evaluated on the FE-Grid. However, to get an
            // average apparent velocity at the vertex, all contributing velocities have to be interpolated.
            GlobalPosition velocity(0.);
393
394
395
396
397
398
399
            ElementVolumeVariables elemVolVars;

            elemVolVars.update(this->problem_(),
                               *elemIt,
                               fvElemGeom,
                               false /* isOldSol? */);

400
            // loop over the phases
401
            for (int faceIdx = 0; faceIdx< fvElemGeom.numEdges; faceIdx++)
402
403
            {
                //prepare the flux calculations (set up and prepare geometry, FE gradients)
404
                FluxVariables fluxDat(this->problem_(),
405
406
                                 *elemIt,
                                 fvElemGeom,
407
                                 faceIdx,
408
                                 elemVolVars);
409
410
411
412
413

                // choose phase of interest. Alternatively, a loop over all phases would be possible.
                int phaseIdx = gPhaseIdx;

                // get darcy velocity
414
                velocity = fluxDat.Kmvp(phaseIdx); // mind the sign: vDarcy = kf grad p
415
416

                // up+downstream mobility
417
418
                const VolumeVariables &up = elemVolVars[fluxDat.upstreamIdx(phaseIdx)];
                const VolumeVariables &down = elemVolVars[fluxDat.downstreamIdx(phaseIdx)];
419
420
421
                Scalar scvfArea = fluxDat.face().normal.two_norm(); //get surface area to weight velocity at the IP with the surface area
                velocity *= (mobilityUpwindAlpha*up.mobility(phaseIdx) + (1-mobilityUpwindAlpha)*down.mobility(phaseIdx))* scvfArea;

422
                int vertIIdx = this->problem_().vertexMapper().map(*elemIt,
423
424
                        fluxDat.face().i,
                        dim);
425
                int vertJIdx = this->problem_().vertexMapper().map(*elemIt,
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
                        fluxDat.face().j,
                        dim);
                // add surface area for weighting purposes
                boxSurface[vertIIdx] += scvfArea;
                boxSurface[vertJIdx] += scvfArea;

                // Add velocity to upstream and downstream vertex.
                // Beware: velocity has to be substracted because of the (wrong) sign of vDarcy
                (*velocityX)[vertIIdx] -= velocity[0];
                (*velocityX)[vertJIdx] -= velocity[0];
                if (dim >= 2)
                {
                    (*velocityY)[vertIIdx] -= velocity[1];
                    (*velocityY)[vertJIdx] -= velocity[1];
                }
                if (dim == 3)
                {
                    (*velocityZ)[vertIIdx] -= velocity[2];
                    (*velocityZ)[vertJIdx] -= velocity[2];
                }
            }
#endif
        }

450
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
        // normalize the velocities at the vertices
        for (int i = 0; i < numVertices; ++i)
        {
            (*velocityX)[i] /= boxSurface[i];
            if (dim >= 2)
            (*velocityY)[i] /= boxSurface[i];
            if (dim == 3)
            (*velocityZ)[i] /= boxSurface[i];
        }
#endif

        writer.addVertexData(Sg, "Sg");
        writer.addVertexData(Sl, "Sl");
        writer.addVertexData(pg, "pg");
        writer.addVertexData(pl, "pl");
        writer.addVertexData(pc, "pc");
        writer.addVertexData(rhoL, "rhoL");
        writer.addVertexData(rhoG, "rhoG");
        writer.addVertexData(mobL, "mobL");
        writer.addVertexData(mobG, "mobG");
        for (int i = 0; i < numPhases; ++i)
        {
            for (int j = 0; j < numComponents; ++j)
            {
                std::string name = (boost::format("X_%s%s")
                        % ((i == lPhaseIdx) ? "l" : "g")
                        % FluidSystem::componentName(j)).str();
                writer.addVertexData(massFrac[i][j], name.c_str());
            }
        }
481
        writer.addVertexData(poro, "porosity");
482
483
484
        writer.addVertexData(temperature, "temperature");
        writer.addVertexData(phasePresence, "phase presence");

485
#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
486
487
488
489
490
491
492
493
494
495
        writer.addVertexData(velocityX, "Vx");
        if (dim >= 2)
        writer.addVertexData(velocityY, "Vy");
        if (dim == 3)
        writer.addVertexData(velocityZ, "Vz");
#endif
        writer.addCellData(rank, "process rank");
    }

    /*!
496
     * \brief Write the current solution to a restart file.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
497
498
499
     *
     * \param outStream The output stream of one vertex for the restart file
     * \param vert The vertex
500
     */
501
    void serializeEntity(std::ostream &outStream, const Vertex &vert)
502
    {
503
504
        // write primary variables
        ParentType::serializeEntity(outStream, vert);
505

506
507
508
        int vertIdx = this->dofMapper().map(vert);
        if (!outStream.good())
            DUNE_THROW(Dune::IOError, "Could not serialize vertex " << vertIdx);
509

510
        outStream << staticVertexDat_[vertIdx].phasePresence << " ";
511
512
513
514
515
    }

    /*!
     * \brief Reads the current solution for a vertex from a restart
     *        file.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
516
517
518
     *
     * \param inStream The input stream of one vertex from the restart file
     * \param vert The vertex
519
520
521
     */
    void deserializeEntity(std::istream &inStream, const Vertex &vert)
    {
522
523
        // read primary variables
        ParentType::deserializeEntity(inStream, vert);
524
525

        // read phase presence
526
        int vertIdx = this->dofMapper().map(vert);
527
        if (!inStream.good())
528
529
            DUNE_THROW(Dune::IOError,
                       "Could not deserialize vertex " << vertIdx);
530
531
532
533

        inStream >> staticVertexDat_[vertIdx].phasePresence;
        staticVertexDat_[vertIdx].oldPhasePresence
                = staticVertexDat_[vertIdx].phasePresence;
534

535
536
537
    }

    /*!
538
     * \brief Update the static data of all vertices in the grid.
Melanie Darcis's avatar
doc    
Melanie Darcis committed
539
540
541
     *
     * \param curGlobalSol The current global solution
     * \param oldGlobalSol The previous global solution
542
     */
543
    void updateStaticData(SolutionVector &curGlobalSol,
544
                          const SolutionVector &oldGlobalSol)
545
    {
546
547
548
549
550
551
        bool wasSwitched = false;

        for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
            staticVertexDat_[i].visited = false;

        FVElementGeometry fvElemGeom;
552
        static VolumeVariables volVars;
553
554
555
        ElementIterator it = this->gridView_().template begin<0> ();
        const ElementIterator &endit = this->gridView_().template end<0> ();
        for (; it != endit; ++it)
556
        {
557
558
559
560
561
562
563
564
565
            fvElemGeom.update(this->gridView_(), *it);
            for (int i = 0; i < fvElemGeom.numVertices; ++i)
            {
                int globalIdx = this->vertexMapper().map(*it, i, dim);

                if (staticVertexDat_[globalIdx].visited)
                    continue;

                staticVertexDat_[globalIdx].visited = true;
566
                volVars.update(curGlobalSol[globalIdx],
567
568
569
570
571
572
                               this->problem_(),
                               *it,
                               fvElemGeom,
                               i,
                               false);
                const GlobalPosition &global = it->geometry().corner(i);
573
574
                if (primaryVarSwitch_(curGlobalSol,
                                      volVars,
575
                                      globalIdx,
576
577
                                      global))
                    wasSwitched = true;
578
            }
579
580
        }

581
582
583
        // make sure that if there was a variable switch in an
        // other partition we will also set the switch flag
        // for our partition.
Bernd Flemisch's avatar
Bernd Flemisch committed
584
585
        if (this->gridView_().comm().size() > 1)
        	wasSwitched = this->gridView_().comm().max(wasSwitched);
586
587

        setSwitched_(wasSwitched);
588
589
590
    }

protected:
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
    /*!
     * \brief Data which is attached to each vertex and is not only
     *        stored locally.
     */
    struct StaticVars
    {
        int phasePresence;
        bool wasSwitched;

        int oldPhasePresence;
        bool visited;
    };

    /*!
     * \brief Reset the current phase presence of all vertices to the old one.
     *
     * This is done after an update failed.
     */
    void resetPhasePresence_()
610
    {
611
612
613
614
615
616
617
        int numVertices = this->gridView_().size(dim);
        for (int i = 0; i < numVertices; ++i)
        {
            staticVertexDat_[i].phasePresence
                    = staticVertexDat_[i].oldPhasePresence;
            staticVertexDat_[i].wasSwitched = false;
        }
618
    }
619
620
621
622
623

    /*!
     * \brief Set the old phase of all verts state to the current one.
     */
    void updateOldPhasePresence_()
624
    {
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
        int numVertices = this->gridView_().size(dim);
        for (int i = 0; i < numVertices; ++i)
        {
            staticVertexDat_[i].oldPhasePresence
                    = staticVertexDat_[i].phasePresence;
            staticVertexDat_[i].wasSwitched = false;
        }
    }

    /*!
     * \brief Set whether there was a primary variable switch after in
     *        the last timestep.
     */
    void setSwitched_(bool yesno)
    {
        switchFlag_ = yesno;
641
642
643
644
645
    }

    //  perform variable switch at a vertex; Returns true if a
    //  variable switch was performed.
    bool primaryVarSwitch_(SolutionVector &globalSol,
646
                           const VolumeVariables &volVars, int globalIdx,
647
                           const GlobalPosition &globalPos)
648
649
650
651
652
653
654
655
656
657
    {
        // evaluate primary variable switch
        bool wouldSwitch = false;
        int phasePresence = staticVertexDat_[globalIdx].phasePresence;
        int newPhasePresence = phasePresence;

        // check if a primary var switch is necessary
        if (phasePresence == gPhaseOnly)
        {
            // calculate mole fraction in the hypothetic liquid phase
658
659
            Scalar xll = volVars.fluidState().moleFrac(lPhaseIdx, lCompIdx);
            Scalar xlg = volVars.fluidState().moleFrac(lPhaseIdx, gCompIdx);
660
661
662
663
664
665
666
667
668
669
670
671

            Scalar xlMax = 1.0;
            if (xll + xlg > xlMax)
                wouldSwitch = true;
            if (staticVertexDat_[globalIdx].wasSwitched)
                xlMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, liquid phase appears
            if (xll + xlg > xlMax)
            {
                // liquid phase appears
672
                std::cout << "liquid phase appears at vertex " << globalIdx
673
674
675
676
                        << ", coordinates: " << globalPos << ", xll + xlg: "
                        << xll + xlg << std::endl;
                newPhasePresence = bothPhases;
                if (formulation == pgSl)
677
                    globalSol[globalIdx][switchIdx] = 0.0;
678
                else if (formulation == plSg)
679
                    globalSol[globalIdx][switchIdx] = 1.0;
680
681
682
683
684
685
            };
        }
        else if (phasePresence == lPhaseOnly)
        {
            // calculate fractions of the partial pressures in the
            // hypothetic gas phase
686
687
            Scalar xgl = volVars.fluidState().moleFrac(gPhaseIdx, lCompIdx);
            Scalar xgg = volVars.fluidState().moleFrac(gPhaseIdx, gCompIdx);
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

            Scalar xgMax = 1.0;
            if (xgl + xgg > xgMax)
                wouldSwitch = true;
            if (staticVertexDat_[globalIdx].wasSwitched)
                xgMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, gas phase appears
            if (xgl + xgg > xgMax)
            {
                // gas phase appears
                std::cout << "gas phase appears at vertex " << globalIdx
                        << ", coordinates: " << globalPos << ", x_gl + x_gg: "
                        << xgl + xgg << std::endl;
                newPhasePresence = bothPhases;
                if (formulation == pgSl)
                    globalSol[globalIdx][switchIdx] = 0.999;
                else if (formulation == plSg)
                    globalSol[globalIdx][switchIdx] = 0.001;
            }
        }
        else if (phasePresence == bothPhases)
        {
            Scalar Smin = 0.0;
            if (staticVertexDat_[globalIdx].wasSwitched)
                Smin = -0.01;

716
            if (volVars.saturation(gPhaseIdx) <= Smin)
717
718
719
720
721
            {
                wouldSwitch = true;
                // gas phase disappears
                std::cout << "Gas phase disappears at vertex " << globalIdx
                        << ", coordinates: " << globalPos << ", Sg: "
722
                        << volVars.saturation(gPhaseIdx) << std::endl;
723
724
725
                newPhasePresence = lPhaseOnly;

                globalSol[globalIdx][switchIdx]
726
                        = volVars.fluidState().massFrac(lPhaseIdx, gCompIdx);
727
            }
728
            else if (volVars.saturation(lPhaseIdx) <= Smin)
729
730
731
732
733
            {
                wouldSwitch = true;
                // liquid phase disappears
                std::cout << "Liquid phase disappears at vertex " << globalIdx
                        << ", coordinates: " << globalPos << ", Sl: "
734
                        << volVars.saturation(lPhaseIdx) << std::endl;
735
736
737
                newPhasePresence = gPhaseOnly;

                globalSol[globalIdx][switchIdx]
738
                        = volVars.fluidState().massFrac(gPhaseIdx, lCompIdx);
739
740
741
742
743
744
745
746
747
            }
        }

        staticVertexDat_[globalIdx].phasePresence = newPhasePresence;
        staticVertexDat_[globalIdx].wasSwitched = wouldSwitch;
        return phasePresence != newPhasePresence;
    }

    // parameters given in constructor
748
    std::vector<StaticVars> staticVertexDat_;
749
750
751
    bool switchFlag_;
};

752
}
753

754
755
#include "2p2cpropertydefaults.hh"

756
#endif