2p2ccouplinglocalresidual.hh 13.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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
21
22
 * \brief Extending the TwoPTwoCLocalResidual by the required functions for
 *        a coupled application.
23
24
25
26
27
 */
#ifndef DUMUX_2P2C_COUPLING_LOCAL_RESIDUAL_HH
#define DUMUX_2P2C_COUPLING_LOCAL_RESIDUAL_HH

#include <dumux/implicit/2p2c/2p2clocalresidual.hh>
Natalie Schroeder's avatar
Natalie Schroeder committed
28
#include <dumux/implicit/2p2c/2p2cproperties.hh>
29
30
31
32
33

namespace Dumux
{

/*!
Thomas Fetzer's avatar
Thomas Fetzer committed
34
35
 * \ingroup ImplicitLocalResidual
 * \ingroup TwoPTwoCStokesTwoCModel
36
37
 * \brief Extending the TwoPTwoCLocalResidual by the required functions for
 *        a coupled application.
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
 */
template<class TypeTag>
class TwoPTwoCCouplingLocalResidual : public TwoPTwoCLocalResidual<TypeTag>
{
    typedef TwoPTwoCLocalResidual<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
    typedef typename GridView::IntersectionIterator IntersectionIterator;

    enum { dim = GridView::dimension };
    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum {
        pressureIdx = Indices::pressureIdx
    };
    enum {
        wPhaseIdx = Indices::wPhaseIdx,
        nPhaseIdx = Indices::nPhaseIdx
    };
    enum {
        massBalanceIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
        contiWEqIdx = Indices::contiWEqIdx
    };
    enum {
        wCompIdx = Indices::wCompIdx,
        nCompIdx = Indices::nCompIdx
    };
    enum { phaseIdx = nPhaseIdx }; // index of the phase for the phase flux calculation
    enum { compIdx = wCompIdx}; // index of the component for the phase flux calculation

    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef Dune::BlockVector<Dune::FieldVector<Scalar,1> > ElementFluxVector;

