3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
cclocalassembler.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * 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 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
26#define DUMUX_CC_LOCAL_ASSEMBLER_HH
27
28#include <dune/common/reservedvector.hh>
29#include <dune/grid/common/gridenums.hh> // for GhostEntity
30#include <dune/istl/matrixindexset.hh>
31
43
44namespace Dumux {
45
55template<class TypeTag, class Assembler, class Implementation, bool implicit>
56class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
57{
63 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
65
66 static_assert(!Assembler::Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");
67
68public:
69
70 using ParentType::ParentType;
71
76 template <class PartialReassembler = DefaultPartialReassembler>
77 void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables,
78 const PartialReassembler* partialReassembler)
79 {
80 this->asImp_().bindLocalViews();
81 const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
82 if (partialReassembler
83 && partialReassembler->elementColor(globalI) == EntityColor::green)
84 {
85 res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
86 }
87 else
88 {
89 res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
90 }
91 }
92
97 void assembleJacobian(JacobianMatrix& jac, GridVariables& gridVariables)
98 {
99 this->asImp_().bindLocalViews();
100 this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
101 }
102
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
122template<class TypeTag, class Assembler, DiffMethod diffMethod = DiffMethod::numeric, bool implicit = true>
124
130template<class TypeTag, class Assembler>
131class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
132: public CCLocalAssemblerBase<TypeTag, Assembler,
133 CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>, true >
134{
139 using Element = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::template Codim<0>::Entity;
141 using FVElementGeometry = typename GridGeometry::LocalView;
144
147
149 static constexpr int maxElementStencilSize = GridGeometry::maxElementStencilSize;
150 static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
151
152public:
153
155
162 NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
163 {
165 // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
166 // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
167 // actual element. In the actual element we evaluate the derivative of the entire residual. //
169
170 // get some aliases for convenience
171 const auto& element = this->element();
172 const auto& fvGeometry = this->fvGeometry();
173 const auto& gridGeometry = this->assembler().gridGeometry();
174 auto&& curElemVolVars = this->curElemVolVars();
175 auto&& elemFluxVarsCache = this->elemFluxVarsCache();
176
177 // get stencil informations
178 const auto globalI = gridGeometry.elementMapper().index(element);
179 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
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 // if the neighbor is a ghost we don't want to add anything to their residual
193 // so we return 0 and omit computing the flux
194 auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
195 {
196 if (neighbor.partitionType() == Dune::GhostEntity)
197 return NumEqVector(0.0);
198 else
199 return this->localResidual().evalFlux(this->problem(),
200 neighbor,
201 this->fvGeometry(),
202 this->curElemVolVars(),
203 this->elemFluxVarsCache(), scvf);
204 };
205
206 // get the elements in which we need to evaluate the fluxes
207 // 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 // reference to the element's scv (needed later) and corresponding vol vars
219 const auto& scv = fvGeometry.scv(globalI);
220 auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
221
222 // save a copy of the original privars and vol vars in order
223 // to restore the original solution after deflection
224 const auto& curSol = this->curSol();
225 const auto origPriVars = curSol[globalI];
226 const auto origVolVars = curVolVars;
227
228 // 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 // in index 0 we save the derivative of the element residual with respect to it's own dofs
233 Residuals partialDerivs(numNeighbors + 1);
234
235 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
236 {
237 partialDerivs = 0.0;
238
239 auto evalResiduals = [&](Scalar priVar)
240 {
241 Residuals partialDerivsTmp(numNeighbors + 1);
242 partialDerivsTmp = 0.0;
243 // update the volume variables and the flux var cache
244 elemSol[0][pvIdx] = priVar;
245 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 // calculate the residual with the deflected primary variables
251 partialDerivsTmp[0] = this->evalLocalResidual()[0];
252
253 // calculate the fluxes in the neighbors with the deflected primary variables
254 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 return partialDerivsTmp;
259 };
260
261 // derive the residuals numerically
262 static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
263 static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
264 NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals,
265 eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
266
267 // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the
268 // current primary variable and a 0 elsewhere. As we always solve for a delta of the
269 // solution with repect to the initial one, this results in a delta of 0 for ghosts.
270 if (this->elementIsGhost())
271 {
272 partialDerivs[0] = 0.0;
273 partialDerivs[0][pvIdx] = 1.0;
274 }
275
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 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
279 {
280 // 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 // restore the original state of the scv's volume variables
290 curVolVars = origVolVars;
291
292 // restore the current element solution
293 elemSol[0][pvIdx] = origPriVars[pvIdx];
294 }
295
296 // restore original state of the flux vars cache in case of global caching.
297 // This has to be done in order to guarantee that everything is in an undeflected
298 // state before the assembly of another element is called. In the case of local caching
299 // this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
300 // We only have to do this for the last primary variable, for all others the flux var cache
301 // is updated with the correct element volume variables before residual evaluations
302 if (enableGridFluxVarsCache)
303 gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
304
305 // return the original residual
306 return origResiduals[0];
307 }
308};
309
310
316template<class TypeTag, class Assembler>
317class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
318: public CCLocalAssemblerBase<TypeTag, Assembler,
319 CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>, false>
320{
325 using Element = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::template Codim<0>::Entity;
328
330
331public:
333
340 NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
341 {
342 if (this->assembler().isStationaryProblem())
343 DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
344
345 // assemble the undeflected residual
346 const auto residual = this->evalLocalResidual()[0];
347 const auto storageResidual = this->evalLocalStorageResidual()[0];
348
350 // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
351 // neighboring elements all derivatives are zero. For the assembled element only the storage //
352 // derivatives are non-zero. //
354
355 // get some aliases for convenience
356 const auto& element = this->element();
357 const auto& fvGeometry = this->fvGeometry();
358 const auto& gridGeometry = this->assembler().gridGeometry();
359 auto&& curElemVolVars = this->curElemVolVars();
360
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 // save a copy of the original privars and vol vars in order
367 // to restore the original solution after deflection
368 const auto& curSol = this->curSol();
369 const auto origPriVars = curSol[globalI];
370 const auto origVolVars = curVolVars;
371
372 // 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 // update the volume variables and calculate
386 // the residual with the deflected primary variables
387 elemSol[0][pvIdx] = priVar;
388 curVolVars.update(elemSol, this->problem(), element, scv);
389 return this->evalLocalStorageResidual()[0];
390 };
391
392 // for non-ghosts compute the derivative numerically
393 if (!this->elementIsGhost())
394 {
395 static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
396 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 // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
402 // as we always solve for a delta of the solution with repect to the initial solution this
403 // results in a delta of zero for ghosts
404 else partialDeriv[pvIdx] = 1.0;
405
406 // add the current partial derivatives to the global jacobian matrix
407 // only diagonal entries for explicit jacobians
408 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
409 A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
410
411 // restore the original state of the scv's volume variables
412 curVolVars = origVolVars;
413
414 // restore the current element solution
415 elemSol[0][pvIdx] = origPriVars[pvIdx];
416 }
417
418 // return the original residual
419 return residual;
420 }
421};
422
428template<class TypeTag, class Assembler>
429class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
430: public CCLocalAssemblerBase<TypeTag, Assembler,
431 CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>, true>
432{
438
440
441public:
443
450 NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
451 {
452 // treat ghost separately, we always want zero update for ghosts
453 if (this->elementIsGhost())
454 {
455 const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
456 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
457 A[globalI][globalI][pvIdx][pvIdx] = 1.0;
458
459 // return zero residual
460 return NumEqVector(0.0);
461 }
462
463 // assemble the undeflected residual
464 const auto residual = this->evalLocalResidual()[0];
465
466 // get some aliases for convenience
467 const auto& problem = this->problem();
468 const auto& element = this->element();
469 const auto& fvGeometry = this->fvGeometry();
470 const auto& curElemVolVars = this->curElemVolVars();
471 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 const auto& scv = fvGeometry.scv(globalI);
476 const auto& volVars = curElemVolVars[scv];
477
478 // if the problem is instationary, add derivative of storage term
479 if (!this->assembler().isStationaryProblem())
480 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 // add flux derivatives for each scvf
486 for (const auto& scvf : scvfs(fvGeometry))
487 {
488 // inner faces
489 if (!scvf.boundary())
490 this->localResidual().addFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
491
492 // boundary faces
493 else
494 {
495 const auto& bcTypes = problem.boundaryTypes(element, scvf);
496
497 // add Dirichlet boundary flux derivatives
498 if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
499 this->localResidual().addCCDirichletFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
500
501 // add Robin ("solution dependent Neumann") boundary flux derivatives
502 else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
503 this->localResidual().addRobinFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
504
505 else
506 DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
507 }
508 }
509
510 // return element residual
511 return residual;
512 }
513};
514
520template<class TypeTag, class Assembler>
521class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/false>
522: public CCLocalAssemblerBase<TypeTag, Assembler,
523 CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>, false>
524{
530
532
533public:
535
542 NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
543 {
544 // treat ghost separately, we always want zero update for ghosts
545 if (this->elementIsGhost())
546 {
547 const auto globalI = this->assembler().gridGeometry().elementMapper().index(this->element());
548 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
549 A[globalI][globalI][pvIdx][pvIdx] = 1.0;
550
551 // return zero residual
552 return NumEqVector(0.0);
553 }
554
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 const auto& scv = this->fvGeometry().scv(globalI);
561 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 // return the original residual
567 return residual;
568 }
569};
570
571} // end namespace Dumux
572
573#endif
An enum class to define various differentiation methods available in order to compute the derivatives...
An enum class to define the colors of elements and vertices required for partial Jacobian reassembly.
An adapter class for local assemblers using numeric differentiation.
Detects which entries in the Jacobian have to be recomputed.
A class for numeric differentiation.
A arithmetic block vector type based on DUNE's reserved vector.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The flux stencil specialized for different discretization schemes.
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:115
DiffMethod
Differentiation methods in order to compute the derivatives of the residual i.e. the entries in the j...
Definition: diffmethod.hh:37
@ green
does not need to be reassembled
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
A base class for all local cell-centered assemblers.
Definition: cclocalassembler.hh:57
void assembleResidual(SolutionVector &res)
Assemble the residual only.
Definition: cclocalassembler.hh:106
void assembleJacobian(JacobianMatrix &jac, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: cclocalassembler.hh:97
void assembleJacobianAndResidual(JacobianMatrix &jac, SolutionVector &res, GridVariables &gridVariables, const PartialReassembler *partialReassembler)
Computes the derivatives with respect to the given element and adds them to the global matrix....
Definition: cclocalassembler.hh:77
An assembler for Jacobian and residual contribution per element (cell-centered methods)
Definition: cclocalassembler.hh:123
Cell-centered scheme local assembler using numeric differentiation and implicit time discretization.
Definition: cclocalassembler.hh:134
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix &A, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: cclocalassembler.hh:162
Cell-centered scheme local assembler using numeric differentiation and explicit time discretization.
Definition: cclocalassembler.hh:320
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix &A, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: cclocalassembler.hh:340
Cell-centered scheme local assembler using analytic (hand-coded) differentiation and implicit time di...
Definition: cclocalassembler.hh:432
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix &A, const GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: cclocalassembler.hh:450
Cell-centered scheme local assembler using analytic (hand-coded) differentiation and explicit time di...
Definition: cclocalassembler.hh:524
NumEqVector assembleJacobianAndResidualImpl(JacobianMatrix &A, const GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: cclocalassembler.hh:542
A base class for all local assemblers.
Definition: fvlocalassemblerbase.hh:48
Implementation & asImp_()
Definition: fvlocalassemblerbase.hh:307
const Assembler & assembler() const
The assembler.
Definition: fvlocalassemblerbase.hh:243
const Element & element() const
The current element.
Definition: fvlocalassemblerbase.hh:247
A helper class for local assemblers using numeric differentiation to determine the epsilon.
Definition: numericepsilon.hh:41
detects which entries in the Jacobian have to be recomputed
Definition: partialreassembler.hh:424
EntityColor elementColor(size_t idx) const
Definition: partialreassembler.hh:496
static void partialDerivative(const Function &function, Scalar x0, FunctionEvalType &derivative, const FunctionEvalType &fx0, const int numericDifferenceMethod=1)
Computes the derivative of a function with repect to a function parameter.
Definition: numericdifferentiation.hh:61
A arithmetic block vector type based on DUNE's reserved vector.
Definition: reservedblockvector.hh:38
The flux stencil specialized for different discretization schemes.
Definition: fluxstencil.hh:45
Declares all properties used in Dumux.
The local element solution class for cell-centered methods.