biominproblem.hh 12 KB
Newer Older
Simon Scholz's avatar
Simon Scholz committed
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
// -*- 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
 *
 * \brief Tutorial problem for a fully coupled TwoPNCMineralization CC Tpfa model.
 */
#ifndef DUMUX_EXERCISE_FOUR_PROBLEM_HH
#define DUMUX_EXERCISE_FOUR_PROBLEM_HH

27
#include <dumux/common/properties.hh>
Simon Scholz's avatar
Simon Scholz committed
28
29
30
31
32
33
34
35
36
37
38
#include <dumux/porousmediumflow/problem.hh>

// TODO: dumux-course-task
// include chemistry file here

namespace Dumux {

/*!
 * \brief Problem biomineralization (MICP) in an experimental setup.
 */
template <class TypeTag>
Katharina Heck's avatar
Katharina Heck committed
39
class BioMinProblem : public PorousMediumFlowProblem<TypeTag>
Simon Scholz's avatar
Simon Scholz committed
40
41
42
{
    using ParentType = PorousMediumFlowProblem<TypeTag>;

43
44
45
46
47
48
49
50
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
51
52
    using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
Simon Scholz's avatar
Simon Scholz committed
53
54
55
56
57
58
59
60
61

    // Grid dimension
    enum
    {
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };

    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
62
63
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
Simon Scholz's avatar
Simon Scholz committed
64
    using Element = typename GridView::template Codim<0>::Entity;
65
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
Simon Scholz's avatar
Simon Scholz committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    // TODO: dumux-course-task
    // set the chemistry TypeTag

    enum
    {
        numComponents = FluidSystem::numComponents,

        pressureIdx = Indices::pressureIdx,
        switchIdx   = Indices::switchIdx, //Saturation

        //Indices of the components
        H2OIdx  = FluidSystem::H2OIdx,
        CO2Idx  = FluidSystem::CO2Idx,
        CaIdx   = FluidSystem::CaIdx,
        UreaIdx = FluidSystem::UreaIdx,

        BiofilmIdx = SolidSystem::BiofilmIdx+numComponents,
        CalciteIdx = SolidSystem::CalciteIdx+numComponents,

        //Index of the primary component of wetting and nonwetting phase
        conti0EqIdx = Indices::conti0EqIdx,

        // Phase State
        liquidPhaseOnly = Indices::firstPhaseOnly,
    };

    /*!
     * \brief The constructor
     *
     * \param fvGridGeometry The finite volume grid geometry
     */
public:
Katharina Heck's avatar
Katharina Heck committed
99
    BioMinProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
Simon Scholz's avatar
Simon Scholz committed
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
    : ParentType(fvGridGeometry)
    {
        name_           = getParam<std::string>("Problem.Name");
        //initial values
        densityW_       = getParam<Scalar>("Initial.InitDensityW");
        initPressure_   = getParam<Scalar>("Initial.InitPressure");
        initBiofilm_    = getParam<Scalar>("Initial.InitBiofilm");

        //injection values
        injVolumeflux_  = getParam<Scalar>("Injection.InjVolumeflux");
        injBioTime_     = getParam<Scalar>("Injection.InjBioTime")*86400;
        injCO2_         = getParam<Scalar>("Injection.InjCO2");
        concCa_         = getParam<Scalar>("Injection.ConcCa");
        concUrea_       = getParam<Scalar>("Injection.ConcUrea");

115
        unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethod::box ? dim : 0;
Simon Scholz's avatar
Simon Scholz committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
        Kxx_.resize(fvGridGeometry->gridView().size(codim));
        Kyy_.resize(fvGridGeometry->gridView().size(codim));

        //initialize the fluidsystem
        FluidSystem::init(/*startTemp=*/288,
                          /*endTemp=*/ 303,
                          /*tempSteps=*/5,
                          /*startPressure=*/1e4,
                          /*endPressure=*/1e6,
                          /*pressureSteps=*/500);
    }

    void setTime(Scalar time)
    {
        time_ = time;
    }

    /*!
     * \name Problem parameters
     */

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

    /*!
     * \brief Returns the temperature within the domain.
     *
     * This problem assumes a temperature of 300 degrees Kelvin.
     */
    Scalar temperature() const
    {
        return 300; // [K]
    }

    // \}

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

    /*!
     * \brief Specifies which kind of boundary condition should be
     *        used for which equation on a given boundary segment.
     */
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes bcTypes;

        // we don't set any BCs for the solid phases
171
        Scalar xmax = this->gridGeometry().bBoxMax()[0];
Simon Scholz's avatar
Simon Scholz committed
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
        // set all Neumann
        bcTypes.setAllNeumann();

        // set right boundary to Dirichlet
        if(globalPos[0] > xmax - eps_)
            bcTypes.setAllDirichlet();

        return bcTypes;
    }

    /*!
     * \brief Evaluates the boundary conditions for a Dirichlet
     *        boundary segment
     *
     * \param globalPos The global position
     */
    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
    {
        // recycle initialAtPos() for Dirichlet values
        return initialAtPos(globalPos);
    }

