test_diffusion_problem.hh 10.2 KB
Newer Older
Bernd Flemisch's avatar
Bernd Flemisch committed
1
2
// $Id: test_diffusion_problem.hh 3655 2010-05-26 17:13:50Z bernd $
/*****************************************************************************
Bernd Flemisch's avatar
Bernd Flemisch committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
*   Copyright (C) 2007-2008 by Klaus Mosthaf                                *
*   Copyright (C) 2007-2008 by Bernd Flemisch                               *
*   Copyright (C) 2008-2009 by Andreas Lauser                               *
*   Institute of Hydraulic Engineering                                      *
*   University of Stuttgart, Germany                                        *
*   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
*                                                                           *
*   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, as long as this copyright notice    *
*   is included in its original form.                                       *
*                                                                           *
*   This program is distributed WITHOUT ANY WARRANTY.                       *
*****************************************************************************/
Bernd Flemisch's avatar
Bernd Flemisch committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#ifndef DUMUX_TEST_2P_PROBLEM_HH
#define DUMUX_TEST_2P_PROBLEM_HH

#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif

#include <dune/grid/yaspgrid.hh>
#include <dune/grid/sgrid.hh>

#include <dumux/material/components/unit.hh>

#include <dumux/decoupled/2p/diffusion/fvmpfa/mpfaproperties.hh>
#include <dumux/decoupled/2p/diffusion/diffusionproblem2p.hh>
#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
33
#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh>
Bernd Flemisch's avatar
Bernd Flemisch committed
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh>

#include "test_diffusion_spatialparams.hh"

namespace Dumux
{

struct FileNames
{
    enum
    {
        TPFAName = 0, MPFAName = 1, MimeticName = 2
    };
};

template<class TypeTag>
class TestDiffusionProblem;

//////////
// Specify the properties
//////////
namespace Properties
{
NEW_TYPE_TAG(DiffusionTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties))
Bernd Flemisch's avatar
Bernd Flemisch committed
58
        ;
Bernd Flemisch's avatar
Bernd Flemisch committed
59
60
61
62

// Set the grid type
SET_PROP(DiffusionTestProblem, Grid)
{//    typedef Dune::YaspGrid<2> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
63
    typedef Dune::SGrid<2, 2> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
64
65
66
67
68
69
70
71
};

SET_PROP(DiffusionTestProblem, PressurePreconditioner)
{
private:
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) Vector;
public:
Bernd Flemisch's avatar
Bernd Flemisch committed
72
    typedef Dune::SeqILUn<Matrix, Vector, Vector> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
73
74
75
76
77
78
};

//SET_INT_PROP(DiffusionTestProblem, VelocityFormulation,
//        GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);

SET_INT_PROP(DiffusionTestProblem, PressureFormulation,
Bernd Flemisch's avatar
Bernd Flemisch committed
79
        GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal);
Bernd Flemisch's avatar
Bernd Flemisch committed
80
81
82
83
84

// Set the wetting phase
SET_PROP(DiffusionTestProblem, WettingPhase)
{
private:
Bernd Flemisch's avatar
Bernd Flemisch committed
85
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
Bernd Flemisch's avatar
Bernd Flemisch committed
86
public:
Bernd Flemisch's avatar
Bernd Flemisch committed
87
    typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
Bernd Flemisch's avatar
Bernd Flemisch committed
88
89
90
91
92
93
};

// Set the non-wetting phase
SET_PROP(DiffusionTestProblem, NonwettingPhase)
{
private:
Bernd Flemisch's avatar
Bernd Flemisch committed
94
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
Bernd Flemisch's avatar
Bernd Flemisch committed
95
public:
Bernd Flemisch's avatar
Bernd Flemisch committed
96
    typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
Bernd Flemisch's avatar
Bernd Flemisch committed
97
98
};

99
// Set the spatial parameters
Bernd Flemisch's avatar
Bernd Flemisch committed
100
101
102
SET_PROP(DiffusionTestProblem, SpatialParameters)
{
private:
Bernd Flemisch's avatar
Bernd Flemisch committed
103
104
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
Bernd Flemisch's avatar
Bernd Flemisch committed
105
106

public:
Bernd Flemisch's avatar
Bernd Flemisch committed
107
    typedef Dumux::TestDiffusionSpatialParams<TypeTag> type;
Bernd Flemisch's avatar
Bernd Flemisch committed
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
};

// Enable gravity
SET_BOOL_PROP(DiffusionTestProblem, EnableGravity, false);

