3p3cmodel.hh 41.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// -*- 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/>.   *
 *****************************************************************************/
/*!
 * \file
 *
22
23
 * \brief Adaption of the fully implicit scheme to the three-phase three-component
 *        flow model.
24
25
26
27
28
29
30
 *
 * The model is designed for simulating three fluid phases with water, gas, and
 * a liquid contaminant (NAPL - non-aqueous phase liquid)
 */
#ifndef DUMUX_3P3C_MODEL_HH
#define DUMUX_3P3C_MODEL_HH

31
#include <dune/common/version.hh>
32
#include <dumux/implicit/common/implicitvelocityoutput.hh>
33
34
35
36
37
38
#include "3p3cproperties.hh"

namespace Dumux
{
/*!
 * \ingroup ThreePThreeCModel
39
40
 * \brief Adaption of the fully implicit scheme to the three-phase three-component
 *        flow model.
41
42
43
44
45
46
 *
 * This model implements three-phase three-component flow of three fluid phases
 * \f$\alpha \in \{ water, gas, NAPL \}\f$ each composed of up to three components
 * \f$\kappa \in \{ water, air, contaminant \}\f$. The standard multiphase Darcy
 * approach is used as the equation for the conservation of momentum:
 * \f[
47
 v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
Thomas Fetzer's avatar
Thomas Fetzer committed
48
 \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
49
50
51
52
53
 * \f]
 *
 * By inserting this into the equations for the conservation of the
 * components, one transport equation for each component is obtained as
 * \f{eqnarray*}
Thomas Fetzer's avatar
Thomas Fetzer committed
54
 && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa
55
56
 S_\alpha )}{\partial t}
 - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha}
57
 \varrho_\alpha x_\alpha^\kappa \mathbf{K}
Thomas Fetzer's avatar
Thomas Fetzer committed
58
 (\textbf{grad}\, p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\}
59
60
 \nonumber \\
 \nonumber \\
Thomas Fetzer's avatar
Thomas Fetzer committed
61
62
 && - \sum\limits_\alpha \text{div} \left\{ D_\text{pm}^\kappa \varrho_\alpha \frac{M^\kappa}{M_\alpha}
 \textbf{grad} x^\kappa_{\alpha} \right\}
63
64
65
66
67
 - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha
 \f}
 *
 * Note that these balance equations are molar.
 *
Thomas Fetzer's avatar
Thomas Fetzer committed
68
69
70
 * All equations are discretized using a vertex-centered finite volume (box)
 * or cell-centered finite volume scheme as spatial
 * and the implicit Euler method as time discretization.
71
72
73
74
75
76
77
78
 *
 * The model uses commonly applied auxiliary conditions like
 * \f$S_w + S_n + S_g = 1\f$ for the saturations and
 * \f$x^w_\alpha + x^a_\alpha + x^c_\alpha = 1\f$ for the mole fractions.
 * Furthermore, the phase pressures are related to each other via
 * capillary pressures between the fluid phases, which are functions of
 * the saturation, e.g. according to the approach of Parker et al.
 *
79
 * The used primary variables are dependent on the locally present fluid phases.
80
81
82
 * An adaptive primary variable switch is included. The phase state is stored for all nodes
 * of the system. The following cases can be distinguished:
 * <ul>
83
84
85
86
87
 *  <li> All three phases are present: Primary variables are two saturations \f$(S_w\f$ and \f$S_n)\f$,
 *       and a pressure, in this case \f$p_g\f$. </li>
 *  <li> Only the water phase is present: Primary variables are now the mole fractions of air and 
 *       contaminant in the water phase \f$(x_w^a\f$ and \f$x_w^c)\f$, as well as the gas pressure, which is,
 *       of course, in a case where only the water phase is present, just the same as the water pressure. </li>
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
 *  <li> Gas and NAPL phases are present: Primary variables \f$(S_n\f$, \f$x_g^w\f$, \f$p_g)\f$. </li>
 *  <li> Water and NAPL phases are present: Primary variables \f$(S_n\f$, \f$x_w^a\f$, \f$p_g)\f$. </li>
 *  <li> Only gas phase is present: Primary variables \f$(x_g^w\f$, \f$x_g^c\f$, \f$p_g)\f$. </li>
 *  <li> Water and gas phases are present: Primary variables \f$(S_w\f$, \f$x_w^g\f$, \f$p_g)\f$. </li>
 * </ul>
 */
template<class TypeTag>
class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
{
    typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
107
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
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
    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

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

        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),

        switch1Idx = Indices::switch1Idx,
        switch2Idx = Indices::switch2Idx,


        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx,
        gPhaseIdx = Indices::gPhaseIdx,