    typedef Dune::FieldVector<Scalar, dim> DimVector;

public:
81
82
83
    /*!
     * \brief Implementation of the boundary evaluation
     */
84
85
86
87
    void evalBoundary_()
    {
        ParentType::evalBoundary_();

88
89
        typedef Dune::ReferenceElements<Scalar, dim> ReferenceElements;
        typedef Dune::ReferenceElement<Scalar, dim> ReferenceElement;
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
        const ReferenceElement &refElement = ReferenceElements::general(this->element_().geometry().type());

        // evaluate Dirichlet-like coupling conditions
        for (int idx = 0; idx < this->fvGeometry_().numScv; idx++)
        {
            // evaluate boundary conditions for the intersections of
            // the current element
            IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
            const IntersectionIterator &endIt = this->gridView_().iend(this->element_());
            for (; isIt != endIt; ++isIt)
            {
                // handle only intersections on the boundary
                if (!isIt->boundary())
                    continue;

                // assemble the boundary for all vertices of the current face
106
107
                const int fIdx = isIt->indexInInside();
                const int numFaceVertices = refElement.size(fIdx, 1, dim);
108
109

                // loop over the single vertices on the current face
110
                for (int faceVertexIdx = 0; faceVertexIdx < numFaceVertices; ++faceVertexIdx)
111
                {
112
113
                    const int boundaryFaceIdx = this->fvGeometry_().boundaryFaceIndex(fIdx, faceVertexIdx);
                    const int vIdx = refElement.subEntity(fIdx, 1, faceVertexIdx, dim);
114
115
                    // only evaluate, if we consider the same face vertex as in the outer
                    // loop over the element vertices
116
                    if (vIdx != idx)
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
                        continue;

                    //for the corner points, the boundary flux across the vertical non-coupling boundary faces
                    //has to be calculated to fulfill the mass balance
                    //convert suddomain intersection into multidomain intersection and check whether it is an outer boundary
                    if(!GridView::Grid::multiDomainIntersection(*isIt).neighbor()
                            && this->boundaryHasMortarCoupling_(this->bcTypes_(idx)))//this->boundaryHasNeumann_(this->bcTypes_(idx)))
                    {
                        const DimVector& globalPos = this->fvGeometry_().subContVol[idx].global;
                        //problem specific function, in problem orientation of interface is known
                        if(this->problem_().isInterfaceCornerPoint(globalPos))
                        {
                            //check whether the second boundary node on the vertical edge is Neumann-Null
                            const int numVertices = refElement.size(dim);
                            bool evalBoundaryFlux = false;
                            for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
                            {
                                for(int i= 0; i < numVertices; i++)
                                {
                                    //if vertex is on boundary and not the coupling vertex: check whether an outflow condition is set
137
                                    if(this->model_().onBoundary(this->element_(), i) && i!=vIdx)
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
                                        if (!this->bcTypes_(i).isOutflow(equationIdx))
                                            evalBoundaryFlux = true;
                                }

                                PrimaryVariables values(0.0);
                                //calculate the actual boundary fluxes and add to residual
                                this->asImp_()->computeFlux(values, boundaryFaceIdx, true /*on boundary*/);
                                if(evalBoundaryFlux)
                                    this->residual_[idx][equationIdx] += values[equationIdx];
                            }

                        }
                    }

                    if (boundaryHasCoupling_(this->bcTypes_(idx)))
                        evalCouplingVertex_(idx);
                }
            }
        }
    }

159
160
161
    /*!
     * \brief Evaluates the phase storage
     */
162
163
164
    Scalar evalPhaseStorage(const int scvIdx) const
    {
        Scalar phaseStorage = computePhaseStorage(scvIdx, false);
165
        Scalar oldPhaseStorage = computePhaseStorage(scvIdx, true);
166
167
168
169
170
171
172
173
174
175
176
177
178
        Valgrind::CheckDefined(phaseStorage);
        Valgrind::CheckDefined(oldPhaseStorage);

        phaseStorage -= oldPhaseStorage;
        phaseStorage *= this->fvGeometry_().subContVol[scvIdx].volume
            / this->problem_().timeManager().timeStepSize()
            * this->curVolVars_(scvIdx).extrusionFactor();
        Valgrind::CheckDefined(phaseStorage);

        return phaseStorage;
    }

    /*!
179
     * \brief Compute storage term of all components within all phases
180
181
182
183
184
185
186
187
188
     */
    Scalar computePhaseStorage(const int scvIdx, bool usePrevSol) const
    {
        const ElementVolumeVariables &elemVolVars = 
            usePrevSol ? this->prevVolVars_() : this->curVolVars_();
        const VolumeVariables &volVars = elemVolVars[scvIdx];

        return volVars.density(phaseIdx)
            * volVars.saturation(phaseIdx)
189
            * volVars.massFraction(phaseIdx, compIdx)
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
            * volVars.porosity();
    }

    /*!
     * \brief Compute the fluxes within the different fluid phases. This is
     *        merely for the computation of flux output.
     */
    void evalPhaseFluxes()
    {
        elementFluxes_.resize(this->fvGeometry_().numScv);
        elementFluxes_ = 0.;

        Scalar flux(0.);

        // calculate the mass flux over the faces and subtract
        // it from the local rates
206
        for (int fIdx = 0; fIdx < this->fvGeometry_().numScvf; fIdx++)
207
208
209
210
        {
            FluxVariables vars(this->problem_(),
                               this->element_(),
                               this->fvGeometry_(),
211
                               fIdx,
212
213
                               this->curVolVars_());

214
215
            int i = this->fvGeometry_().subContVolFace[fIdx].i;
            int j = this->fvGeometry_().subContVolFace[fIdx].j;
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

            const Scalar extrusionFactor =
                (this->curVolVars_(i).extrusionFactor()
                 + this->curVolVars_(j).extrusionFactor())
                / 2;
            flux = computeAdvectivePhaseFlux(vars) +
                computeDiffusivePhaseFlux(vars);
            flux *= extrusionFactor;

            elementFluxes_[i] += flux;
            elementFluxes_[j] -= flux;
        }
    }

