cclocalassembler.hh 26.6 KB
Newer Older
1
2
3
4
5
6
7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
23
 * \ingroup Assembly
 * \ingroup CCDiscretization
 * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
24
 */
25
26
#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
#define DUMUX_CC_LOCAL_ASSEMBLER_HH
27

28
#include <dune/common/reservedvector.hh>
Timo Koch's avatar
Timo Koch committed
29
#include <dune/grid/common/gridenums.hh> // for GhostEntity
30
31
#include <dune/istl/matrixindexset.hh>

32
#include <dumux/common/reservedblockvector.hh>
33
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
34
#include <dumux/common/parameters.hh>
35
#include <dumux/common/numericdifferentiation.hh>
36
#include <dumux/assembly/numericepsilon.hh>
37
#include <dumux/assembly/diffmethod.hh>
Kilian Weishaupt's avatar
Kilian Weishaupt committed
38
#include <dumux/assembly/fvlocalassemblerbase.hh>
39
40
#include <dumux/assembly/entitycolor.hh>
#include <dumux/assembly/partialreassembler.hh>
41
#include <dumux/discretization/fluxstencil.hh>
42
#include <dumux/discretization/cellcentered/elementsolution.hh>
43

44
45
46
namespace Dumux {

/*!
47
48
 * \ingroup Assembly
 * \ingroup CCDiscretization
Kilian Weishaupt's avatar
Kilian Weishaupt committed
49
50
51
52
 * \brief A base class for all local cell-centered assemblers
 * \tparam TypeTag The TypeTag
 * \tparam Assembler The assembler type
 * \tparam Implementation The actual implementation
53
 * \tparam implicit Specifies whether the time discretization is implicit or not (i.e. explicit)
54
 */
55
56
template<class TypeTag, class Assembler, class Implementation, bool implicit>
class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
57
{
58
    using ParentType = FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>;
59
60
61
62
63
64
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
65

66
67
    static_assert(!Assembler::Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");

68
69
public:

70
71
    using ParentType::ParentType;

72
73
74
75
    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix. The element residual is written into the right hand side.
     */
76
77
78
    template <class PartialReassembler = DefaultPartialReassembler>
    void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables,
                                     const PartialReassembler* partialReassembler)
79
    {
80
        this->asImp_().bindLocalViews();
81
        const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
82
83
84
85
86
87
88
89
90
        if (partialReassembler
            && partialReassembler->elementColor(globalI) == EntityColor::green)
        {
            res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
        }
        else
        {
            res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
        }
91
92
93
94
95
96
    }

    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix.
     */
97
    void assembleJacobian(JacobianMatrix& jac, GridVariables& gridVariables)
98
    {
99
100
        this->asImp_().bindLocalViews();
        this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
101
102
103
104
105
    }

    /*!
     * \brief Assemble the residual only
     */
106
    void assembleResidual(SolutionVector& res)
107
    {
108
        this->asImp_().bindLocalViews();
109
        const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
110
        res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
111
    }
112
113
114
115
116
117
};

/*!
 * \ingroup Assembly
 * \ingroup CCDiscretization
 * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
Kilian Weishaupt's avatar
Kilian Weishaupt committed
118
 * \tparam TypeTag The TypeTag
119
 * \tparam diffMethod The differentiation method to residual compute derivatives
Kilian Weishaupt's avatar
Kilian Weishaupt committed
120
 * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
121
 */
122
template<class TypeTag, class Assembler, DiffMethod diffMethod = DiffMethod::numeric, bool implicit = true>
123
class CCLocalAssembler;
124

125
126
127
128
129
130
131
/*!
 * \ingroup Assembly
 * \ingroup CCDiscretization
 * \brief Cell-centered scheme local assembler using numeric differentiation and implicit time discretization
 */
template<class TypeTag, class Assembler>
class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
132
: public CCLocalAssemblerBase<TypeTag, Assembler,
133
                              CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>, true >
134
{
135
    using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>;
136
    using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
137
138
139
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using Element = typename GetPropType<TypeTag, Properties::GridView>::template Codim<0>::Entity;
140
141
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using FVElementGeometry = typename GridGeometry::LocalView;
142
143
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
144

145
146
    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
    enum { dim = GetPropType<TypeTag, Properties::GridView>::dimension };
147

148
    using FluxStencil = Dumux::FluxStencil<FVElementGeometry>;
149
    static constexpr int maxElementStencilSize = GridGeometry::maxElementStencilSize;
150
    static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
151

152
public:
153

154
    using ParentType::ParentType;
155

156
157
158
159
160
161
    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix.
     *
     * \return The element residual at the current solution.
     */
162
    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
163
    {
164
165
166
167
168
169
        //////////////////////////////////////////////////////////////////////////////////////////////////
        // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
        // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
        // actual element. In the actual element we evaluate the derivative of the entire residual.     //
        //////////////////////////////////////////////////////////////////////////////////////////////////

170
171
172
        // get some aliases for convenience
        const auto& element = this->element();
        const auto& fvGeometry = this->fvGeometry();
173
        const auto& gridGeometry = this->assembler().gridGeometry();
174
175
176
        auto&& curElemVolVars = this->curElemVolVars();
        auto&& elemFluxVarsCache = this->elemFluxVarsCache();

177
        // get stencil informations
178
179
        const auto globalI = gridGeometry.elementMapper().index(element);
        const auto& connectivityMap = gridGeometry.connectivityMap();
180
        const auto numNeighbors = connectivityMap[globalI].size();
181
182

        // container to store the neighboring elements
183
        Dune::ReservedVector<Element, maxElementStencilSize> neighborElements;
184
        neighborElements.resize(numNeighbors);
185

186
        // assemble the undeflected residual
187
        using Residuals = ReservedBlockVector<NumEqVector, maxElementStencilSize>;
188
        Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
189
        origResiduals[0] = this->evalLocalResidual()[0];
190

191
        // lambda for convenient evaluation of the fluxes across scvfs in the neighbors
192
193
        // if the neighbor is a ghost we don't want to add anything to their residual
        // so we return 0 and omit computing the flux
194
195
        auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
        {
196
197
198
199
200
201
202
203
            if (neighbor.partitionType() == Dune::GhostEntity)
                return NumEqVector(0.0);
            else
                return this->localResidual().evalFlux(this->problem(),
                                                      neighbor,
                                                      this->fvGeometry(),
                                                      this->curElemVolVars(),
                                                      this->elemFluxVarsCache(), scvf);
204
205
        };

206
207
        // get the elements in which we need to evaluate the fluxes
        // and calculate these in the undeflected state
208
        unsigned int j = 1;
209
        for (const auto& dataJ : connectivityMap[globalI])
210
        {
211
            neighborElements[j-1] = gridGeometry.element(dataJ.globalJ);
212
            for (const auto scvfIdx : dataJ.scvfsJ)
213
                origResiduals[j] += evalNeighborFlux(neighborElements[j-1], fvGeometry.scvf(scvfIdx));
214
215

            ++j;
216
217
218
219
        }

        // reference to the element's scv (needed later) and corresponding vol vars
        const auto& scv = fvGeometry.scv(globalI);
220
        auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
221
222
223

        // save a copy of the original privars and vol vars in order
        // to restore the original solution after deflection
224
        const auto& curSol = this->curSol();
225
226
227
228
        const auto origPriVars = curSol[globalI];
        const auto origVolVars = curVolVars;

        // element solution container to be deflected
229
        auto elemSol = elementSolution(element, curSol, gridGeometry);
230
231

        // derivatives in the neighbors with repect to the current elements
232
233
234
        // in index 0 we save the derivative of the element residual with respect to it's own dofs
        Residuals partialDerivs(numNeighbors + 1);

235
        for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
236
        {
237
            partialDerivs = 0.0;
238

239
            auto evalResiduals = [&](Scalar priVar)
240
            {
241
242
                Residuals partialDerivsTmp(numNeighbors + 1);
                partialDerivsTmp = 0.0;
243
                // update the volume variables and the flux var cache
244
245
                elemSol[0][pvIdx] = priVar;
                curVolVars.update(elemSol, this->problem(), element, scv);
246
                elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
247
                if (enableGridFluxVarsCache)
248
                    gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
249

250
251
                // calculate the residual with the deflected primary variables
                partialDerivsTmp[0] = this->evalLocalResidual()[0];
252
253
254

                // calculate the fluxes in the neighbors with the deflected primary variables
                for (std::size_t k = 0; k < numNeighbors; ++k)
255
                    for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
256
                        partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k], fvGeometry.scvf(scvfIdx));
257

258
259
                return partialDerivsTmp;
            };
260

261
            // derive the residuals numerically
262
263
            static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
            static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
264
265
            NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals,
                                                      eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
266

267
268
269
270
271
272
273
274
275
            // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the
            // current primary variable and a 0 elsewhere. As we always solve for a delta of the
            // solution with repect to the initial one, this results in a delta of 0 for ghosts.
            if (this->elementIsGhost())
            {
                partialDerivs[0] = 0.0;
                partialDerivs[0][pvIdx] = 1.0;
            }

276
            // add the current partial derivatives to the global jacobian matrix
277
            // no special treatment is needed if globalJ is a ghost because then derivatives have been assembled to 0 above
278
279
280
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
            {
                // the diagonal entries
281
                A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
282
283

                // off-diagonal entries
284
                j = 1;
285
                for (const auto& dataJ : connectivityMap[globalI])
286
                    A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
287
            }
288
289
290
291
292
293

            // restore the original state of the scv's volume variables
            curVolVars = origVolVars;

            // restore the current element solution
            elemSol[0][pvIdx] = origPriVars[pvIdx];
294
295
        }

296
        // restore original state of the flux vars cache in case of global caching.
297
298
299
        // This has to be done in order to guarantee that everything is in an undeflected
        // state before the assembly of another element is called. In the case of local caching
        // this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
300
301
        // We only have to do this for the last primary variable, for all others the flux var cache
        // is updated with the correct element volume variables before residual evaluations
302
303
304
        if (enableGridFluxVarsCache)
            gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);