        wCompIdx = Indices::wCompIdx,
        nCompIdx = Indices::nCompIdx,
        gCompIdx = Indices::gCompIdx,

        threePhases = Indices::threePhases,
        wPhaseOnly  = Indices::wPhaseOnly,
        gnPhaseOnly = Indices::gnPhaseOnly,
        wnPhaseOnly = Indices::wnPhaseOnly,
        gPhaseOnly  = Indices::gPhaseOnly,
        wgPhaseOnly = Indices::wgPhaseOnly

    };


    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;

    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
144
    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
145
    enum { dofCodim = isBox ? dim : 0 };
146

147
148
149
150
151
152
153
154
155
156
public:
    /*!
     * \brief Initialize the static data with the initial solution.
     *
     * \param problem The problem to be solved
     */
    void init(Problem &problem)
    {
        ParentType::init(problem);

157
        staticDat_.resize(this->numDofs());
158
159
160

        setSwitched_(false);

161
        if (isBox)
162
        {
163
164
165
            VertexIterator vIt = this->gridView_().template begin<dim> ();
            const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
            for (; vIt != vEndIt; ++vIt)
166
            {
167
168
169
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                int vIdxGlobal = this->dofMapper().index(*vIt);
#else
170
                int vIdxGlobal = this->dofMapper().map(*vIt);
171
#endif
172
                const GlobalPosition &globalPos = vIt->geometry().corner(0);
173

174
                // initialize phase presence
175
176
                staticDat_[vIdxGlobal].phasePresence
                    = this->problem_().initialPhasePresence(*vIt, vIdxGlobal,
177
                                                        globalPos);
178
                staticDat_[vIdxGlobal].wasSwitched = false;
179

180
181
                staticDat_[vIdxGlobal].oldPhasePresence
                    = staticDat_[vIdxGlobal].phasePresence;
182
183
184
185
            }
        }
        else
        {
186
187
188
            ElementIterator eIt = this->gridView_().template begin<0> ();
            const ElementIterator &eEndIt = this->gridView_().template end<0> ();
            for (; eIt != eEndIt; ++eIt)
189
            {
190
191
192
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                int eIdxGlobal = this->dofMapper().index(*eIt);
#else
193
                int eIdxGlobal = this->dofMapper().map(*eIt);
194
#endif
195
                const GlobalPosition &globalPos = eIt->geometry().center();
196
197

                // initialize phase presence
198
                staticDat_[eIdxGlobal].phasePresence
199
                    = this->problem_().initialPhasePresence(*this->gridView_().template begin<dim> (), 
200
201
                                                            eIdxGlobal, globalPos);
                staticDat_[eIdxGlobal].wasSwitched = false;
202

203
204
                staticDat_[eIdxGlobal].oldPhasePresence
                    = staticDat_[eIdxGlobal].phasePresence;  
205
            }
206
207
208
209
210
211
212
        }
    }

    /*!
     * \brief Compute the total storage inside one phase of all
     *        conservation quantities.
     *
213
     * \param storage Contains the storage of each component for one phase
214
215
     * \param phaseIdx The phase index
     */
216
    void globalPhaseStorage(PrimaryVariables &storage, const int phaseIdx)
217
    {
218
        storage = 0;
219

220
221
222
223
        ElementIterator eIt = this->gridView_().template begin<0>();
        const ElementIterator eEndIt = this->gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt) {
	 if(eIt->partitionType() == Dune::InteriorEntity)
224
           {
225
            this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
226

227
228
            for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
                storage += this->localResidual().storageTerm()[i];
229
230
           }
	}
231
        if (this->gridView_().comm().size() > 1)
232
            storage = this->gridView_().comm().sum(storage);
233
234
    }

235

236
237
    /*!
     * \brief Called by the update() method if applying the newton
238
     *        method was unsuccessful.
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
264
265
266
267
268
269
270
271
272
273
     */
    void updateFailed()
    {
        ParentType::updateFailed();

        setSwitched_(false);
        resetPhasePresence_();
    };

    /*!
     * \brief Called by the problem if a time integration was
     *        successful, post processing of the solution is done and the
     *        result has been written to disk.
     *
     * This should prepare the model for the next time integration.
     */
    void advanceTimeLevel()
    {
        ParentType::advanceTimeLevel();

        // update the phase state
        updateOldPhasePresence_();
        setSwitched_(false);
    }

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

    /*!
274
     * \brief Returns the phase presence of the current or the old solution of a degree of freedom.
275
     *
276
     * \param dofIdxGlobal The global index of the degree of freedom
277
278
     * \param oldSol Evaluate function with solution of current or previous time step
     */
