volumevariables.hh 23.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// -*- 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
 *
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
#ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH
#define DUMUX_2P2C_VOLUME_VARIABLES_HH

28
#include <dumux/material/fluidstates/compositional.hh>
29
30
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
Timo Koch's avatar
Timo Koch committed
31
#include <dumux/porousmediumflow/volumevariables.hh>
32
#include <dumux/discretization/methods.hh>
33
34
35
36
37
38
39
40
41
42

namespace Dumux
{

/*!
 * \ingroup TwoPTwoCModel
 * \brief Contains the quantities which are constant within a
 *        finite volume in the two-phase two-component model.
 */
template <class TypeTag>
Timo Koch's avatar
Timo Koch committed
43
class TwoPTwoCVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
44
{
Timo Koch's avatar
Timo Koch committed
45
    using ParentType = PorousMediumFlowVolumeVariables<TypeTag>;
46
    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
47
48

    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
49
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
50
51
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
52
    using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
53
54
55
56
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
    using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
57

58
59
60
    // component indices
    enum
    {
61
62
63
64
65
66
67
        wCompIdx = Indices::wCompIdx,
        nCompIdx = Indices::nCompIdx,
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx
    };

    // present phases
68
69
    enum
    {
70
71
72
73
74
75
        wPhaseOnly = Indices::wPhaseOnly,
        nPhaseOnly = Indices::nPhaseOnly,
        bothPhases = Indices::bothPhases
    };

    // formulations
76
77
    enum
    {
78
79
80
81
82
83
        formulation = GET_PROP_VALUE(TypeTag, Formulation),
        pwsn = TwoPTwoCFormulation::pwsn,
        pnsw = TwoPTwoCFormulation::pnsw
    };

    // primary variable indices
84
85
    enum
    {
86
87
88
89
        switchIdx = Indices::switchIdx,
        pressureIdx = Indices::pressureIdx
    };

90
    using Element = typename GridView::template Codim<0>::Entity;
91
92
93
94
95
96
    using PermeabilityType = typename SpatialParams::PermeabilityType;
    using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;

    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
    static constexpr bool useKelvinEquation = GET_PROP_VALUE(TypeTag, UseKelvinEquation);
    static constexpr bool useConstraintSolver = GET_PROP_VALUE(TypeTag, UseConstraintSolver);
97
98
99
    static_assert(useMoles || (!useMoles && useConstraintSolver),
                  "if UseMoles is set false, UseConstraintSolver has to be set to true");

100
101
    using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem, useKelvinEquation>;

102
103
104
    static constexpr int dim = GridView::dimension;
    static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
    static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
105
public:
Timo Koch's avatar
Timo Koch committed
106
    //! The type of the object returned by the fluidState() method
107
    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
108

Timo Koch's avatar
Timo Koch committed
109
110
111
112
    //! Pull member functions of the parent for deriving classes
    using ParentType::temperature;
    using ParentType::enthalpy;

113
    /*!
Timo Koch's avatar
Timo Koch committed
114
     * \copydoc PorousMediumFlowVolumeVariables::update
115
     */
116
    void update(const ElementSolution &elemSol,
117
118
                const Problem &problem,
                const Element &element,
119
                const SubControlVolume& scv)
120
    {
121
        ParentType::update(elemSol, problem, element, scv);
122

Timo Koch's avatar
Timo Koch committed
123
        Implementation::completeFluidState(elemSol, problem, element, scv, fluidState_);
124
125
126
127

        /////////////
        // calculate the remaining quantities
        /////////////
128
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
129
130
131
132
133
134

        // Second instance of a parameter cache.
        // Could be avoided if diffusion coefficients also
        // became part of the fluid state.
        typename FluidSystem::ParameterCache paramCache;
        paramCache.updateAll(fluidState_);
135
136
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
137
138
139
140
141
142
143
144
145
146
147
            // relative permeabilities
            Scalar kr;
            if (phaseIdx == wPhaseIdx)
                kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
            else // ATTENTION: krn requires the wetting phase saturation
                // as parameter!
                kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
            relativePermeability_[phaseIdx] = kr;
            Valgrind::CheckDefined(relativePermeability_[phaseIdx]);

            // binary diffusion coefficients
148
149
150
151
152
            diffCoeff_[phaseIdx] = FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                                           paramCache,
                                                                           phaseIdx,
                                                                           wCompIdx,
                                                                           nCompIdx);
153
154
        }

155
156
157
        // porosity & permeabilty
        porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
        permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
158
159
160
161
162
    }

    /*!
     * \copydoc ImplicitModel::completeFluidState
     */
