test_1pproblem.hh 8.73 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
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
 *   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/>.   *
18
*****************************************************************************/
Bernd Flemisch's avatar
Bernd Flemisch committed
19
20
21
22
23
24
25
/*!
 * \file
 *
 * \brief test problem for the decoupled one-phase model.
 */
#ifndef DUMUX_TEST_1P_PROBLEM_HH
#define DUMUX_TEST_1P_PROBLEM_HH
26
27

#include <dune/grid/yaspgrid.hh>
28
#include <dumux/io/cubegridcreator.hh>
29
30
31
32

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

33
#include <dumux/decoupled/1p/diffusion/fv/fvpressureproperties1p.hh>
34
#include <dumux/decoupled/1p/diffusion/diffusionproblem1p.hh>
35
#include <dumux/decoupled/common/fv/fvvelocity.hh>
36

37
#include "test_1pspatialparams.hh"
38
39
40
41
42
43
44
45
46
47
48
49

namespace Dumux
{

template<class TypeTag>
class TestProblemOneP;

//////////
// Specify the properties
//////////
namespace Properties
{
50
51
52
53
NEW_TYPE_TAG(TestProblemOneP, INHERITS_FROM(FVPressureOneP));

// set the GridCreator property
SET_TYPE_PROP(TestProblemOneP, GridCreator, CubeGridCreator<TypeTag>);
54
55

// Set the grid type
56
SET_TYPE_PROP(TestProblemOneP, Grid, Dune::YaspGrid<2>);
57
58
59
60
61

// Set the wetting phase
SET_PROP(TestProblemOneP, Fluid)
{
private:
62
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
63
64
65
66
67
public:
    typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
};

// Set the spatial parameters
68
SET_TYPE_PROP(TestProblemOneP, SpatialParams, Dumux::TestOnePSpatialParams<TypeTag>);
69
70

//Set the problem
71
SET_TYPE_PROP(TestProblemOneP, Problem, Dumux::TestProblemOneP<TypeTag>);
72
73
74
}

/*!
Benjamin Faigle's avatar
Benjamin Faigle committed
75
 * \ingroup IMPETtests
76
77
78
 *
 * \brief test problem for the decoupled one-phase model.
 */
79
template<class TypeTag>
80
class TestProblemOneP: public DiffusionProblem1P<TypeTag >
81
{
82
    typedef DiffusionProblem1P<TypeTag> ParentType;
83
84
    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;

85
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
86

87
    typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
88

89
90
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
91

92
93
94
95
96
    enum
    {
        dim = GridView::dimension, dimWorld = GridView::dimensionworld
    };

97
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
98
99
100
101

    typedef typename GridView::Traits::template Codim<0>::Entity Element;
    typedef typename GridView::Intersection Intersection;
    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
102
    typedef Dune::FieldVector<Scalar, dim> LocalPosition;
103
    typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
104
    typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
105

106
107

public:
108
109
    TestProblemOneP(TimeManager &timeManager, const GridView &gridView) :
        ParentType(gridView), velocity_(*this)
110
    {
111
        delta_ = 1e-3 ;
112
      
113
            if (ParameterTree::tree().hasKey("Problem.Delta"))
114
            delta_       = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, Delta);
115
            int numRefine;
116
            numRefine = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Grid, NumRefine);
117
118
            GridCreator::grid().globalRefine(numRefine);

Bernd Flemisch's avatar
Bernd Flemisch committed
119
        this->spatialParams().initialize(delta_);
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
    }

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

    /*!
    * \brief The problem name.
    *
    * This is used as a prefix for files generated by the simulation.
    */
    const char *name() const
    {
        return "test_1p";
    }

    bool shouldWriteRestartFile() const
    { return false; }

140
141
142
143
144
145
    void addOutputVtkFields()
    {
        velocity_.calculateVelocity();
        velocity_.addOutputVtkFields(this->resultWriter());
    }

146
147
148
149
150
    /*!
    * \brief Returns the temperature within the domain.
    *
    * This problem assumes a temperature of 10 degrees Celsius.
    */
151
    Scalar temperatureAtPos(const GlobalPosition& globalPos) const
152
153
154
155
156
157
    {
        return 273.15 + 10; // -> 10°C
    }

    // \}

Markus Wolff's avatar
Markus Wolff committed
158
    //! Returns the reference pressure for evaluation of constitutive relations
159
    Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
160
161
162
163
    {
        return 1e5; // -> 10°C
    }

Markus Wolff's avatar
Markus Wolff committed
164
    //!source term [kg/(m^3 s)]