279
    int phasePresence(int dofIdxGlobal, bool oldSol) const
280
281
282
    {
        return
            oldSol
283
284
            ? staticDat_[dofIdxGlobal].oldPhasePresence
            : staticDat_[dofIdxGlobal].phasePresence;
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    }

    /*!
     * \brief Append all quantities of interest which can be derived
     *        from the solution of the current time step to the VTK
     *        writer.
     *
     * \param sol The solution vector
     * \param writer The writer for multi-file VTK datasets
     */
    template<class MultiWriter>
    void addOutputVtkFields(const SolutionVector &sol,
                            MultiWriter &writer)
    {
        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
300
        typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
301

302
303
        // get the number of degrees of freedom
        unsigned numDofs = this->numDofs();
304

305
        // create the required scalar fields
306
307
308
309
310
        ScalarField *saturation[numPhases];
        ScalarField *pressure[numPhases];
        ScalarField *density[numPhases];

        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
311
312
313
            saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
            pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
            density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
314
315
        }

316
        ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
317
318
319
        ScalarField *moleFraction[numPhases][numComponents];
        for (int i = 0; i < numPhases; ++i)
            for (int j = 0; j < numComponents; ++j)
320
321
322
                moleFraction[i][j] = writer.allocateManagedBuffer (numDofs);
        ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
        ScalarField *poro = writer.allocateManagedBuffer(numDofs);
323
324
325
        VectorField *velocityN = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
        VectorField *velocityW = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
        VectorField *velocityG = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
326
        ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
327

328
329
330
331
332
333
334
335
336
337
        if (velocityOutput.enableOutput()) // check if velocity output is demanded
        {
            // initialize velocity fields
            for (unsigned int i = 0; i < numDofs; ++i)
            {
                (*velocityN)[i] = Scalar(0);
                (*velocityW)[i] = Scalar(0);
                (*velocityG)[i] = Scalar(0);
            }
        }
338

339
340
        unsigned numElements = this->gridView_().size(0);
        ScalarField *rank = writer.allocateManagedBuffer (numElements);
341

342
343
344
        ElementIterator eIt = this->gridView_().template begin<0>();
        ElementIterator eEndIt = this->gridView_().template end<0>();
        for (; eIt != eEndIt; ++eIt)
345
        {
346
347
            if(eIt->partitionType() == Dune::InteriorEntity)
            {
348
349
350
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                int eIdx = this->problem_().elementMapper().index(*eIt);
#else
351
                int eIdx = this->problem_().elementMapper().map(*eIt);
352
#endif
353
                (*rank)[eIdx] = this->gridView_().comm().rank();
354

355
356
                FVElementGeometry fvGeometry;
                fvGeometry.update(this->gridView_(), *eIt);
357

358

359
360
361
362
363
                ElementVolumeVariables elemVolVars;
                elemVolVars.update(this->problem_(),
                                   *eIt,
                                   fvGeometry,
                                   false /* oldSol? */);
364

365
366
                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
                {
367
368
369
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                    int dofIdxGlobal = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
#else
370
                    int dofIdxGlobal = this->dofMapper().map(*eIt, scvIdx, dofCodim);
371
#endif
372

373
                    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
374
375
376
                        (*saturation[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].saturation(phaseIdx);
                        (*pressure[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
                        (*density[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
377
                    }
378

379
380
                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
381
                            (*moleFraction[phaseIdx][compIdx])[dofIdxGlobal] =
382
383
                                elemVolVars[scvIdx].moleFraction(phaseIdx,
                                                                  compIdx);
384

385
                            Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[dofIdxGlobal]);
386
                        }
387
                    }
388

389
390
391
                    (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
                    (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
                    (*phasePresence)[dofIdxGlobal] = staticDat_[dofIdxGlobal].phasePresence;
392
393
                }

394
395
396
397
                // velocity output
                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
398
399
400
            }
        }

401
402
403
        writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
        writer.attachDofData(*saturation[nPhaseIdx], "Sn", isBox);
        writer.attachDofData(*saturation[gPhaseIdx], "Sg", isBox);
404
405
406
407
408
409
410
411
        writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
        writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
        writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
        writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
        writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
        writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);

        for (int i = 0; i < numPhases; ++i)
412
        {
413
            for (int j = 0; j < numComponents; ++j)
414
            {
415
416
417
418
419
420
                std::ostringstream oss;
                oss << "x^"
                    << FluidSystem::phaseName(i)
                    << "_"
                    << FluidSystem::componentName(j);
                writer.attachDofData(*moleFraction[i][j], oss.str().c_str(), isBox);
421
422
            }
        }
423
424
425
        writer.attachDofData(*poro, "porosity", isBox);
        writer.attachDofData(*temperature, "temperature", isBox);
        writer.attachDofData(*phasePresence, "phase presence", isBox);
426
427
428
429
430
431
432
433

        if (velocityOutput.enableOutput()) // check if velocity output is demanded
        {
            writer.attachDofData(*velocityW,  "velocityW", isBox, dim);
            writer.attachDofData(*velocityN,  "velocityN", isBox, dim);
            writer.attachDofData(*velocityG,  "velocityG", isBox, dim);
        }

434
435
436
437
438
439
        writer.attachCellData(*rank, "process rank");
    }

    /*!
     * \brief Write the current solution to a restart file.
     *
440
441
     * \param outStream The output stream of one entity for the restart file
     * \param entity The entity, either a vertex or an element
442
     */
443
444
    template<class Entity>
    void serializeEntity(std::ostream &outStream, const Entity &entity)
445
446
    {
        // write primary variables
447
        ParentType::serializeEntity(outStream, entity);
448

449
450
451
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
        int dofIdxGlobal = this->dofMapper().index(entity);
#else
452
        int dofIdxGlobal = this->dofMapper().map(entity);
453
#endif
454
        if (!outStream.good())
455
            DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal);
456

457
        outStream << staticDat_[dofIdxGlobal].phasePresence << " ";
458
459
460
    }

    /*!
461
     * \brief Reads the current solution from a restart file.
462
     *
463
464
     * \param inStream The input stream of one entity from the restart file
     * \param entity The entity, either a vertex or an element
465
     */