163
    static void completeFluidState(const ElementSolution& elemSol,
164
165
                                   const Problem& problem,
                                   const Element& element,
166
167
                                   const SubControlVolume& scv,
                                   FluidState& fluidState)
168
    {
169
        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
170
        fluidState.setTemperature(t);
171

172
        const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
173
        const auto phasePresence = priVars.state();
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199

        /////////////
        // set the saturations
        /////////////
        Scalar sn;
        if (phasePresence == nPhaseOnly)
            sn = 1.0;
        else if (phasePresence == wPhaseOnly) {
            sn = 0.0;
        }
        else if (phasePresence == bothPhases) {
            if (formulation == pwsn)
                sn = priVars[switchIdx];
            else if (formulation == pnsw)
                sn = 1.0 - priVars[switchIdx];
            else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
        }
        else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
        fluidState.setSaturation(wPhaseIdx, 1 - sn);
        fluidState.setSaturation(nPhaseIdx, sn);

        /////////////
        // set the pressures of the fluid phases
        /////////////

        // calculate capillary pressure
200
        const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
        Scalar pc = MaterialLaw::pc(materialParams, 1 - sn);

        if (formulation == pwsn) {
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
        }
        else if (formulation == pnsw) {
            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
        }
        else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");

        /////////////
        // calculate the phase compositions
        /////////////
        typename FluidSystem::ParameterCache paramCache;

        //get the phase pressures and set the fugacity coefficients here if constraintsolver is not used
        Scalar pn = 0;
        Scalar pw = 0;

        if(!useConstraintSolver) {
            if (formulation == pwsn) {
                pw = priVars[pressureIdx];
                pn = pw + pc;
            }
            else {
                pn = priVars[pressureIdx];
                pw = pn - pc;
            }

            for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                assert(FluidSystem::isIdealMixture(phaseIdx));

                for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
                    Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
                    fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
                }
            }
        }

        // now comes the tricky part: calculate phase compositions
        if (phasePresence == bothPhases) {
            // both phases are present, phase compositions are a
245
            // result of the nonwetting <-> wetting equilibrium. This is
246
247
248
249
250
251
            // the job of the "MiscibleMultiPhaseComposition"
            // constraint solver
            if(useConstraintSolver) {
                MiscibleMultiPhaseComposition::solve(fluidState,
                                                     paramCache,
                                                     /*setViscosity=*/true,
252
                                                     /*setEnthalpy=*/false);
253
254
255
256
257
258
259
260
261
            }
            // ... or calculated explicitly this way ...
            else {
                //get the partial pressure of the main component of the the wetting phase ("H20") within the nonwetting (gas) phase == vapor pressure due to equilibrium
                //note that in this case the fugacityCoefficient * pw is the vapor pressure (see implementation in respective fluidsystem)
                Scalar partPressH2O = FluidSystem::fugacityCoefficient(fluidState,
                                                                       wPhaseIdx,
                                                                       wCompIdx) * pw;

262
263
264
                if (useKelvinEquation)
                    partPressH2O = FluidSystem::kelvinVaporPressure(fluidState, wPhaseIdx, wCompIdx);

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
                // get the partial pressure of the main component of the the nonwetting (gas) phase ("Air")
                Scalar partPressAir = pn - partPressH2O;

                //calculate the mole fractions of the components within the nonwetting phase
                Scalar xnn = partPressAir/pn;
                Scalar xnw = partPressH2O/pn;

                // calculate the mole fractions of the components within the wetting phase
                //note that in this case the fugacityCoefficient * pw is the Henry Coefficient (see implementation in respective fluidsystem)
                Scalar xwn = partPressAir
                  / (FluidSystem::fugacityCoefficient(fluidState,
                                                      wPhaseIdx,nCompIdx)
                  * pw);

                Scalar xww = 1.0 -xwn;

                //set all mole fractions
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xww);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
            }
        }