NEW_PROP_TAG( FileName );

// set the types for the 2PFA FV method
NEW_TYPE_TAG(FVVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
SET_INT_PROP(FVVelocity2PTestProblem, FileName, FileNames::TPFAName);
SET_INT_PROP(FVVelocity2PTestProblem, IterationNumberPreconditioner, 0);
SET_TYPE_PROP(FVVelocity2PTestProblem, Model, Dumux::FVVelocity2P<TTAG(FVVelocity2PTestProblem)>);
SET_TYPE_PROP(FVVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVVelocity2PTestProblem)>);

// set the types for the MPFA-O FV method
NEW_TYPE_TAG(FVMPFAOVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
SET_INT_PROP(FVMPFAOVelocity2PTestProblem, FileName, FileNames::MPFAName);
SET_INT_PROP(FVMPFAOVelocity2PTestProblem, IterationNumberPreconditioner, 1);
126
SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Model, Dumux::FVMPFAOVelocity2P<TypeTag>);
Bernd Flemisch's avatar
Bernd Flemisch committed
127
128
129
130
131
SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVMPFAOVelocity2PTestProblem)>);

// set the types for the mimetic FD method
NEW_TYPE_TAG(MimeticPressure2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
SET_INT_PROP(MimeticPressure2PTestProblem, FileName, FileNames::MimeticName);
Bernd Flemisch's avatar
Bernd Flemisch committed
132
SET_INT_PROP(MimeticPressure2PTestProblem, IterationNumberPreconditioner, 0);
Bernd Flemisch's avatar
Bernd Flemisch committed
133
134
135
136
137
SET_TYPE_PROP(MimeticPressure2PTestProblem, Model, Dumux::MimeticPressure2P<TTAG(MimeticPressure2PTestProblem)>);
SET_TYPE_PROP(MimeticPressure2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(MimeticPressure2PTestProblem)>);
}

/*!
Bernd Flemisch's avatar
Bernd Flemisch committed
138
139
* \ingroup DecoupledProblems
*/
Bernd Flemisch's avatar
Bernd Flemisch committed
140
141
142
template<class TypeTag = TTAG(DiffusionTestProblem)>
class TestDiffusionProblem: public DiffusionProblem2P<TypeTag, TestDiffusionProblem<TypeTag> >
{
Bernd Flemisch's avatar
Bernd Flemisch committed
143
144
145
    typedef TestDiffusionProblem<TypeTag> ThisType;
    typedef DiffusionProblem2P<TypeTag, ThisType> ParentType;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
Bernd Flemisch's avatar
Bernd Flemisch committed
146

Bernd Flemisch's avatar
Bernd Flemisch committed
147
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
Bernd Flemisch's avatar
Bernd Flemisch committed
148

Bernd Flemisch's avatar
Bernd Flemisch committed
149
150
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
Bernd Flemisch's avatar
Bernd Flemisch committed
151

Bernd Flemisch's avatar
Bernd Flemisch committed
152
153
154
155
    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };
Bernd Flemisch's avatar
Bernd Flemisch committed
156

Bernd Flemisch's avatar
Bernd Flemisch committed
157
158
159
160
    enum
    {
        wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
    };
Bernd Flemisch's avatar
Bernd Flemisch committed
161

Bernd Flemisch's avatar
Bernd Flemisch committed
162
    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
Bernd Flemisch's avatar
Bernd Flemisch committed
163

Bernd Flemisch's avatar
Bernd Flemisch committed
164
165
166
167
    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::Intersection Intersection;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
    typedef Dune::FieldVector<Scalar, dim> LocalPosition;
Bernd Flemisch's avatar
Bernd Flemisch committed
168
169

public:
Bernd Flemisch's avatar
Bernd Flemisch committed
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
    TestDiffusionProblem(const GridView &gridView, const double delta = 1.0) :
        ParentType(gridView), delta_(delta)
    {
        this->variables().saturation() = 1.0;
        this->spatialParameters().setDelta(delta_);
    }

    /*!
    * \name Problem parameters
    */
    // \{

    /*!
    * \brief The problem name.
    *
    * This is used as a prefix for files generated by the simulation.
    */
    const char *name() const
    {
        switch (GET_PROP_VALUE(TypeTag, PTAG(FileName)))
        {
        case FileNames::TPFAName:
            return "fvdiffusion";
        case FileNames::MPFAName:
            return "fvmpfaodiffusion";
        case FileNames::MimeticName:
            return "mimeticdiffusion";
        }
        return "test_diffusion";
    }