466
467
    template<class Entity>
    void deserializeEntity(std::istream &inStream, const Entity &entity)
468
469
    {
        // read primary variables
470
        ParentType::deserializeEntity(inStream, entity);
471
472

        // read phase presence
473
474
475
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
        int dofIdxGlobal = this->dofMapper().index(entity);
#else
476
        int dofIdxGlobal = this->dofMapper().map(entity);
477
#endif
478
479
        if (!inStream.good())
            DUNE_THROW(Dune::IOError,
480
                       "Could not deserialize entity " << dofIdxGlobal);
481

482
483
484
        inStream >> staticDat_[dofIdxGlobal].phasePresence;
        staticDat_[dofIdxGlobal].oldPhasePresence
            = staticDat_[dofIdxGlobal].phasePresence;
485
486
487
488
489
490
491
492
493
494
495
496
497

    }

    /*!
     * \brief Update the static data of all vertices in the grid.
     *
     * \param curGlobalSol The current global solution
     * \param oldGlobalSol The previous global solution
     */
    void updateStaticData(SolutionVector &curGlobalSol,
                          const SolutionVector &oldGlobalSol)
    {
        bool wasSwitched = false;
498
499
        int succeeded;
        try {
500

501
502
            for (unsigned i = 0; i < staticDat_.size(); ++i)
                staticDat_[i].visited = false;
503

504
505
506
507
508
            FVElementGeometry fvGeometry;
            static VolumeVariables volVars;
            ElementIterator eIt = this->gridView_().template begin<0> ();
            const ElementIterator &eEndIt = this->gridView_().template end<0> ();
            for (; eIt != eEndIt; ++eIt)
509
            {
510
511
512
                fvGeometry.update(this->gridView_(), *eIt);
                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
                {
513
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
514
                    int dofIdxGlobal = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
515
#else
516
                    int dofIdxGlobal = this->dofMapper().map(*eIt, scvIdx, dofCodim);
517
#endif
518

519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
                    if (staticDat_[dofIdxGlobal].visited)
                        continue;

                    staticDat_[dofIdxGlobal].visited = true;
                    volVars.update(curGlobalSol[dofIdxGlobal],
                            this->problem_(),
                            *eIt,
                            fvGeometry,
                            scvIdx,
                            false);
                    const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
                    if (primaryVarSwitch_(curGlobalSol,
                            volVars,
                            dofIdxGlobal,
                            globalPos))
                    {
                        this->jacobianAssembler().markDofRed(dofIdxGlobal);
                        wasSwitched = true;
                    }
538
                }
539
            }
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
            succeeded = 1;
        }
        catch (Dumux::NumericalProblem &e)
        {
            std::cout << "\n"
                      << "Rank " << this->problem_().gridView().comm().rank()
                      << " caught an exception while updating the static data." << e.what()
                      << "\n";
            succeeded = 0;
        }
        //make sure that all processes succeeded. If not throw a NumericalProblem to decrease the time step size.
        if (this->gridView_().comm().size() > 1)
            succeeded = this->gridView_().comm().min(succeeded);

        if (!succeeded) {
            if(this->problem_().gridView().comm().rank() == 0)
                DUNE_THROW(NumericalProblem,
                        "A process did not succeed in updating the static data.");
            return;
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
        }

        // make sure that if there was a variable switch in an
        // other partition we will also set the switch flag
        // for our partition.
        if (this->gridView_().comm().size() > 1)
            wasSwitched = this->gridView_().comm().max(wasSwitched);

        setSwitched_(wasSwitched);
    }