305
        // return the original residual
306
        return origResiduals[0];
307
308
309
    }
};

310
311
312
313

/*!
 * \ingroup Assembly
 * \ingroup CCDiscretization
314
 * \brief Cell-centered scheme local assembler using numeric differentiation and explicit time discretization
315
 */
316
317
template<class TypeTag, class Assembler>
class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
318
: public CCLocalAssemblerBase<TypeTag, Assembler,
319
            CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>, false>
320
{
321
    using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>;
322
    using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
323
324
325
326
327
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using Element = typename GetPropType<TypeTag, Properties::GridView>::template Codim<0>::Entity;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
328

329
    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
330
331

public:
332
    using ParentType::ParentType;
333
334
335
336
337
338
339

    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix.
     *
     * \return The element residual at the current solution.
     */
340
    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
341
    {
342
        if (this->assembler().isStationaryProblem())
343
344
            DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");

345
        // assemble the undeflected residual
346
        const auto residual = this->evalLocalResidual()[0];
347
        const auto storageResidual = this->evalLocalStorageResidual();
348
349
350
351
352
353
354

        //////////////////////////////////////////////////////////////////////////////////////////////////
        // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
        // neighboring elements all derivatives are zero. For the assembled element only the storage    //
        // derivatives are non-zero.                                                                    //
        //////////////////////////////////////////////////////////////////////////////////////////////////

355
356
357
        // get some aliases for convenience
        const auto& element = this->element();
        const auto& fvGeometry = this->fvGeometry();
358
        const auto& gridGeometry = this->assembler().gridGeometry();
359
360
        auto&& curElemVolVars = this->curElemVolVars();

361
        // reference to the element's scv (needed later) and corresponding vol vars
362
        const auto globalI = gridGeometry.elementMapper().index(element);
363
        const auto& scv = fvGeometry.scv(globalI);
364
        auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
365
366
367

        // save a copy of the original privars and vol vars in order
        // to restore the original solution after deflection
368
        const auto& curSol = this->curSol();
369
370
371
372
        const auto origPriVars = curSol[globalI];
        const auto origVolVars = curVolVars;

        // element solution container to be deflected
373
        auto elemSol = elementSolution(element, curSol, gridGeometry);
374

375
        NumEqVector partialDeriv;
376
377

        // derivatives in the neighbors with repect to the current elements
378
        for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
379
380
        {
            // reset derivatives of element dof with respect to itself
381
            partialDeriv = 0.0;
382

383
            auto evalStorage = [&](Scalar priVar)
384
            {
385
386
387
388
                // update the volume variables and calculate
                // the residual with the deflected primary variables
                elemSol[0][pvIdx] = priVar;
                curVolVars.update(elemSol, this->problem(), element, scv);
389
                return this->evalLocalStorageResidual();
390
            };
391

392
393
            // for non-ghosts compute the derivative numerically
            if (!this->elementIsGhost())
394
            {
395
396
                static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
                static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
397
                NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, storageResidual,
398
                                                          eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
399
            }
400

401
402
403
404
            // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
            // as we always solve for a delta of the solution with repect to the initial solution this
            // results in a delta of zero for ghosts
            else partialDeriv[pvIdx] = 1.0;
405
406
407
408
409

            // add the current partial derivatives to the global jacobian matrix
            // only diagonal entries for explicit jacobians
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
                A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
410
411
412
413
414
415
416
417
418
419
420
421
422

            // restore the original state of the scv's volume variables
            curVolVars = origVolVars;

            // restore the current element solution
            elemSol[0][pvIdx] = origPriVars[pvIdx];
        }

        // return the original residual
        return residual;
    }
};

