henry1p2cproblem.hh 7.59 KB
Newer Older
Lena Walter's avatar
Lena Walter committed
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.                      *
 *                                                                           *
Lena Walter's avatar
Lena Walter committed
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            *
Lena Walter's avatar
Lena Walter committed
14
15
16
17
18
19
20
 *   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
Lena Walter's avatar
Lena Walter committed
21
22
 * \brief Definition of the 1p2c Henry problem
  *
Lena Walter's avatar
Lena Walter committed
23
24
25
26
 */
#ifndef DUMUX_HENRY1P2C_PROBLEM_HH
#define DUMUX_HENRY1P2C_PROBLEM_HH

27
28
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
29
#include <dumux/common/boundarytypes.hh>
Kai Wendel's avatar
Kai Wendel committed
30
#include <dumux/porousmediumflow/problem.hh>
Lena Walter's avatar
Lena Walter committed
31

Kai Wendel's avatar
Kai Wendel committed
32
namespace Dumux {
Lena Walter's avatar
Lena Walter committed
33
34

template <class TypeTag>
Kai Wendel's avatar
Kai Wendel committed
35
class Henry1p2cProblem : public PorousMediumFlowProblem<TypeTag>
Lena Walter's avatar
Lena Walter committed
36
{
37
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
38
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
39
40
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
41
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
42
    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
43
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
Kai Wendel's avatar
Kai Wendel committed
44
45
    using Element = typename GridView::template Codim<0>::Entity;
    using ElementIterator = typename GridView::template Codim<0>::Iterator;
46
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
Kai Wendel's avatar
Kai Wendel committed
47
    using ParentType = PorousMediumFlowProblem<TypeTag>;
48
49
50
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
    using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
51
52
53
54
55
56
57
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;

     // copy some indices for convenience
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
      // indices of the equations
    static constexpr int conti0EqIdx = Indices::conti0EqIdx;

Kai Wendel's avatar
Kai Wendel committed
58
59
60
    // Grid and world dimension
    static constexpr int dim = GridView::dimension;
    static constexpr int dimWorld = GridView::dimensionworld;
61

Kai Wendel's avatar
Kai Wendel committed
62
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
Lena Walter's avatar
Lena Walter committed
63
64

public:
Kai Wendel's avatar
Kai Wendel committed
65
66
67
    Henry1p2cProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
    : ParentType(fvGridGeometry)
    {
68
      FluidSystem::init();
Kai Wendel's avatar
Kai Wendel committed
69
70
      freshWaterFluxRate_= getParam<Scalar>("Problem.freshWaterFluxRate");
    }
Lena Walter's avatar
Lena Walter committed
71
72
73
74
75
76
77
78
79
80
81
82

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

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 36 degrees Celsius.
     */
    Scalar temperature() const
Kai Wendel's avatar
Kai Wendel committed
83
84
85
    {
        return 273.15 + 20; // in [K]
    };
Lena Walter's avatar
Lena Walter committed
86
87
88
89
90
91
92
93
94
95
96
97

    // \}

    /*!
     * \name Boundary conditions
     */
    // \{

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     */
Kai Wendel's avatar
Kai Wendel committed
98
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
Lena Walter's avatar
Lena Walter committed
99
    {
Kai Wendel's avatar
Kai Wendel committed
100
        BoundaryTypes values;
Lena Walter's avatar
Lena Walter committed
101
102
103

        values.setAllDirichlet();

Kai Wendel's avatar
Kai Wendel committed
104
        if (globalPos[0]<eps_  || globalPos[1]<eps_ || globalPos[1]>1-eps_)
Lena Walter's avatar
Lena Walter committed
105
106
107
108
        {
            values.setAllNeumann();
        }

Kai Wendel's avatar
Kai Wendel committed
109
        return values;
Lena Walter's avatar
Lena Walter committed
110
111
112
113
114
115
116
117
    }

    /*!
     * \brief Evaluate the boundary conditions for a dirichlet
     *        boundary segment.
     *
     * For this method, the \a values parameter stores primary variables.
     */
Kai Wendel's avatar
Kai Wendel committed
118
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
Lena Walter's avatar
Lena Walter committed
119
    {
Kai Wendel's avatar
Kai Wendel committed
120
        PrimaryVariables values;
Lena Walter's avatar
Lena Walter committed
121
122

        initial_(values, globalPos);
Kai Wendel's avatar
Kai Wendel committed
123
        if (globalPos[0]>2-eps_ && globalPos[1]<0.8+eps_)
Lena Walter's avatar
Lena Walter committed
124
        {
125
            values[FluidSystem::SaltIdx] =  0.03922;
Lena Walter's avatar
Lena Walter committed
126
127
        }

Kai Wendel's avatar
Kai Wendel committed
128
        return values;
Lena Walter's avatar
Lena Walter committed
129
130
131
132
133
134
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
135
136
137
138
139
140
141
142
     * \param values The neumann values for the conservation equations in units of
     *                 \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
     * \param elemFluxVarsCache Flux variables caches for all faces in stencil
     * \param scvf The sub control volume face
     *
Lena Walter's avatar
Lena Walter committed
143
144
145
146
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each component. Negative values mean
     * influx.
     */
147
148
149
    NumEqVector neumann(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
150
                        const ElementFluxVariablesCache& elemFluxVarsCache,
151
                        const SubControlVolumeFace& scvf) const
Lena Walter's avatar
Lena Walter committed
152
    {
153
        const auto globalPos = element.geometry().corner(scvf.insideScvIdx());
Lena Walter's avatar
Lena Walter committed
154

155
        NumEqVector values(0.0);
Kai Wendel's avatar
Kai Wendel committed
156
157
        if (globalPos[0]<eps_)
        {
158
            values[conti0EqIdx] = -freshWaterFluxRate_; //-6.6E-2; //6.6e-5*m2/s*1000kg/m^3
Kai Wendel's avatar
Kai Wendel committed
159
160
        }
        return values;
Lena Walter's avatar
Lena Walter committed
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
    }
    // \}

    /*!
     * \name Volume terms
     */
    // \{

    /*!
     * \brief Evaluate the source term for all phases within a given
     *        sub-control-volume.
     *
     * For this method, the \a values parameter stores the rate mass
     * of a component is generated or annihilate per volume
     * unit. Positive values mean that mass is created, negative ones
     * mean that it vanishes.
     */
Kai Wendel's avatar
Kai Wendel committed
178
    NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
Lena Walter's avatar
Lena Walter committed
179
    {
Kai Wendel's avatar
Kai Wendel committed
180
        return NumEqVector(0.0);
Lena Walter's avatar
Lena Walter committed
181
182
183
184
185
186
187
188
    }

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
Kai Wendel's avatar
Kai Wendel committed
189
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
Lena Walter's avatar
Lena Walter committed
190
    {
Kai Wendel's avatar
Kai Wendel committed
191
        PrimaryVariables values;
Lena Walter's avatar
Lena Walter committed
192
193
194

        initial_(values, globalPos);

Kai Wendel's avatar
Kai Wendel committed
195
196
        return values;
    }
Lena Walter's avatar
Lena Walter committed
197
198
199
200
    // \}

private:
    // the internal method for the initial condition
Kai Wendel's avatar
Kai Wendel committed
201
    void initial_(PrimaryVariables &values, const GlobalPosition &globalPos) const
Lena Walter's avatar
Lena Walter committed
202
203
    {
        Scalar densityW = 1025;
Kai Wendel's avatar
Kai Wendel committed
204
205
        values[Indices::pressureIdx] = 1.0133e5+(depthBOR_-globalPos[1])*densityW*9.81; //initial condition for the pressure
        values[FluidSystem::SaltIdx] = 0.0; //initial condition for the salt molefraction
Lena Walter's avatar
Lena Walter committed
206
207
208
    }

    Scalar freshWaterFluxRate_;
209
210
    static constexpr Scalar eps_ = 1e-6;
    static constexpr Scalar depthBOR_ = 1.0; // [m]
Kai Wendel's avatar
Kai Wendel committed
211
212
213
214
};

} // end namespace Dumux

215
#endif