protected:
    /*!
     * \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_()
    {
591
        for (unsigned int i = 0; i < this->numDofs(); ++i)
592
        {
593
594
595
            staticDat_[i].phasePresence
                = staticDat_[i].oldPhasePresence;
            staticDat_[i].wasSwitched = false;
596
597
598
599
600
601
602
603
        }
    }

    /*!
     * \brief Set the old phase of all verts state to the current one.
     */
    void updateOldPhasePresence_()
    {
604
        for (unsigned int i = 0; i < this->numDofs(); ++i)
605
        {
606
607
608
            staticDat_[i].oldPhasePresence
                = staticDat_[i].phasePresence;
            staticDat_[i].wasSwitched = false;
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
        }
    }

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

    //  perform variable switch at a vertex; Returns true if a
    //  variable switch was performed.
    bool primaryVarSwitch_(SolutionVector &globalSol,
                           const VolumeVariables &volVars,
625
                           int dofIdxGlobal,
626
627
628
629
                           const GlobalPosition &globalPos)
    {
        // evaluate primary variable switch
        bool wouldSwitch = false;
630
        int phasePresence = staticDat_[dofIdxGlobal].phasePresence;
631
632
633
634
635
636
        int newPhasePresence = phasePresence;

        // check if a primary var switch is necessary
        if (phasePresence == threePhases)
        {
            Scalar Smin = 0;
637
            if (staticDat_[dofIdxGlobal].wasSwitched)
638
639
640
641
642
643
                Smin = -0.01;

            if (volVars.saturation(gPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // gas phase disappears
644
                std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
645
                          << ", coordinates: " << globalPos << ", sg: "
646
647
648
                          << volVars.saturation(gPhaseIdx) << std::endl;
                newPhasePresence = wnPhaseOnly;

649
                globalSol[dofIdxGlobal][switch1Idx]
650
                    = volVars.moleFraction(wPhaseIdx, gCompIdx);
651
652
653
654
655
            }
            else if (volVars.saturation(wPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // water phase disappears
656
                std::cout << "Water phase disappears at vertex " << dofIdxGlobal
657
                          << ", coordinates: " << globalPos << ", sw: "
658
659
660
                          << volVars.saturation(wPhaseIdx) << std::endl;
                newPhasePresence = gnPhaseOnly;

661
                globalSol[dofIdxGlobal][switch1Idx]
662
                    = volVars.moleFraction(gPhaseIdx, wCompIdx);
663
664
665
666
667
            }
            else if (volVars.saturation(nPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // NAPL phase disappears
668
                std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
669
                          << ", coordinates: " << globalPos << ", sn: "
670
671
672
                          << volVars.saturation(nPhaseIdx) << std::endl;
                newPhasePresence = wgPhaseOnly;

673
                globalSol[dofIdxGlobal][switch2Idx]
674
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
675
676
677
678
679
680
681
682
            }
        }
        else if (phasePresence == wPhaseOnly)
        {
            bool gasFlag = 0;
            bool nonwettingFlag = 0;
            // calculate fractions of the partial pressures in the
            // hypothetical gas phase
683
684
685
            Scalar xwg = volVars.moleFraction(gPhaseIdx, wCompIdx);
            Scalar xgg = volVars.moleFraction(gPhaseIdx, gCompIdx);
            Scalar xng = volVars.moleFraction(gPhaseIdx, nCompIdx);
686
687
688
689
690
691
692
693
694
            /* take care:
               for xgg in case wPhaseOnly we compute xgg=henry_air*x2w
               for xwg in case wPhaseOnly we compute xwg=pwsat
               for xng in case wPhaseOnly we compute xng=henry_NAPL*x1w
            */

            Scalar xgMax = 1.0;
            if (xwg + xgg + xng > xgMax)
                wouldSwitch = true;
695
            if (staticDat_[dofIdxGlobal].wasSwitched)
696
697
698
699
700
701
702
                xgMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, gas phase appears
            if (xwg + xgg + xng > xgMax)
            {
                // gas phase appears
703
                std::cout << "gas phase appears at vertex " << dofIdxGlobal
704
705
706
707
708
709
                          << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
                          << xwg + xgg + xng << std::endl;
                gasFlag = 1;
            }

            // calculate fractions in the hypothetical NAPL phase
710
            Scalar xnn = volVars.moleFraction(nPhaseIdx, nCompIdx);
711
712
713
714
715
716
717
718
719
720
            /* take care:
               for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
               where a hypothetical gas pressure is assumed for the Henry
               x0n is set to NULL  (all NAPL phase is dirty)
               x2n is set to NULL  (all NAPL phase is dirty)
            */

            Scalar xnMax = 1.0;
            if (xnn > xnMax)
                wouldSwitch = true;
721
            if (staticDat_[dofIdxGlobal].wasSwitched)
722
723
724
725
726
727
728
                xnMax *= 1.02;

            // if the sum of the hypothetical mole fractions would be larger than
            // 100%, NAPL phase appears
            if (xnn > xnMax)
            {
                // NAPL phase appears
729
                std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
730
731
732
733
734
735
736
737
                          << ", coordinates: " << globalPos << ", xnn: "
                          << xnn << std::endl;
                nonwettingFlag = 1;
            }

            if ((gasFlag == 1) && (nonwettingFlag == 0))
            {
                newPhasePresence = wgPhaseOnly;
738
739
                globalSol[dofIdxGlobal][switch1Idx] = 0.9999;
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
740
741
742
743
            }
            else if ((gasFlag == 1) && (nonwettingFlag == 1))
            {
                newPhasePresence = threePhases;
744
745
                globalSol[dofIdxGlobal][switch1Idx] = 0.9999;
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
746
747
748
749
            }
            else if ((gasFlag == 0) && (nonwettingFlag == 1))
            {
                newPhasePresence = wnPhaseOnly;
750
                globalSol[dofIdxGlobal][switch1Idx]
751
                    = volVars.moleFraction(wPhaseIdx, gCompIdx);
752
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
753
754
755
756
757
758
759
760
            }
        }
        else if (phasePresence == gnPhaseOnly)
        {
            bool nonwettingFlag = 0;
            bool wettingFlag = 0;

            Scalar Smin = 0.0;
761
            if (staticDat_[dofIdxGlobal].wasSwitched)
762
763
764
765
766
767
                Smin = -0.01;

            if (volVars.saturation(nPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // NAPL phase disappears
768
                std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
769
                          << ", coordinates: " << globalPos << ", sn: "
770
771
772
773
774
775
                          << volVars.saturation(nPhaseIdx) << std::endl;
                nonwettingFlag = 1;
            }


            // calculate fractions of the hypothetical water phase
776
            Scalar xww = volVars.moleFraction(wPhaseIdx, wCompIdx);
777
778
779
780
781
782
783
            /*
              take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
              If this is larger than 1, then water appears
            */
            Scalar xwMax = 1.0;
            if (xww > xwMax)
                wouldSwitch = true;
784
            if (staticDat_[dofIdxGlobal].wasSwitched)
785
786
787
788
789
790
791
                xwMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, gas phase appears
            if (xww > xwMax)
            {
                // water phase appears
792
                std::cout << "water phase appears at vertex " << dofIdxGlobal
793
794
795
796
797
798
799
800
                          << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
                          << xww << std::endl;
                wettingFlag = 1;
            }

            if ((wettingFlag == 1) && (nonwettingFlag == 0))
            {
                newPhasePresence = threePhases;
801
802
                globalSol[dofIdxGlobal][switch1Idx] = 0.0001;
                globalSol[dofIdxGlobal][switch2Idx] = volVars.saturation(nPhaseIdx);
803
804
805
806
            }
            else if ((wettingFlag == 1) && (nonwettingFlag == 1))
            {
                newPhasePresence = wgPhaseOnly;
807
808
                globalSol[dofIdxGlobal][switch1Idx] = 0.0001;
                globalSol[dofIdxGlobal][switch2Idx]
809
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
810
811
812
813
            }
            else if ((wettingFlag == 0) && (nonwettingFlag == 1))
            {
                newPhasePresence = gPhaseOnly;
814
                globalSol[dofIdxGlobal][switch1Idx]
815
                    = volVars.moleFraction(gPhaseIdx, wCompIdx);
816
                globalSol[dofIdxGlobal][switch2Idx]
817
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
818
819
820
821
822
823
824
825
            }
        }
        else if (phasePresence == wnPhaseOnly)
        {
            bool nonwettingFlag = 0;
            bool gasFlag = 0;

            Scalar Smin = 0.0;
826
            if (staticDat_[dofIdxGlobal].wasSwitched)
827
828
829
830
831
832
                Smin = -0.01;

            if (volVars.saturation(nPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // NAPL phase disappears
833
                std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
834
                          << ", coordinates: " << globalPos << ", sn: "
835
836
837
838
839
840
                          << volVars.saturation(nPhaseIdx) << std::endl;
                nonwettingFlag = 1;
            }

            // calculate fractions of the partial pressures in the
            // hypothetical gas phase
841
842
843
            Scalar xwg = volVars.moleFraction(gPhaseIdx, wCompIdx);
            Scalar xgg = volVars.moleFraction(gPhaseIdx, gCompIdx);
            Scalar xng = volVars.moleFraction(gPhaseIdx, nCompIdx);
844
845
846
847
848
849
850
851
            /* take care:
               for xgg in case wPhaseOnly we compute xgg=henry_air*x2w
               for xwg in case wPhaseOnly we compute xwg=pwsat
               for xng in case wPhaseOnly we compute xng=henry_NAPL*x1w
            */
            Scalar xgMax = 1.0;
            if (xwg + xgg + xng > xgMax)
                wouldSwitch = true;
852
            if (staticDat_[dofIdxGlobal].wasSwitched)
853
854
855
856
857
858
859
                xgMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, gas phase appears
            if (xwg + xgg + xng > xgMax)
            {
                // gas phase appears
860
                std::cout << "gas phase appears at vertex " << dofIdxGlobal
861
862
863
864
865
866
867
868
                          << ", coordinates: " << globalPos << ", xwg + xgg + xng: "
                          << xwg + xgg + xng << std::endl;
                gasFlag = 1;
            }

            if ((gasFlag == 1) && (nonwettingFlag == 0))
            {
                newPhasePresence = threePhases;
869
870
                globalSol[dofIdxGlobal][switch1Idx] = volVars.saturation(wPhaseIdx);
                globalSol[dofIdxGlobal][switch2Idx] = volVars.saturation(nPhaseIdx);
871
872
873
874
            }
            else if ((gasFlag == 1) && (nonwettingFlag == 1))
            {
                newPhasePresence = wgPhaseOnly;
875
876
                globalSol[dofIdxGlobal][switch1Idx] = volVars.saturation(wPhaseIdx);
                globalSol[dofIdxGlobal][switch2Idx]
877
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
878
879
880
881
            }
            else if ((gasFlag == 0) && (nonwettingFlag == 1))
            {
                newPhasePresence = wPhaseOnly;
882
                globalSol[dofIdxGlobal][switch1Idx]
883
                    = volVars.moleFraction(wPhaseIdx, gCompIdx);
884
                globalSol[dofIdxGlobal][switch2Idx]
885
                    = volVars.moleFraction(wPhaseIdx, nCompIdx);
886
887
888
889
890
891
892
893
            }
        }
        else if (phasePresence == gPhaseOnly)
        {
            bool nonwettingFlag = 0;
            bool wettingFlag = 0;

            // calculate fractions in the hypothetical NAPL phase
894
            Scalar xnn = volVars.moleFraction(nPhaseIdx, nCompIdx);
895
896
897
898
899
900
901
902
            /*
              take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
              if this is larger than 1, then NAPL appears
            */

            Scalar xnMax = 1.0;
            if (xnn > xnMax)
                wouldSwitch = true;
903
            if (staticDat_[dofIdxGlobal].wasSwitched)
904
905
906
907
908
909
910
                xnMax *= 1.02;

            // if the sum of the hypothetical mole fraction would be larger than
            // 100%, NAPL phase appears
            if (xnn > xnMax)
            {
                // NAPL phase appears
911
                std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
912
913
914
915
916
                          << ", coordinates: " << globalPos << ", xnn: "
                          << xnn << std::endl;
                nonwettingFlag = 1;
            }
            // calculate fractions of the hypothetical water phase
917
            Scalar xww = volVars.moleFraction(wPhaseIdx, wCompIdx);
918
919
920
921
922
923
924
            /*
              take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
              If this is larger than 1, then water appears
            */
            Scalar xwMax = 1.0;
            if (xww > xwMax)
                wouldSwitch = true;
925
            if (staticDat_[dofIdxGlobal].wasSwitched)
926
927
928
929
930
931
932
                xwMax *= 1.02;

            // if the sum of the mole fractions would be larger than
            // 100%, gas phase appears
            if (xww > xwMax)
            {
                // water phase appears
933
                std::cout << "water phase appears at vertex " << dofIdxGlobal
934
935
936
937
938
939
940
                          << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
                          << xww << std::endl;
                wettingFlag = 1;
            }
            if ((wettingFlag == 1) && (nonwettingFlag == 0))
            {
                newPhasePresence = wgPhaseOnly;
941
942
                globalSol[dofIdxGlobal][switch1Idx] = 0.0001;
                globalSol[dofIdxGlobal][switch2Idx]
943
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
944
945
946
947
            }
            else if ((wettingFlag == 1) && (nonwettingFlag == 1))
            {
                newPhasePresence = threePhases;
948
949
                globalSol[dofIdxGlobal][switch1Idx] = 0.0001;
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
950
951
952
953
            }
            else if ((wettingFlag == 0) && (nonwettingFlag == 1))
            {
                newPhasePresence = gnPhaseOnly;
954
                globalSol[dofIdxGlobal][switch1Idx]
955
                    = volVars.moleFraction(gPhaseIdx, wCompIdx);
956
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
957
958
959
960
961
962
963
964
965
            }
        }
        else if (phasePresence == wgPhaseOnly)
        {
            bool nonwettingFlag = 0;
            bool gasFlag = 0;
            bool wettingFlag = 0;

            // get the fractions in the hypothetical NAPL phase
966
            Scalar xnn = volVars.moleFraction(nPhaseIdx, nCompIdx);
967
968
969
970
971
972
973

            // take care: if the NAPL phase is not present, take
            // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
            // appears
            Scalar xnMax = 1.0;
            if (xnn > xnMax)
                wouldSwitch = true;
974
            if (staticDat_[dofIdxGlobal].wasSwitched)
975
976
977
978
979
980
981
                xnMax *= 1.02;

            // if the sum of the hypothetical mole fraction would be larger than
            // 100%, NAPL phase appears
            if (xnn > xnMax)
            {
                // NAPL phase appears
982
                std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
983
984
985
986
987
988
                          << ", coordinates: " << globalPos << ", xnn: "
                          << xnn << std::endl;
                nonwettingFlag = 1;
            }

            Scalar Smin = -1.e-6;
989
            if (staticDat_[dofIdxGlobal].wasSwitched)
990
991
992
993
994
995
                Smin = -0.01;

            if (volVars.saturation(gPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // gas phase disappears
996
                std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
997
                          << ", coordinates: " << globalPos << ", sg: "
998
999
1000
1001
1002
                          << volVars.saturation(gPhaseIdx) << std::endl;
                gasFlag = 1;
            }

            Smin = 0.0;
1003
            if (staticDat_[dofIdxGlobal].wasSwitched)
1004
1005
1006
1007
1008
1009
                Smin = -0.01;

            if (volVars.saturation(wPhaseIdx) <= Smin)
            {
                wouldSwitch = true;
                // gas phase disappears
1010
                std::cout << "Water phase disappears at vertex " << dofIdxGlobal
1011
                          << ", coordinates: " << globalPos << ", sw: "
1012
1013
1014
1015
1016
1017
1018
                          << volVars.saturation(wPhaseIdx) << std::endl;
                wettingFlag = 1;
            }

            if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
            {
                newPhasePresence = gnPhaseOnly;
1019
                globalSol[dofIdxGlobal][switch1Idx]
1020
                    = volVars.moleFraction(gPhaseIdx, wCompIdx);
1021
                globalSol[dofIdxGlobal][switch2Idx] = 0.0001;
1022
1023
1024
1025
            }
            else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
            {
                newPhasePresence = threePhases;
1026
1027
                globalSol[dofIdxGlobal][switch1Idx] = volVars.saturation(wPhaseIdx);
                globalSol[dofIdxGlobal][switch2Idx] = 0.0;
1028
1029
1030
1031
            }
            else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
            {
                newPhasePresence = wPhaseOnly;
1032
                globalSol[dofIdxGlobal][switch1Idx]
1033
                    = volVars.moleFraction(wPhaseIdx, gCompIdx);
1034
                globalSol[dofIdxGlobal][switch2Idx]
1035
                    = volVars.moleFraction(wPhaseIdx, nCompIdx);
1036
1037
1038
1039
            }
            else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
            {
                newPhasePresence = gPhaseOnly;
1040
                globalSol[dofIdxGlobal][switch1Idx]
1041
                    = volVars.moleFraction(gPhaseIdx, wCompIdx);
1042
                globalSol[dofIdxGlobal][switch2Idx]
1043
                    = volVars.moleFraction(gPhaseIdx, nCompIdx);
1044
1045
1046
            }
        }

1047
1048
        staticDat_[dofIdxGlobal].phasePresence = newPhasePresence;
        staticDat_[dofIdxGlobal].wasSwitched = wouldSwitch;
1049
1050
1051
1052
        return phasePresence != newPhasePresence;
    }

    // parameters given in constructor
1053
    std::vector<StaticVars> staticDat_;
1054
1055
1056
1057
1058
1059
1060
1061
    bool switchFlag_;
};

}

#include "3p3cpropertydefaults.hh"

#endif