423
424
425
426
427
/*!
 * \ingroup Assembly
 * \ingroup CCDiscretization
 * \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and implicit time discretization
 */
428
429
template<class TypeTag, class Assembler>
class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
430
: public CCLocalAssemblerBase<TypeTag, Assembler,
431
            CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>, true>
432
{
433
    using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>;
434
    using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
435
436
437
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
438

439
440
    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };

441
public:
442
    using ParentType::ParentType;
443
444
445
446
447
448
449

    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix.
     *
     * \return The element residual at the current solution.
     */
450
    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
451
    {
452
453
454
455
456
457
458
459
460
461
462
        // treat ghost separately, we always want zero update for ghosts
        if (this->elementIsGhost())
        {
            const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
            for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
                A[globalI][globalI][pvIdx][pvIdx] = 1.0;

            // return zero residual
            return NumEqVector(0.0);
        }

463
        // assemble the undeflected residual
464
        const auto residual = this->evalLocalResidual()[0];
465

466
467
468
469
        // get some aliases for convenience
        const auto& problem = this->problem();
        const auto& element = this->element();
        const auto& fvGeometry = this->fvGeometry();
470
471
        const auto& curElemVolVars = this->curElemVolVars();
        const auto& elemFluxVarsCache = this->elemFluxVarsCache();
472
473

        // get reference to the element's current vol vars
474
        const auto globalI = this->assembler().gridGeometry().elementMapper().index(element);
475
476
        const auto& scv = fvGeometry.scv(globalI);
        const auto& volVars = curElemVolVars[scv];
477
478

        // if the problem is instationary, add derivative of storage term
479
480
        if (!this->assembler().isStationaryProblem())
            this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
481
482

        // add source term derivatives
483
        this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
484
485
486
487

        // add flux derivatives for each scvf
        for (const auto& scvf : scvfs(fvGeometry))
        {
488
            // inner faces
489
            if (!scvf.boundary())
490
491
492
                this->localResidual().addFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);

            // boundary faces
493
494
495
496
497
498
            else
            {
                const auto& bcTypes = problem.boundaryTypes(element, scvf);

                // add Dirichlet boundary flux derivatives
                if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
499
500
                    this->localResidual().addCCDirichletFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);

501
502
                // add Robin ("solution dependent Neumann") boundary flux derivatives
                else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
503
504
                    this->localResidual().addRobinFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);

