buckleyleverettanalyticsolution.hh 11.8 KB
Newer Older
1
2
3
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
4
5
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
6
7
8
9
10
11
12
 *   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          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
18
19
20
21
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH
#define DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH

22
#include <dumux/common/properties.hh>
23
#include <dumux/porousmediumflow/2p/sequential/properties.hh>
24
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
25
26
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>

27
28
29
30
/*!
 * \file
 * \brief  Analytical solution of the buckley-leverett problem
 * \author Markus Wolff
31
32
 */

Kai Wendel's avatar
Kai Wendel committed
33
namespace Dumux {
34
35
/*!
 * \brief IMplicit Pressure Explicit Saturation (IMPES) scheme for the solution of
36
37
38
 * the Buckley-Leverett problem
 */

39
40
template<typename Scalar, typename Law>
struct CheckMaterialLaw
41
42
43
44
45
46
47
{
    static bool isLinear()
    {
        return false;
    }
};

48
49
50
51
52
53
54
55
56
57
58
template<typename Scalar>
struct CheckMaterialLaw<Scalar, LinearMaterial<Scalar> >
{
    static bool isLinear()
    {
        return true;
    }
};

template<typename Scalar>
struct CheckMaterialLaw<Scalar, EffToAbsLaw< LinearMaterial<Scalar> > >
59
60
61
62
63
64
65
66
67
{
    static bool isLinear()
    {
        return true;
    }
};

template<class TypeTag> class BuckleyLeverettAnalytic
{
68
69
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
70
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
71
72
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
Kai Wendel's avatar
Kai Wendel committed
73
74
    using MaterialLaw = typename SpatialParams::MaterialLaw;
    using MaterialLawParams = typename MaterialLaw::Params;
75
76
77
78
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using FluidState = GetPropType<TypeTag, Properties::FluidState>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using CellData = GetPropType<TypeTag, Properties::CellData>;
Kai Wendel's avatar
Kai Wendel committed
79
80
81
82
83
84
85
86
87
    static constexpr int dim = GridView::dimension;
    static constexpr int dimworld = GridView::dimensionworld;
    static constexpr int wPhaseIdx = Indices::wPhaseIdx;
    static constexpr int nPhaseIdx = Indices::nPhaseIdx;
    static constexpr int satEqIdx = Indices::satEqIdx;
    using BlockVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1> >;
    using Element = typename GridView::Traits::template Codim<0>::Entity;
    using ElementIterator = typename GridView::template Codim<0>::Iterator;
    using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
88
89
90
91
92
93
94
95
96
97
98
99
100
101

public:

    // functions needed for analytical solution

    void initializeAnalytic()
    {
        analyticSolution_.resize(size_);
        analyticSolution_ = 0;
        error_.resize(size_);
        error_ = 0;
        elementVolume_.resize(size_);
        elementVolume_ = 0;
    }
Kai Wendel's avatar
Kai Wendel committed
102

Nicolas Schwenck's avatar
Nicolas Schwenck committed
103
104
105
106
    /*!
     * \brief DOC ME!
     * \param materialLawParams DOC ME!
     */
107
108
    void prepareAnalytic()
    {
109
        const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(dummyElement_));
110

111
112
        swr_ = materialLawParams.swr();
        snr_ = materialLawParams.snr();
113
        porosity_ = problem_.spatialParams().porosity(dummyElement_);
114
115
        time_ = 0;
        satVec_ = swr_;
Kai Wendel's avatar
Kai Wendel committed
116

117
118
119
120
121
122
        for (int i = 1; i < pointNum_; i++)
        {
            satVec_[i] = satVec_[i - 1] + (1 - snr_ - swr_) / intervalNum_;
        }

        FluidState fluidState;
123
        fluidState.setTemperature(problem_.temperatureAtPos(dummyGlobal_));
124
125
        fluidState.setPressure(wPhaseIdx, problem_.referencePressureAtPos(dummyGlobal_));
        fluidState.setPressure(nPhaseIdx, problem_.referencePressureAtPos(dummyGlobal_));
126
127
        Scalar viscosityW = FluidSystem::viscosity(fluidState, wPhaseIdx);
        Scalar viscosityNW = FluidSystem::viscosity(fluidState, nPhaseIdx);
128
129
130
131
132
133
134
135

        for (int i = 0; i < pointNum_; i++)
        {
            fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
            fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
        }

        dfwdsw_ = 0;
Kai Wendel's avatar
Kai Wendel committed
136

137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
        for (int i = 1; i < intervalNum_; i++)
        {
            dfwdsw_[i] = (fractionalW_[i + 1] - fractionalW_[i - 1]) / (satVec_[i + 1] - satVec_[i - 1]);
        }

        for (int i = 0; i < pointNum_; i++)
        {
            if (dfwdsw_[i] > dfwdsw_[i + 1])
            {
                dfwdswmax_ = i;
                break;
            }
        }
    }

152
    void calcSatError()
153
154
155
156
    {
        error_ = 0;
        elementVolume_ = 0;
        ElementIterator eItEnd = problem_.gridView().template end<0> ();
Kai Wendel's avatar
Kai Wendel committed
157

158
159
160
161
162
163
164
165
166
167
        for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
        {
            // get entity
            const Element& element = *eIt;
            int index = problem_.variables().index(*eIt);
            elementVolume_[index] = element.geometry().volume();
        }

        for (int i = 0; i < size_; i++)
        {
168
            error_[i] = analyticSolution_[i] - problem_.variables().cellData(i).saturation(wPhaseIdx);
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
        }
    }

