matrixproblem.hh 10.7 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
28
29
30
31
32
// -*- 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
 * \ingroup MultiDomain
 * \ingroup MultiDomainFacet
 * \ingroup TwoPTests
 * \brief The sub-problem for the matrix domain in the exercise on two-phase flow in fractured porous media.
 */
#ifndef DUMUX_COURSE_FRACTURESEXERCISE_MATRIX_PROBLEM_HH
#define DUMUX_COURSE_FRACTURESEXERCISE_MATRIX_PROBLEM_HH

// we need this in this test in order to define the domain
// id of the fracture problem (see function interiorBoundaryTypes())
#include <dune/common/indices.hh>

33
// include the base problem and properties we inherit from
34
#include <dumux/porousmediumflow/problem.hh>
35
#include <dumux/common/properties.hh>
36
37
38
39
40
41
42
43
44
45
46
47
48
49

namespace Dumux {

/*!
 * \ingroup MultiDomain
 * \ingroup MultiDomainFacet
 * \ingroup TwoPTests
 * \brief The sub-problem for the matrix domain in the exercise on two-phase flow in fractured porous media.
 */
template<class TypeTag>
class MatrixSubProblem : public PorousMediumFlowProblem<TypeTag>
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;

50
    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
51
52
53
    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
54
55
56
57
58
59
60
61
62
63
64
65
    using PrimaryVariables = typename GridVariables::PrimaryVariables;
    using Scalar = typename GridVariables::Scalar;

    using FVGridGeometry = typename GridVariables::GridGeometry;
    using FVElementGeometry = typename FVGridGeometry::LocalView;
    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
    using GridView = typename FVGridGeometry::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

    // some indices for convenience
66
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
    enum
    {
        pressureIdx = Indices::pressureIdx,
        saturationIdx = Indices::saturationIdx
    };

public:
    //! The constructor
    MatrixSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
                     std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
                     const std::string& paramGroup = "")
    : ParentType(fvGridGeometry, spatialParams, paramGroup)
    , boundaryOverPressure_(getParamFromGroup<Scalar>(paramGroup, "Problem.BoundaryOverPressure"))
    , boundarySaturation_(getParamFromGroup<Scalar>(paramGroup, "Problem.BoundarySaturation"))
    , isExercisePartA_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartA"))
    , isExercisePartB_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartB"))
    , isExercisePartC_(getParamFromGroup<bool>(paramGroup, "Problem.IsExercisePartC"))
    {
        // initialize the fluid system, i.e. the tabulation
        // of water properties. Use the default p/T ranges.
87
        using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
        FluidSystem::init();

        // there can only be one exercise part active
        if (isExercisePartA_ + isExercisePartB_ + isExercisePartC_ > 1)
            DUNE_THROW(Dune::InvalidStateException, "Invalid combination of activated exercise parts");
    }

    //! Specifies the type of boundary condition at a given position
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
    {
        BoundaryTypes values;

        // in the unmodified state, we consider buoancy-driven upwards migration
        // of nitrogen and set Dirichlet BCs on the top and bottom boundary
        if (!isExercisePartA_ && !isExercisePartB_ && !isExercisePartC_)
        {
            values.setAllNeumann();
105
            if (globalPos[1] < 1e-6 || globalPos[1] > this->gridGeometry().bBoxMax()[1] - 1e-6)
106
107
108
109
110
111
112
                values.setAllDirichlet();
        }
        else
        {
            // for exercise part a,b & c we use no-flow boundaries everywhere
            // except for the lower quarter of the left and the upper quarter of the right boundary
            values.setAllNeumann();
113
            if ( (globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-6 && globalPos[1] > 75.0) ||
114
115
116
117
118
119
120
                 (globalPos[0] < 1e-6 && globalPos[1] < 25.0) )
                values.setAllDirichlet();
        }

        return values;
    }

121
    //! Specifies the type of interior boundary condition to be used on a sub-control volume face