165
    void source(PrimaryVariables &values, const Element& element) const
166
    {
167
168
        values = 0;
        values = integratedSource_(element, 4);
169
    }
170

Markus Wolff's avatar
Markus Wolff committed
171
172
173
174
175
    /*!
    * \brief Returns the type of boundary condition.
    *
    * BC can be dirichlet (pressure) or neumann (flux).
    */
176
177
    void boundaryTypes(BoundaryTypes &bcType,
            const Intersection& intersection) const
178
    {
179
        bcType.setAllDirichlet();
180
181
    }

Markus Wolff's avatar
Markus Wolff committed
182
    //! return dirichlet condition  (pressure, [Pa])
183
184
    void dirichletAtPos(PrimaryVariables &values,
                        const GlobalPosition &globalPos) const
185
    {
186
        values = exact(globalPos);
187
188
    }

189

Markus Wolff's avatar
Markus Wolff committed
190
    //! return neumann condition  (flux, [kg/(m^2 s)])
191
    void neumann(PrimaryVariables &values, const Intersection& intersection) const
192
    {
193
        values = 0;
194
    }
195

Markus Wolff's avatar
Markus Wolff committed
196
private:
197
198
199
200
201
202
203
204
    Scalar exact (const GlobalPosition& globalPos) const
    {
        double pi = 4.0*atan(1.0);

        return (sin(pi*globalPos[0])*sin(pi*globalPos[1]));
    }

    Dune::FieldVector<Scalar,dim> exactGrad (const GlobalPosition& globalPos) const
205
    {
206
207
208
209
210
211
        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]);

        return grad;
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
    Scalar integratedSource_(const Element& element, int integrationPoints) const
    {
        Scalar source = 0.;
        LocalPosition localPos(0.0);
        GlobalPosition globalPos(0.0);
        Scalar halfInterval = 1.0/double(integrationPoints)/2.;
        for (int i = 1; i <= integrationPoints; i++)
        {
            for (int j = 1; j <= integrationPoints; j++)
            {
                localPos[0] = double(i)/double(integrationPoints) - halfInterval;
                localPos[1] = double(j)/double(integrationPoints) - halfInterval;
                globalPos = element.geometry().global(localPos);
                source += 1./(integrationPoints*integrationPoints) * evaluateSource_(globalPos);
            }
        }

        return source;
    }

    Scalar evaluateSource_(const GlobalPosition& globalPos) const
    {
        Scalar temp = temperatureAtPos(globalPos);
        Scalar referencePress = referencePressureAtPos(globalPos);

        Scalar pi = 4.0 * atan(1.0);
        Scalar x = globalPos[0];
        Scalar y = globalPos[1];

        Scalar dpdx = pi * cos(pi * x) * sin(pi * y);
        Scalar dpdy = pi * sin(pi * x) * cos(pi * y);
        Scalar dppdxx = -pi * pi * sin(pi * x) * sin(pi * y);
        Scalar dppdxy = pi * pi * cos(pi * x) * cos(pi * y);
        Scalar dppdyx = dppdxy;
        Scalar dppdyy = dppdxx;
        Scalar kxx = (delta_* x*x + y*y)/(x*x + y*y);
        Scalar kxy = -(1.0 - delta_) * x * y / (x*x + y*y);
        Scalar kyy = (x*x + delta_*y*y)/(x*x + y*y);
        Scalar dkxxdx = 2 * x * y*y * (delta_ - 1.0)/((x*x + y*y) * (x*x + y*y));
        Scalar dkyydy = 2 * x*x * y * (delta_ - 1.0)/((x*x + y*y) * (x*x + y*y));
        Scalar dkxydx = (1.0 - delta_) * y * (x*x - y*y) /((x*x + y*y) * (x*x + y*y));
        Scalar dkxydy = (1.0 - delta_) * x * (y*y - x*x) /((x*x + y*y) * (x*x + y*y));

        Scalar fx = dkxxdx * dpdx + kxx * dppdxx + dkxydx * dpdy + kxy * dppdyx;
        Scalar fy = dkxydy * dpdx + kxy * dppdxy + dkyydy * dpdy + kyy * dppdyy;

        return -(fx + fy) / Fluid::viscosity(temp, referencePress) * Fluid::density(temp, referencePress);
    }

263
    double delta_;
264
    Dumux::FVVelocity<TypeTag, typename GET_PROP_TYPE(TypeTag, Velocity) > velocity_;
265
266
267
268
};
} //end namespace

#endif