288
289
        else if (phasePresence == nPhaseOnly)
        {
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
            // only the nonwetting phase is present, i.e. nonwetting phase
            // composition is stored explicitly.
            if(useMoles)
            {
                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, 1 - priVars[switchIdx]);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
            }
            else
            {
                // setMassFraction() has only to be called 1-numComponents times
                fluidState.setMassFraction(nPhaseIdx, wCompIdx, priVars[switchIdx]);
            }

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 nPhaseIdx,
                                                 /*setViscosity=*/true,
311
                                                 /*setEnthalpy=*/false);
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
            }
            // ... or calculated explicitly this way ...
            else {
                // note that the water phase is actually not existing!
                // thus, this is used as phase switch criterion
                Scalar xnw = priVars[switchIdx];
                Scalar xnn = 1.0 -xnw;

                //first, xww:
                // xnw * pn = "actual" (hypothetical) vapor pressure
                // fugacityCoefficient * pw = vapor pressure given by thermodynamic conditions
                // Here, xww is not actually the mole fraction of water in the wetting phase
                // xww is only the ratio of "actual" vapor pressure / "thermodynamic" vapor pressure
                // If xww > 1 : gas is over-saturated with water vapor,
                // condensation takes place (see switch criterion in model)
                Scalar xww = xnw * pn
                  / (FluidSystem::fugacityCoefficient(fluidState,
                                                      wPhaseIdx,wCompIdx)
                     * pw);

                // now, xwn:
                //partialPressure / xwn = Henry
                //partialPressure = xnn * pn
                //xwn = xnn * pn / Henry
                // Henry = fugacityCoefficient * pw
                Scalar xwn = xnn * pn / (FluidSystem::fugacityCoefficient(fluidState,
                                                                          wPhaseIdx,nCompIdx)
                                         * pw);

                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, xww);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
            }
        }
345
346
        else if (phasePresence == wPhaseOnly)
        {
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
            // only the wetting phase is present, i.e. wetting phase
            // composition is stored explicitly.
            if(useMoles) // mole-fraction formulation
            {
                fluidState.setMoleFraction(wPhaseIdx, wCompIdx, 1-priVars[switchIdx]);
                fluidState.setMoleFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
            }
            else // mass-fraction formulation
            {
                // setMassFraction() has only to be called 1-numComponents times
                fluidState.setMassFraction(wPhaseIdx, nCompIdx, priVars[switchIdx]);
            }

            // calculate the composition of the remaining phases (as
            // well as the densities of all phases). This is the job
            // of the "ComputeFromReferencePhase" constraint solver
            if (useConstraintSolver) {
                ComputeFromReferencePhase::solve(fluidState,
                                                 paramCache,
                                                 wPhaseIdx,
                                                 /*setViscosity=*/true,
368
                                                 /*setEnthalpy=*/false);
369
370
            }
            // ... or calculated explicitly this way ...
371
372
            else
            {
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
                // note that the gas phase is actually not existing!
                // thus, this is used as phase switch criterion
                Scalar xwn = priVars[switchIdx];

                //first, xnw:
                //psteam = xnw * pn = partial pressure of water in gas phase
                //psteam = fugacityCoefficient * pw
                Scalar xnw = (FluidSystem::fugacityCoefficient(fluidState,
                                                               wPhaseIdx,wCompIdx)
                              * pw) / pn ;

                //now, xnn:
                // xwn = partialPressure / Henry
                // partialPressure = pn * xnn
                // xwn = pn * xnn / Henry
                // xnn = xwn * Henry / pn
                // Henry = fugacityCoefficient * pw
                Scalar xnn = xwn * (FluidSystem::fugacityCoefficient(fluidState,
                                                                     wPhaseIdx,nCompIdx)
                                    * pw) / pn ;

                fluidState.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
                fluidState.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
            }
        }