122
123
124
125
126
127
128
129
    BoundaryTypes interiorBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
    {
        BoundaryTypes values;

        // Here we set the type of condition to be used on faces that coincide
        // with a fracture. If Neumann is specified, a flux continuity condition
        // on the basis of the normal fracture permeability is evaluated. If this
        // permeability is lower than that of the matrix, this approach is able to
130
        // represent the resulting pressure jump across the fracture. If Dirichlet is set,
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
        // the pressure jump across the fracture is neglected and the pressure inside
        // the fracture is directly applied at the interface between fracture and matrix.
        // This assumption is justified for highly-permeable fractures, but lead to erroneous
        // results for low-permeable fractures.
        // Here, we consider "open" fractures for which we cannot define a normal permeability
        // and for which the pressure jump across the fracture is neglectable. Thus, we set
        // the interior boundary conditions to Dirichlet.
        // IMPORTANT: Note that you will never be asked to set any values at the interior boundaries!
        //            This simply chooses a different interface condition!
        values.setAllDirichlet();

        if (isExercisePartB_)
            values.setAllNeumann();
        else if (isExercisePartC_)
        {
            // we need to obtain the domain marker of the fracture element that is coupled to this face
            // therefore we first get the fracture problem from the coupling manager. For this test, we
            // know that the fracture domain id is 1 (see main file)
149
150
            static constexpr auto fractureDomainId = Dune::index_constant<1>();
            const auto& fractureProblem = couplingManager().problem(fractureDomainId);
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

            // use helper function in coupling manager to obtain the element this face couples to
            const auto fractureElement = couplingManager().getLowDimElement(element, scvf);

            // obtain marker from the spatial params of the fracture problem
            const auto fractureElementMarker = fractureProblem.spatialParams().getElementDomainMarker(fractureElement);

            // now define Neumann coupling on those elements that are defined as blocking
            // fractures, i.e. marker == 2, see .geo file in grids folder
            if (fractureElementMarker == 2)
                values.setAllNeumann();
        }

        return values;
    }

    //! evaluates the Dirichlet boundary condition for a given position
    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
    {
        // initialize values with the initial conditions
        auto values = initialAtPos(globalPos);

        // In the unmodified state of the exercise nitrogen is in contact
        // with the domain on the center half of the lower boundary
        if (!isExercisePartA_ && !isExercisePartB_ && !isExercisePartC_)
        {
            if (globalPos[1] < 1e-6 && globalPos[0] > 25.0 && globalPos[0] < 75.0)
                values[saturationIdx] = boundarySaturation_;
        }
        // exercise parts a, b & c
        else
        {
            // apply overpressure on the right Dirichlet boundary segment
184
            if (globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-6)
185
186
187
188
189
190
191
192
193
194
195
196
197
198
            {
                values[pressureIdx] += boundaryOverPressure_;
                values[saturationIdx] = boundarySaturation_;
            }
        }

        return values;
    }

    //! evaluate the initial conditions
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
    {
        // For the grid used here, the height of the domain is equal
        // to the maximum y-coordinate
199
        const auto domainHeight = this->gridGeometry().bBoxMax()[1];
200
201

        // we assume a constant water density of 1000 for initial conditions!
202
        const auto& g = this->spatialParams().gravity(globalPos);
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
        PrimaryVariables values;
        Scalar densityW = 1000.0;
        values[pressureIdx] = 1e5 - (domainHeight - globalPos[1])*densityW*g[1];
        values[saturationIdx] = 0.0;
        return values;
    }

    //! returns the temperature in \f$\mathrm{[K]}\f$ in the domain
    Scalar temperature() const
    { return 283.15; /*10°*/ }

    //! sets the pointer to the coupling manager.
    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
    { couplingManagerPtr_ = cm; }

    //! returns reference to the coupling manager.
    const CouplingManager& couplingManager() const
    { return *couplingManagerPtr_; }

private:
    std::shared_ptr<CouplingManager> couplingManagerPtr_;

    Scalar boundaryOverPressure_;
    Scalar boundarySaturation_;
    bool isExercisePartA_;
    bool isExercisePartB_;
    bool isExercisePartC_;
};

} // end namespace Dumux

234
#endif // DUMUX_COURSE_FRACTURESEXERCISE_MATRIX_PROBLEM_HH