    /*!
231
     * \brief Returns the advective fluxes within the different phases.
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
     */
    Scalar computeAdvectivePhaseFlux(const FluxVariables &fluxVars) const
    {
        Scalar advectivePhaseFlux = 0.;
        const Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);

        // data attached to upstream and the downstream vertices
        // of the current phase
        const VolumeVariables &up =
                this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
        const VolumeVariables &dn =
                this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
        if (massUpwindWeight > 0.0)
            // upstream vertex
            advectivePhaseFlux +=
                fluxVars.volumeFlux(phaseIdx)
                * massUpwindWeight
                * up.density(phaseIdx)
250
                * up.massFraction(phaseIdx, compIdx);
251
252
253
254
255
256
        if (massUpwindWeight < 1.0)
            // downstream vertex
            advectivePhaseFlux +=
                fluxVars.volumeFlux(phaseIdx)
                * (1 - massUpwindWeight)
                * dn.density(phaseIdx)
257
                * dn.massFraction(phaseIdx, compIdx);
258
259
260
261
262

        return advectivePhaseFlux;
    }

    /*!
263
     * \brief Returns the diffusive fluxes within the different phases.
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
     */
    Scalar computeDiffusivePhaseFlux(const FluxVariables &fluxVars) const
    {
        // add diffusive flux of gas component in liquid phase
        Scalar diffusivePhaseFlux = fluxVars.moleFractionGrad(phaseIdx)*fluxVars.face().normal;

        return -1.0
                * fluxVars.porousDiffCoeff(phaseIdx)
                * fluxVars.molarDensity(phaseIdx)
                * diffusivePhaseFlux
                * FluidSystem::molarMass(compIdx);
    }

    /*!
     * \brief Set the Dirichlet-like conditions for the coupling
     *        and replace the existing residual
     *
281
     * \param scvIdx Sub control vertex index for the coupling condition
282
283
284
285
286
287
288
289
290
     */
    void evalCouplingVertex_(const int scvIdx)
    {
        const VolumeVariables &volVars = this->curVolVars_()[scvIdx];

        // set pressure as part of the momentum coupling
        if (this->bcTypes_(scvIdx).isCouplingOutflow(massBalanceIdx))
            this->residual_[scvIdx][massBalanceIdx] = volVars.pressure(nPhaseIdx);

Thomas Fetzer's avatar
Thomas Fetzer committed
291
        // set mass fraction;
292
        if (this->bcTypes_(scvIdx).isCouplingOutflow(contiWEqIdx))
293
            this->residual_[scvIdx][contiWEqIdx] = volVars.massFraction(nPhaseIdx, wCompIdx);
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
    }

    /*!
     * \brief Check if one of the boundary conditions is coupling.
     */
    bool boundaryHasCoupling_(const BoundaryTypes& bcTypes) const
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            if (bcTypes.isCouplingInflow(eqIdx) || bcTypes.isCouplingOutflow(eqIdx))
                return true;
        return false;
    }

    /*!
     * \brief Check if one of the boundary conditions is mortar coupling.
     */
    bool boundaryHasMortarCoupling_(const BoundaryTypes& bcTypes) const
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
            if (bcTypes.isMortarCoupling(eqIdx))
                return true;
        return false;
    }

    /*!
     * \brief Check if one of the boundary conditions is Neumann.
     */
    bool boundaryHasNeumann_(const BoundaryTypes& bcTypes) const
    {
        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
        {
            if (bcTypes.isNeumann(eqIdx))
                return true;
        }
        return false;
    }

331
332
333
    /*!
     * \brief Returns the fluxes of the specified sub control volume
     */
334
335
336
337
338
339
340
341
342
    Scalar elementFluxes(const int scvIdx)
    {
        return elementFluxes_[scvIdx];
    }

protected:
    ElementFluxVector elementFluxes_;
};

343
} // namespace Dumux
344

345
#endif // DUMUX_2P2C_COUPLING_LOCAL_RESIDUAL_HH