399
400
        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
        {
Timo Koch's avatar
Timo Koch committed
401
            // set the viscosity and desity here if constraintsolver is not used
402
403
            if(!useConstraintSolver)
            {
Timo Koch's avatar
Timo Koch committed
404
405
406
                paramCache.updateComposition(fluidState, phaseIdx);
                const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
                fluidState.setDensity(phaseIdx, rho);
407
                const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
408
409
410
                fluidState.setViscosity(phaseIdx,mu);
            }
            // compute and set the enthalpy
411
            Scalar h = ParentType::enthalpy(fluidState, paramCache, phaseIdx);
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
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
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
            fluidState.setEnthalpy(phaseIdx, h);
        }
   }


    /*!
     * \brief Returns the phase state within the control volume.
     */
    const FluidState &fluidState() const
    { return fluidState_; }

    /*!
     * \brief Returns the saturation of a given phase within
     *        the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar saturation(const int phaseIdx) const
    { return fluidState_.saturation(phaseIdx); }

    /*!
     * \brief Returns the mass fraction of a given component in a
     *        given phase within the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar massFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.massFraction(phaseIdx, compIdx); }

    /*!
     * \brief Returns the mole fraction of a given component in a
     *        given phase within the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     * \param compIdx The component index
     */
    Scalar moleFraction(const int phaseIdx, const int compIdx) const
    { return fluidState_.moleFraction(phaseIdx, compIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[kg/m^3]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar density(const int phaseIdx) const
    { return fluidState_.density(phaseIdx); }

    /*!
     * \brief Returns the dynamic viscosity of the fluid within the
     *        control volume in \f$\mathrm{[Pa s]}\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar viscosity(const int phaseIdx) const
    { return fluidState_.viscosity(phaseIdx); }

    /*!
     * \brief Returns the mass density of a given phase within the
     *        control volume in \f$[mol/m^3]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar molarDensity(const int phaseIdx) const
    { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }

    /*!
     * \brief Returns the effective pressure of a given phase within
     *        the control volume in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar pressure(const int phaseIdx) const
    { return fluidState_.pressure(phaseIdx); }

    /*!
     * \brief Returns temperature within the control volume in \f$[K]\f$.
     *
     * Note that we assume thermodynamic equilibrium, i.e. the
     * temperature of the rock matrix and of all fluid phases are
     * identical.
     */
    Scalar temperature() const
    { return fluidState_.temperature(/*phaseIdx=*/0); }

    /*!
     * \brief Returns the relative permeability of a given phase within
     *        the control volume in \f$[-]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar relativePermeability(const int phaseIdx) const
    {
        return relativePermeability_[phaseIdx];
    }

    /*!
     * \brief Returns the effective mobility of a given phase within
     *        the control volume in \f$[s*m/kg]\f$.
     *
     * \param phaseIdx The phase index
     */
    Scalar mobility(const int phaseIdx) const
    {
        return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
    }

    /*!
     * \brief Returns the effective capillary pressure within the control volume
     *        in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
     */
    Scalar capillaryPressure() const
    { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }

    /*!
     * \brief Returns the average porosity within the control volume in \f$[-]\f$.
     */
    Scalar porosity() const
    { return porosity_; }

533
534
535
    /*!
     * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
     */
536
    const PermeabilityType& permeability() const
537
538
    { return permeability_; }

539
540
541
    /*!
     * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
     */
542
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
543
544
545
546
547
548
    {
        if(phaseIdx == compIdx)
            DUNE_THROW(Dune::InvalidStateException, "Diffusion coeffiecient called for phaseIdx = compIdx");
        else
            return diffCoeff_[phaseIdx];
    }
549
550
551
552
553


protected:

    Scalar porosity_; //!< Effective porosity within the control volume
554
    PermeabilityType permeability_; //!< Effective permeability within the control volume
555
556
557
558
559
560
561
562
563
564
565
566
    Scalar relativePermeability_[numPhases]; //!< Relative permeability within the control volume
    Scalar diffCoeff_[numPhases]; //!< Binary diffusion coefficients for the phases
    FluidState fluidState_;

private:
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }

    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }
};

567
} // end namespace Dumux
568
569

#endif