incompressiblelocalresidual.hh 30 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
 * \ingroup TwoPModel
22
23
24
25
26
27
 * \brief Element-wise calculation of the residual and its derivatives
 *        for a two-phase, incompressible test problem.
 */
#ifndef DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH
#define DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH

Timo Koch's avatar
Timo Koch committed
28
29
30
#include <cmath>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
Timo Koch's avatar
Timo Koch committed
31
#include <dumux/discretization/methods.hh>
32
#include <dumux/discretization/elementsolution.hh>
33
34
#include <dumux/porousmediumflow/immiscible/localresidual.hh>

35
namespace Dumux {
36

37
38
39
40
41
/*!
 * \ingroup TwoPModel
 * \brief Element-wise calculation of the residual and its derivatives
 *        for a two-phase, incompressible test problem.
 */
42
43
44
45
46
47
48
template<class TypeTag>
class TwoPIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag>
{
    using ParentType = ImmiscibleLocalResidual<TypeTag>;
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
49
    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
50
    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
51
52
    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
53
54
55
    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
    using FVElementGeometry = typename FVGridGeometry::LocalView;
56
57
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
    using Element = typename GridView::template Codim<0>::Entity;
    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
    // first index for the mass balance
    enum
    {
        contiWEqIdx = Indices::contiWEqIdx,
        contiNEqIdx = Indices::contiNEqIdx,