    /*!
     * \brief Evaluate the boundary conditions for a neumann
     *        boundary segment.
     *
     * \param values Stores the Neumann values for the conservation equations in
     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
     * \param globalPos The position of the integration point of the boundary segment.
     *
     * For this method, the \a values parameter stores the mass flux
     * in normal direction of each phase. Negative values mean influx.
     */
    NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
    {
        NumEqVector values(0.0);

        Scalar waterFlux = injVolumeflux_; // divide by area if area not 1! [m/s]

        // Set values for Ca + urea injection above aquitard.
        // Use negative values for injection.
        if(globalPos[0] < eps_
214
           && globalPos[1] > 11.0 - eps_
Simon Scholz's avatar
Simon Scholz committed
215
216
217
218
219
220
221
222
223
224
225
           && globalPos[1] < 12.0 + eps_
           && time_ < injBioTime_)
        {
            values[conti0EqIdx + H2OIdx]  = - waterFlux * densityW_ / FluidSystem::molarMass(H2OIdx);
            values[conti0EqIdx + CaIdx]   = - waterFlux * concCa_ / FluidSystem::molarMass(CaIdx);
            values[conti0EqIdx + UreaIdx] = - waterFlux * concUrea_ / FluidSystem::molarMass(UreaIdx);
        }
        // TODO: dumux-course-task
        // Set CO2 injection below aquitard after the injBioTime
        else
        {
Timo Koch's avatar
Timo Koch committed
226
            // Scalar gasFlux = injCO2_; // divide by area if area not 1! [m/s]
Simon Scholz's avatar
Simon Scholz committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
            values = 0.0; //mol/m²/s
        }

        return values;
    }

    /*!
     * \brief Evaluate the initial value for a control volume.
     *
     * \param globalPos The global position
     *
     * For this method, the \a values parameter stores primary
     * variables.
     */
    PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
    {
        PrimaryVariables priVars(0.0);
        priVars.setState(liquidPhaseOnly);

246
        Scalar hmax = this->gridGeometry().bBoxMax()[1];
Simon Scholz's avatar
Simon Scholz committed
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
305
306
307
308
309
310
311
312
313
        priVars[pressureIdx] = initPressure_ - (globalPos[1] - hmax)*densityW_*9.81;
        priVars[switchIdx]   = 0.0;
        priVars[CaIdx]       = 0.0;
        priVars[UreaIdx]     = 0.0;
        priVars[CalciteIdx]  = 0.0;
        priVars[BiofilmIdx]  = 0.0;

        // set Biofilm presence only in aquitard zone and below
        if (globalPos[1] < 10 + eps_)
        {
            priVars[switchIdx]  = 0.0;
            priVars[BiofilmIdx] = initBiofilm_; // [m^3/m^3]
        }

        return priVars;
    }

    /*!
     * \brief Evaluate the source term for all phases within a given
     *        sub-control-volume.
     *
     * This is the method for the case where the source term is
     * potentially solution dependent and requires some quantities that
     * are specific to the fully-implicit method.
     *
     * \param element The finite element
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
     * \param scv The sub control volume
     *
     * For this method, the return parameter stores the conserved quantity rate
     * generated or annihilate per volume unit. Positive values mean
     * that the conserved quantity is created, negative ones mean that it vanishes.
     * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
     */
    NumEqVector source(const Element &element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& elemVolVars,
                       const SubControlVolume &scv) const
    {
        NumEqVector source(0.0);

        // TODO: dumux-course-task
        // set Chemistry
        // set VolumeVariables
        // call reactionSource() in chemistry

        return source;
    }


    const std::vector<Scalar>& getKxx()
    {
        return Kxx_;
    }

    const std::vector<Scalar>& getKyy()
    {
        return Kyy_;
    }

    /*!
     * \brief Adds additional VTK output data to the VTKWriter.
     * Function is called by the output module on every write.
     */
    void updateVtkOutput(const SolutionVector& curSol)
    {
314
        for (const auto& element : elements(this->gridGeometry().gridView()))
Simon Scholz's avatar
Simon Scholz committed
315
        {
316
            const auto elemSol = elementSolution(element, curSol, this->gridGeometry());
Simon Scholz's avatar
Simon Scholz committed
317

318
            auto fvGeometry = localView(this->gridGeometry());
Simon Scholz's avatar
Simon Scholz committed
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
            fvGeometry.bindElement(element);

            for (auto&& scv : scvs(fvGeometry))
            {
                VolumeVariables volVars;
                volVars.update(elemSol, *this, element, scv);
                const auto dofIdxGlobal = scv.dofIndex();
                Kxx_[dofIdxGlobal] = volVars.permeability()[0][0];
                Kyy_[dofIdxGlobal] = volVars.permeability()[1][1];
            }
        }
    }

private:
    static constexpr Scalar eps_ = 1e-6;

    Scalar densityW_;
    Scalar initPressure_;
    Scalar initBiofilm_;

    Scalar injVolumeflux_;
    Scalar injBioTime_;
    Scalar injCO2_;
    Scalar concCa_;
    Scalar concUrea_;

    std::vector<double> Kxx_;
    std::vector<double> Kyy_;
    std::string name_;
    Scalar time_;
};

} //end namespace

#endif