    void updateExSol()
    {
        // position of the fluid front
        xf_ = 0;
        for (int i = 0; i < pointNum_; i++)
        {
            xf_[i] = vTot_ * time_ / porosity_ * dfwdsw_[i];
        }
        int xhelp = pointNum_ / 3;
        int xhelpold = 0;
        int xhelpoldold = 0;
        int xfmax = 0;

        // position of maximum xf_
        for (int i = 0; i < pointNum_; i++)
        {
            if (xf_[i] > xf_[i + 1])
            {
                xfmax = i;
                break;
            }
        }

        // balancing of the areas ahead of the front and below the curve
        bool a = true;
        Scalar A1;
        Scalar A2;
        Scalar b;
        int xhelp2 = 0;

        while (a)
        {
            if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
                break;

            A1 = 0;

            for (int i = 0; i <= xhelp - 1; i++)
            {
                A1 += (satVec_[i] - swr_ + satVec_[i + 1] - swr_) * 0.5 * (xf_[i + 1] - xf_[i]);
            }

            A2 = 0;

            for (int i = xhelp; i <= xfmax - 1; i++)
            {
                A2 += (satVec_[xfmax] - satVec_[i] + satVec_[xfmax] - satVec_[i + 1]) * 0.5 * (xf_[i + 1] - xf_[i]);
            }

            b = xf_[xfmax];
            xhelp2 = xfmax;

            while (b > xf_[xhelp])
            {
                xhelp2 += 1;
                b = xf_[xhelp2];
            }

            for (int i = xfmax; i <= xhelp2; i++)
            {
                A2 += (satVec_[i] - satVec_[xfmax] + satVec_[i + 1] - satVec_[xfmax]) * 0.5 * (xf_[i] - xf_[i + 1]);
            }

            xhelpoldold = xhelpold;
            xhelpold = xhelp;

            if (fabs(A1) > fabs(A2))
            {
                xhelp = xhelp - 1;
            }

            if (fabs(A1) < fabs(A2))
            {
                xhelp = xhelp + 1;
            }

            if (xhelp == xhelpoldold)
            {
                a = false;
            }
        }

        // iterate over vertices and get analytic saturation solution
        ElementIterator eItEnd = problem_.gridView().template end<0> ();
        for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
        {
            // get global coordinate of cell center
            const GlobalPosition& globalPos = eIt->geometry().center();

            // find index of current vertex
            int index = problem_.variables().index(*eIt);

            // account for linear material law
            if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
            {
                if (globalPos[0] <= xf_[1])
                {
                    analyticSolution_[index] = 1 - snr_;
                }
                if (globalPos[0] > xf_[1])
                {
                    analyticSolution_[index] = swr_;
                }
            }
            // non-linear material law
            else
            {
                // find x_f next to global coordinate of the vertex
                int xnext = 0;
                for (int i = intervalNum_; i >= 0; i--)
                {
                    if (globalPos[0] < xf_[i])
                    {
                        xnext = i;
                        break;
                    }
                }

                // account for the area not yet reached by the front
                if (globalPos[0] > xf_[xhelp2])
                {
                    analyticSolution_[index] = swr_;
                    continue;
                }

                if (globalPos[0] <= xf_[xhelp2])
                {
                    analyticSolution_[index] = satVec_[xnext];
                    continue;
                }
            }
        }

305
306
        // call function to calculate the saturation error
        calcSatError();
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
    }

    void calculateAnalyticSolution()
    {
        time_ += problem_.timeManager().timeStepSize();
        updateExSol();
    }

    BlockVector AnalyticSolution() const
    {
        return analyticSolution_;
    }

    //Write saturation and pressure into file
    template<class MultiWriter>
    void addOutputVtkFields(MultiWriter &writer)
    {
        BlockVector *analyticSolution = writer.allocateManagedBuffer (size_);
        BlockVector *error = writer.allocateManagedBuffer (size_);

        *analyticSolution = analyticSolution_;
        *error = error_;

        writer.attachCellData(*analyticSolution, "saturation (exact solution)");
        writer.attachCellData(*error, "error_");
    }

    //! Construct an IMPES object.
    BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) :
        problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), vTot_(totalVelocity), dummyElement_(
337
338
339
340
341
                *(problem_.gridView().template begin<0> ()))
    {
        dummyGlobal_ = 0.0;
        dummyGlobal_[0] = 1.0;
    }
Markus Wolff's avatar
Markus Wolff committed
342
343

    void initialize()
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    {
        initializeAnalytic();
        prepareAnalytic();
    }

private:
    Problem& problem_;

    BlockVector analyticSolution_;
    BlockVector error_;
    BlockVector elementVolume_;

    Scalar time_;
    int size_;
    Scalar swr_;
    Scalar snr_;
    Scalar porosity_;
    Scalar vTot_;
Kai Wendel's avatar
Kai Wendel committed
362
    static constexpr int intervalNum_ = 1000, pointNum_ = intervalNum_ + 1;
363
364
365
366
367
368
    Dune::FieldVector<Scalar, pointNum_> satVec_;
    Dune::FieldVector<Scalar, pointNum_> fractionalW_;
    Dune::FieldVector<Scalar, pointNum_> dfwdsw_;
    Dune::FieldVector<Scalar, pointNum_> xf_;
    int dfwdswmax_;
    const Element& dummyElement_;
369
    GlobalPosition dummyGlobal_;
370
371

};
Kai Wendel's avatar
Kai Wendel committed
372
373
374
375

} // end namespace Dumux

#endif // DUMUX_BUCKLEYLEVERETT_ANALYTICAL_HH