        pressureIdx = Indices::pressureIdx,
        saturationIdx = Indices::saturationIdx
    };

public:
    using ParentType::ParentType;

75
76
77
78
79
80
81
82
83
84
85
86
87
    /*!
     * \brief Add storage derivatives for wetting and non-wetting phase
     *
     * Compute storage derivatives for the wetting and the non-wetting phase with respect to \f$p_w\f$
     * and \f$S_n\f$.
     *
     * \param partialDerivatives The partial derivatives
     * \param problem The problem
     * \param element The element
     * \param fvGeometry The finite volume element geometry
     * \param curVolVars The current volume variables
     * \param scv The sub control volume
     */
88
89
90
91
92
    template<class PartialDerivativeMatrix>
    void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
93
94
                               const VolumeVariables& curVolVars,
                               const SubControlVolume& scv) const
95
    {
Dennis Gläser's avatar
Dennis Gläser committed
96
97
98
99
100
        static_assert(!FluidSystem::isCompressible(FluidSystem::wPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");
        static_assert(!FluidSystem::isCompressible(FluidSystem::nPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");

101
102
103
104
105
        // we know that these values are constant throughout the simulation
        static const auto phi = curVolVars.porosity();
        static const auto phi_rho_w = phi*curVolVars.density(FluidSystem::wPhaseIdx);
        static const auto phi_rho_n = phi*curVolVars.density(FluidSystem::nPhaseIdx);

106
        const auto volume = scv.volume();
107
108
109
110
111
112
113
114
115
116
117

        // partial derivative of wetting phase storage term w.r.t. p_w
        partialDerivatives[contiWEqIdx][pressureIdx] += 0.0;
        // partial derivative of wetting phase storage term w.r.t. S_n
        partialDerivatives[contiWEqIdx][saturationIdx] -= volume*phi_rho_w/this->timeLoop().timeStepSize();
        // partial derivative of non-wetting phase storage term w.r.t. p_w
        partialDerivatives[contiNEqIdx][pressureIdx] += 0.0;
        // partial derivative of non-wetting phase storage term w.r.t. S_n
        partialDerivatives[contiNEqIdx][saturationIdx] += volume*phi_rho_n/this->timeLoop().timeStepSize();
    }

118
119
120
121
122
123
124
125
126
127
    /*!
     * \brief Add source derivatives for wetting and non-wetting phase
     *
     * \param partialDerivatives The partial derivatives
     * \param problem The problem
     * \param element The element
     * \param fvGeometry The finite volume element geometry
     * \param curVolVars The current volume variables
     * \param scv The sub control volume
     */
128
129
130
131
132
    template<class PartialDerivativeMatrix>
    void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
                              const Problem& problem,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
133
                              const VolumeVariables& curVolVars,
Dennis Gläser's avatar
Dennis Gläser committed
134
135
                              const SubControlVolume& scv) const
    { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ }
136

137
    /*!
138
     * \brief Add flux derivatives for wetting and non-wetting phase for cell-centered FVM using TPFA
139
140
141
142
143
144
145
146
147
148
149
     *
     * Compute derivatives for the wetting and the non-wetting phase flux with respect to \f$p_w\f$
     * and \f$S_n\f$.
     *
     * \param partialDerivatives The partial derivatives
     * \param problem The problem
     * \param element The element
     * \param fvGeometry The finite volume element geometry
     * \param curVolVars The current volume variables
     * \param scv The sub control volume
     */
150
    template<class PartialDerivativeMatrices, class T = TypeTag>
151
    std::enable_if_t<GET_PROP_TYPE(T, FVGridGeometry)::discMethod == DiscretizationMethod::cctpfa, void>
152
153
154
155
156
157
158
    addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                       const Problem& problem,
                       const Element& element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& curElemVolVars,
                       const ElementFluxVariablesCache& elemFluxVarsCache,
                       const SubControlVolumeFace& scvf) const
159
    {
Dennis Gläser's avatar
Dennis Gläser committed
160
161
162
163
164
165
166
167
168
        static_assert(!FluidSystem::isCompressible(FluidSystem::wPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");
        static_assert(!FluidSystem::isCompressible(FluidSystem::nPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");
        static_assert(FluidSystem::viscosityIsConstant(FluidSystem::wPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");
        static_assert(FluidSystem::viscosityIsConstant(FluidSystem::nPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");

169
170
171
172
        using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);

        // evaluate the current wetting phase Darcy flux and resulting upwind weights
173
        static const Scalar upwindWeight = getParam<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
Dennis Gläser's avatar
Dennis Gläser committed
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
        const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
                                                scvf, FluidSystem::wPhaseIdx, elemFluxVarsCache);
        const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
                                                scvf, FluidSystem::nPhaseIdx, elemFluxVarsCache);
        const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_w = 1.0 - insideWeight_w;
        const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_n = 1.0 - insideWeight_n;

        // get references to the two participating vol vars & parameters
        const auto insideScvIdx = scvf.insideScvIdx();
        const auto outsideScvIdx = scvf.outsideScvIdx();
        const auto outsideElement = fvGeometry.fvGridGeometry().element(outsideScvIdx);
        const auto& insideScv = fvGeometry.scv(insideScvIdx);
        const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
        const auto& insideVolVars = curElemVolVars[insideScvIdx];
        const auto& outsideVolVars = curElemVolVars[outsideScvIdx];
        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
                                                                                     insideScv,
193
                                                                                     elementSolution<FVElementGeometry>(insideVolVars.priVars()));
Dennis Gläser's avatar
Dennis Gläser committed
194
195
        const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(outsideElement,
                                                                                      outsideScv,
196
                                                                                      elementSolution<FVElementGeometry>(outsideVolVars.priVars()));
Dennis Gläser's avatar
Dennis Gläser committed
197
198
199
200

        // get references to the two participating derivative matrices
        auto& dI_dI = derivativeMatrices[insideScvIdx];
        auto& dI_dJ = derivativeMatrices[outsideScvIdx];
201
202
203
204
205
206
207
208
209
210
211
212
213
214

        // some quantities to be reused (rho & mu are constant and thus equal for all cells)
        static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
        static const auto rho_n = insideVolVars.density(FluidSystem::nPhaseIdx);
        static const auto rhow_muw = rho_w/insideVolVars.viscosity(FluidSystem::wPhaseIdx);
        static const auto rhon_mun = rho_n/insideVolVars.viscosity(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(FluidSystem::nPhaseIdx);

        // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
        const auto insideSw = insideVolVars.saturation(FluidSystem::wPhaseIdx);
        const auto outsideSw = outsideVolVars.saturation(FluidSystem::wPhaseIdx);
Dennis Gläser's avatar
Dennis Gläser committed
215
216
217
218
219
220
        const auto dKrw_dSn_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
        const auto dKrw_dSn_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
        const auto dKrn_dSn_inside = MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
        const auto dKrn_dSn_outside = MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
        const auto dpc_dSn_inside = MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
        const auto dpc_dSn_outside = MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
221

Dennis Gläser's avatar
Dennis Gläser committed
222
        const auto tij = elemFluxVarsCache[scvf].advectionTij();
223

Dennis Gläser's avatar
Dennis Gläser committed
224
225
226
227
228
229
230
        // precalculate values
        const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w;
        const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n;
        const auto rho_mu_flux_w = rhow_muw*flux_w;
        const auto rho_mu_flux_n = rhon_mun*flux_n;
        const auto tij_up_w = tij*up_w;
        const auto tij_up_n = tij*up_n;
231
232

        // partial derivative of the wetting phase flux w.r.t. p_w
Dennis Gläser's avatar
Dennis Gläser committed
233
234
        dI_dI[contiWEqIdx][pressureIdx] += tij_up_w;
        dI_dJ[contiWEqIdx][pressureIdx] -= tij_up_w;
235
236

        // partial derivative of the wetting phase flux w.r.t. S_n
Dennis Gläser's avatar
Dennis Gläser committed
237
238
        dI_dI[contiWEqIdx][saturationIdx] -= rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
        dI_dJ[contiWEqIdx][saturationIdx] -= rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
239
240

        // partial derivative of the non-wetting phase flux w.r.t. p_w
Dennis Gläser's avatar
Dennis Gläser committed
241
242
        dI_dI[contiNEqIdx][pressureIdx] += tij_up_n;
        dI_dJ[contiNEqIdx][pressureIdx] -= tij_up_n;
243
244

        // partial derivative of the non-wetting phase flux w.r.t. S_n (relative permeability derivative contribution)
Dennis Gläser's avatar
Dennis Gläser committed
245
246
        dI_dI[contiNEqIdx][saturationIdx] -= rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
        dI_dJ[contiNEqIdx][saturationIdx] -= rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
247
248

        // partial derivative of the non-wetting phase flux w.r.t. S_n (capillary pressure derivative contribution)
Dennis Gläser's avatar
Dennis Gläser committed
249
250
251
252
        dI_dI[contiNEqIdx][saturationIdx] -= tij_up_n*dpc_dSn_inside;
        dI_dJ[contiNEqIdx][saturationIdx] += tij_up_n*dpc_dSn_outside;
    }

253
254
255
256
257
258
259
260
261
262
263
264
265
    /*!
     * \brief Add flux derivatives for wetting and non-wetting phase for box method
     *
     * Compute derivatives for the wetting and the non-wetting phase flux with respect to \f$p_w\f$
     * and \f$S_n\f$.
     *
     * \param partialDerivatives The partial derivatives
     * \param problem The problem
     * \param element The element
     * \param fvGeometry The finite volume element geometry
     * \param curVolVars The current volume variables
     * \param scv The sub control volume
     */
Dennis Gläser's avatar
Dennis Gläser committed
266
    template<class JacobianMatrix, class T = TypeTag>
267
    std::enable_if_t<GET_PROP_TYPE(T, FVGridGeometry)::discMethod == DiscretizationMethod::box, void>
Dennis Gläser's avatar
Dennis Gläser committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
    addFluxDerivatives(JacobianMatrix& A,
                       const Problem& problem,
                       const Element& element,
                       const FVElementGeometry& fvGeometry,
                       const ElementVolumeVariables& curElemVolVars,
                       const ElementFluxVariablesCache& elemFluxVarsCache,
                       const SubControlVolumeFace& scvf) const
    {
        static_assert(!FluidSystem::isCompressible(FluidSystem::wPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");
        static_assert(!FluidSystem::isCompressible(FluidSystem::nPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!");
        static_assert(FluidSystem::viscosityIsConstant(FluidSystem::wPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");
        static_assert(FluidSystem::viscosityIsConstant(FluidSystem::nPhaseIdx),
                      "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!");

        using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);

        // evaluate the current wetting phase Darcy flux and resulting upwind weights
289
        static const Scalar upwindWeight = getParam<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
Dennis Gläser's avatar
Dennis Gläser committed
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
        const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
                                                scvf, FluidSystem::wPhaseIdx, elemFluxVarsCache);
        const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
                                                scvf, FluidSystem::nPhaseIdx, elemFluxVarsCache);
        const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_w = 1.0 - insideWeight_w;
        const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_n = 1.0 - insideWeight_n;

        // get references to the two participating vol vars & parameters
        const auto insideScvIdx = scvf.insideScvIdx();
        const auto outsideScvIdx = scvf.outsideScvIdx();
        const auto& insideScv = fvGeometry.scv(insideScvIdx);
        const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
        const auto& insideVolVars = curElemVolVars[insideScv];
        const auto& outsideVolVars = curElemVolVars[outsideScv];

307
308
        const auto elemSol = elementSolution(element, curElemVolVars, fvGeometry);

Dennis Gläser's avatar
Dennis Gläser committed
309
310
311
312
313
314
315
316
317
318
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
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element, insideScv, elemSol);
        const auto& outsideMaterialParams = problem.spatialParams().materialLawParams(element, outsideScv, elemSol);

        // some quantities to be reused (rho & mu are constant and thus equal for all cells)
        static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
        static const auto rho_n = insideVolVars.density(FluidSystem::nPhaseIdx);
        static const auto rhow_muw = rho_w/insideVolVars.viscosity(FluidSystem::wPhaseIdx);
        static const auto rhon_mun = rho_n/insideVolVars.viscosity(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(FluidSystem::nPhaseIdx);

        // let the Law for the advective fluxes calculate the transmissibilities
        const auto ti = AdvectionType::calculateTransmissibilities(problem,
                                                                   element,
                                                                   fvGeometry,
                                                                   curElemVolVars,
                                                                   scvf,
                                                                   elemFluxVarsCache[scvf]);

        // get the rows of the jacobian matrix for the inside/outside scv
        auto& dI_dJ_inside = A[insideScv.dofIndex()];
        auto& dI_dJ_outside = A[outsideScv.dofIndex()];

        // precalculate values
        const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w;
        const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n;
        const auto rho_mu_flux_w = rhow_muw*flux_w;
        const auto rho_mu_flux_n = rhon_mun*flux_n;

        // add the partial derivatives w.r.t all scvs in the element
        for (const auto& scvJ : scvs(fvGeometry))
        {
            const auto globalJ = scvJ.dofIndex();
            const auto localJ = scvJ.indexInElement();

            // the transmissibily associated with the scvJ
            const auto tj = ti[localJ];

            // partial derivative of the wetting phase flux w.r.t. p_w
            const auto tj_up_w = tj*up_w;
            dI_dJ_inside[globalJ][contiWEqIdx][pressureIdx] += tj_up_w;
            dI_dJ_outside[globalJ][contiWEqIdx][pressureIdx] -= tj_up_w;

            // partial derivative of the non-wetting phase flux w.r.t. p_w
            const auto tj_up_n = tj*up_n;
            dI_dJ_inside[globalJ][contiNEqIdx][pressureIdx] += tj_up_n;
            dI_dJ_outside[globalJ][contiNEqIdx][pressureIdx] -= tj_up_n;

            // partial derivatives w.r.t. S_n (are the negative of those w.r.t sw)
            // relative permeability contributions only for inside/outside
            if (localJ == insideScvIdx)
            {
                // partial derivative of the wetting phase flux w.r.t. S_n
                const auto insideSw = insideVolVars.saturation(FluidSystem::wPhaseIdx);
                const auto dKrw_dSn_inside = MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
                const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w;
                dI_dJ_inside[globalJ][contiWEqIdx][saturationIdx] -= dFluxW_dSnJ;
                dI_dJ_outside[globalJ][contiWEqIdx][saturationIdx] += dFluxW_dSnJ;

                // partial derivative of the non-wetting phase flux w.r.t. S_n (k_rn contribution)
                const auto dKrn_dSn_inside = MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
                const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n;
                dI_dJ_inside[globalJ][contiNEqIdx][saturationIdx] -= dFluxN_dSnJ_krn;
                dI_dJ_outside[globalJ][contiWEqIdx][saturationIdx] += dFluxN_dSnJ_krn;

                // partial derivative of the non-wetting phase flux w.r.t. S_n (p_c contribution)
                const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);
                dI_dJ_inside[globalJ][contiNEqIdx][saturationIdx] -= dFluxN_dSnJ_pc;
                dI_dJ_outside[globalJ][contiNEqIdx][saturationIdx] += dFluxN_dSnJ_pc;
            }
            else if (localJ == outsideScvIdx)
            {
                // see comments for (globalJ == insideScvIdx)
                const auto outsideSw = outsideVolVars.saturation(FluidSystem::wPhaseIdx);
                const auto dKrw_dSn_outside = MaterialLaw::dkrw_dsw(outsideMaterialParams, outsideSw);
                const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w;
                dI_dJ_inside[globalJ][contiWEqIdx][saturationIdx] -= dFluxW_dSnJ;
                dI_dJ_outside[globalJ][contiWEqIdx][saturationIdx] += dFluxW_dSnJ;

                const auto dKrn_dSn_outside = MaterialLaw::dkrn_dsw(outsideMaterialParams, outsideSw);
                const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n;
                dI_dJ_inside[globalJ][contiNEqIdx][saturationIdx] -= dFluxN_dSnJ_krn;
                dI_dJ_outside[globalJ][contiNEqIdx][saturationIdx] += dFluxN_dSnJ_krn;

                const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(outsideMaterialParams, outsideSw);
                dI_dJ_inside[globalJ][contiNEqIdx][saturationIdx] -= dFluxN_dSnJ_pc;
                dI_dJ_outside[globalJ][contiNEqIdx][saturationIdx] += dFluxN_dSnJ_pc;
            }
            else
            {
                const auto& paramsJ = problem.spatialParams().materialLawParams(element, scvJ, elemSol);
                const auto swJ = curElemVolVars[scvJ].saturation(FluidSystem::wPhaseIdx);
                const auto dFluxN_dSnJ_pc = tj_up_n*MaterialLaw::dpc_dsw(paramsJ, swJ);
                dI_dJ_inside[globalJ][contiNEqIdx][saturationIdx] -= dFluxN_dSnJ_pc;
                dI_dJ_outside[globalJ][contiNEqIdx][saturationIdx] += dFluxN_dSnJ_pc;
            }
        }
408
409
    }

410
411
412
413
414
415
416
417
418
419
420
421
422
    /*!
     * \brief Add cell-centered Dirichlet flux derivatives for wetting and non-wetting phase
     *
     * Compute derivatives for the wetting and the non-wetting phase flux with respect to \f$p_w\f$
     * and \f$S_n\f$.
     *
     * \param derivativeMatrices The matrices containing the derivatives
     * \param problem The problem
     * \param element The element
     * \param curElemVolVars The current element volume variables
     * \param elemFluxVarsCache The element flux variables cache
     * \param scvf The sub control volume face
     */
423
    template<class PartialDerivativeMatrices>
424
425
426
427
428
429
430
    void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                       const Problem& problem,
                                       const Element& element,
                                       const FVElementGeometry& fvGeometry,
                                       const ElementVolumeVariables& curElemVolVars,
                                       const ElementFluxVariablesCache& elemFluxVarsCache,
                                       const SubControlVolumeFace& scvf) const
431
432
433
434
435
    {
        using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
        using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);

        // evaluate the current wetting phase Darcy flux and resulting upwind weights
436
        static const Scalar upwindWeight = getParam<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
Dennis Gläser's avatar
Dennis Gläser committed
437
        const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
438
                                                     scvf, FluidSystem::wPhaseIdx, elemFluxVarsCache);
Dennis Gläser's avatar
Dennis Gläser committed
439
        const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars,
440
                                                        scvf, FluidSystem::nPhaseIdx, elemFluxVarsCache);
Dennis Gläser's avatar
Dennis Gläser committed
441
442
443
444
445
446
447
448
449
        const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_w = 1.0 - insideWeight_w;
        const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight;
        const auto outsideWeight_n = 1.0 - insideWeight_n;

        // get references to the two participating vol vars & parameters
        const auto insideScvIdx = scvf.insideScvIdx();
        const auto& insideScv = fvGeometry.scv(insideScvIdx);
        const auto& insideVolVars = curElemVolVars[insideScvIdx];
450
        const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()];
Dennis Gläser's avatar
Dennis Gläser committed
451
452
        const auto& insideMaterialParams = problem.spatialParams().materialLawParams(element,
                                                                                     insideScv,
453
                                                                                     elementSolution<FVElementGeometry>(insideVolVars.priVars()));
454
455
456
457
458
459
460
461
462
463
464

        // some quantities to be reused (rho & mu are constant and thus equal for all cells)
        static const auto rho_w = insideVolVars.density(FluidSystem::wPhaseIdx);
        static const auto rho_n = insideVolVars.density(FluidSystem::nPhaseIdx);
        static const auto rhow_muw = rho_w/insideVolVars.viscosity(FluidSystem::wPhaseIdx);
        static const auto rhon_mun = rho_n/insideVolVars.viscosity(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(FluidSystem::nPhaseIdx);
        const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(FluidSystem::wPhaseIdx);
        const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(FluidSystem::nPhaseIdx);

Dennis Gläser's avatar
Dennis Gläser committed
465
466
        // get reference to the inside derivative matrix
        auto& dI_dI = derivativeMatrices[insideScvIdx];
467
468
469
470
471
472
473
474
475

        // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw
        const auto insideSw = insideVolVars.saturation(FluidSystem::wPhaseIdx);
        const auto dKrw_dSn_inside = -1.0*MaterialLaw::dkrw_dsw(insideMaterialParams, insideSw);
        const auto dKrn_dSn_inside = -1.0*MaterialLaw::dkrn_dsw(insideMaterialParams, insideSw);
        const auto dpc_dSn_inside = -1.0*MaterialLaw::dpc_dsw(insideMaterialParams, insideSw);

        const auto tij = elemFluxVarsCache[scvf].advectionTij();
        // partial derivative of the wetting phase flux w.r.t. p_w
Dennis Gläser's avatar
Dennis Gläser committed
476
477
        const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w;
        dI_dI[contiWEqIdx][pressureIdx] += tij*up_w;
478
479

        // partial derivative of the wetting phase flux w.r.t. S_n
Dennis Gläser's avatar
Dennis Gläser committed
480
        dI_dI[contiWEqIdx][saturationIdx] += rhow_muw*flux_w*dKrw_dSn_inside*insideWeight_w;
481
482

        // partial derivative of the non-wetting phase flux w.r.t. p_w
Dennis Gläser's avatar
Dennis Gläser committed
483
484
        const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n;
        dI_dI[contiNEqIdx][pressureIdx] += tij*up_n;
485
486

        // partial derivative of the non-wetting phase flux w.r.t. S_n (relative permeability derivative contribution)
Dennis Gläser's avatar
Dennis Gläser committed
487
        dI_dI[contiNEqIdx][saturationIdx] += rhon_mun*flux_n*dKrn_dSn_inside*insideWeight_n;
488
489

        // partial derivative of the non-wetting phase flux w.r.t. S_n (capillary pressure derivative contribution)
Dennis Gläser's avatar
Dennis Gläser committed
490
        dI_dI[contiNEqIdx][saturationIdx] += tij*dpc_dSn_inside*up_n;
491
492
    }

493
494
495
496
497
498
499
500
501
502
    /*!
     * \brief Add Robin flux derivatives for wetting and non-wetting phase
     *
     * \param derivativeMatrices The matrices containing the derivatives
     * \param problem The problem
     * \param element The element
     * \param curElemVolVars The current element volume variables
     * \param elemFluxVarsCache The element flux variables cache
     * \param scvf The sub control volume face
     */
503
504
505
506
507
508
509
510
    template<class PartialDerivativeMatrices>
    void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                 const Problem& problem,
                                 const Element& element,
                                 const FVElementGeometry& fvGeometry,
                                 const ElementVolumeVariables& curElemVolVars,
                                 const ElementFluxVariablesCache& elemFluxVarsCache,
                                 const SubControlVolumeFace& scvf) const
Dennis Gläser's avatar
Dennis Gläser committed
511
    { /* TODO maybe forward to problem for the user to implement the Robin derivatives?*/ }
512
513
514
515
516
};

} // end namespace Dumux

#endif