201
    bool shouldWriteRestartFile() const
202
    { return false; }
Bernd Flemisch's avatar
Bernd Flemisch committed
203
204
205
206
207
208
209
210
211
212

    /*!
    * \brief Returns the temperature within the domain.
    *
    * This problem assumes a temperature of 10 degrees Celsius.
    */
    Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
    {
        return 273.15 + 10; // -> 10°C
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
213

Bernd Flemisch's avatar
Bernd Flemisch committed
214
    // \}
Bernd Flemisch's avatar
Bernd Flemisch committed
215

Bernd Flemisch's avatar
Bernd Flemisch committed
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 referencePressure(const GlobalPosition& globalPos, const Element& element) const
    {
        return 1e5; // -> 10°C
    }

    std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
        {
        double pi = 4.0*atan(1.0);
        double rt = globalPos[0]*globalPos[0]+globalPos[1]*globalPos[1];
        double ux = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
        double uy = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
        double kxx = (delta_*globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])/rt;
        double kxy = -(1.0 - delta_)*globalPos[0]*globalPos[1]/rt;
        double kyy = (globalPos[0]*globalPos[0] + delta_*globalPos[1]*globalPos[1])/rt;
        double f0 = sin(pi*globalPos[0])*sin(pi*globalPos[1])*pi*pi*(1.0 + delta_)*(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])
        + cos(pi*globalPos[0])*sin(pi*globalPos[1])*pi*(1.0 - 3.0*delta_)*globalPos[0]
                                                                                    + cos(pi*globalPos[1])*sin(pi*globalPos[0])*pi*(1.0 - 3.0*delta_)*globalPos[1]
                                                                                                                                                                + cos(pi*globalPos[1])*cos(pi*globalPos[0])*2.0*pi*pi*(1.0 - delta_)*globalPos[0]*globalPos[1];

        std::vector<double> result(2, 0.0);
        result[wPhaseIdx]=(f0 + 2.0*(globalPos[0]*(kxx*ux + kxy*uy) + globalPos[1]*(kxy*ux + kyy*uy)))/rt;

        return (result);
        }

    typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
    {
        return BoundaryConditions::dirichlet;
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
245

Bernd Flemisch's avatar
Bernd Flemisch committed
246
247
248
249
    Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
    {
        return (exact(globalPos));
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
250

Bernd Flemisch's avatar
Bernd Flemisch committed
251
252
253
254
    typename BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
    {
        return BoundaryConditions::dirichlet;
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
255

Bernd Flemisch's avatar
Bernd Flemisch committed
256
257
258
259
    Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
    {
        return 1.0;
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
260

Bernd Flemisch's avatar
Bernd Flemisch committed
261
262
263
    std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
        {
        std::vector<Scalar> neumannFlux(2, 0.0);
Bernd Flemisch's avatar
Bernd Flemisch committed
264

Bernd Flemisch's avatar
Bernd Flemisch committed
265
266
        return neumannFlux;
        }
Bernd Flemisch's avatar
Bernd Flemisch committed
267

Bernd Flemisch's avatar
Bernd Flemisch committed
268
269
270
    Scalar exact (const GlobalPosition& globalPos) const
    {
        double pi = 4.0*atan(1.0);
Bernd Flemisch's avatar
Bernd Flemisch committed
271

Bernd Flemisch's avatar
Bernd Flemisch committed
272
273
        return (sin(pi*globalPos[0])*sin(pi*globalPos[1]));
    }
Bernd Flemisch's avatar
Bernd Flemisch committed
274

Bernd Flemisch's avatar
Bernd Flemisch committed
275
276
277
278
279
280
    Dune::FieldVector<Scalar,dim> exactGrad (const GlobalPosition& globalPos) const
        {
        Dune::FieldVector<Scalar,dim> grad(0);
        double pi = 4.0*atan(1.0);
        grad[0] = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
        grad[1] = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
Bernd Flemisch's avatar
Bernd Flemisch committed
281

Bernd Flemisch's avatar
Bernd Flemisch committed
282
283
        return grad;
        }
Bernd Flemisch's avatar
Bernd Flemisch committed
284
285

private:
Bernd Flemisch's avatar
Bernd Flemisch committed
286
    double delta_;
Bernd Flemisch's avatar
Bernd Flemisch committed
287
288
289
290
};
} //end namespace

#endif