505
506
507
508
509
510
511
512
513
514
                else
                    DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
            }
        }

        // return element residual
        return residual;
    }
};

515
516
517
518
519
/*!
 * \ingroup Assembly
 * \ingroup CCDiscretization
 * \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and explicit time discretization
 */
520
521
template<class TypeTag, class Assembler>
class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/false>
522
: public CCLocalAssemblerBase<TypeTag, Assembler,
523
            CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>, false>
524
{
525
    using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>;
526
    using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
527
528
529
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
530

531
532
    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };

533
public:
534
    using ParentType::ParentType;
535
536
537
538
539
540
541

    /*!
     * \brief Computes the derivatives with respect to the given element and adds them
     *        to the global matrix.
     *
     * \return The element residual at the current solution.
     */
542
    NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
543
    {
544
545
546
547
548
549
550
551
552
553
554
        // treat ghost separately, we always want zero update for ghosts
        if (this->elementIsGhost())
        {
            const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
            for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
                A[globalI][globalI][pvIdx][pvIdx] = 1.0;

            // return zero residual
            return NumEqVector(0.0);
        }

555
        // assemble the undeflected residual
556
        const auto residual = this->evalLocalResidual()[0];
557
558

        // get reference to the element's current vol vars
559
        const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
560
561
        const auto& scv = this->fvGeometry().scv(globalI);
        const auto& volVars = this->curElemVolVars()[scv];
562

563
        // add hand-code derivative of storage term
564
        this->localResidual().addStorageDerivatives(A[globalI][globalI], this->problem(), this->element(), this->fvGeometry(), volVars, scv);
565
566
567
568
569
570

        // return the original residual
        return residual;
    }
};

571
572
573
} // end namespace Dumux

#endif