3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
subdomaincclocalassembler.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 *****************************************************************************/
26#ifndef DUMUX_MULTIDOMAIN_CC_LOCAL_ASSEMBLER_HH
27#define DUMUX_MULTIDOMAIN_CC_LOCAL_ASSEMBLER_HH
28
29#include <dune/common/indices.hh>
30#include <dune/common/reservedvector.hh>
31#include <dune/common/hybridutilities.hh>
32#include <dune/grid/common/gridenums.hh> // for GhostEntity
33#include <dune/istl/matrixindexset.hh>
34
43
44namespace Dumux {
45
57template<std::size_t id, class TypeTag, class Assembler, class Implementation, bool implicit>
58class SubDomainCCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler,Implementation, implicit>
59{
61
63 using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
65 using SolutionVector = typename Assembler::SolutionVector;
68
70 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
71 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
72 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
73 using Scalar = typename GridVariables::Scalar;
74
75 using GridGeometry = typename GridVariables::GridGeometry;
76 using FVElementGeometry = typename GridGeometry::LocalView;
77 using SubControlVolume = typename GridGeometry::SubControlVolume;
78 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
79 using GridView = typename GridGeometry::GridView;
80 using Element = typename GridView::template Codim<0>::Entity;
81
82 using CouplingManager = typename Assembler::CouplingManager;
83 using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
84
85 static_assert(!Problem::enableInternalDirichletConstraints(), "Internal Dirichlet constraints are currently not implemented for cc-methods!");
86
87public:
89 static constexpr auto domainId = typename Dune::index_constant<id>();
93 using ParentType::ParentType;
94
96 explicit SubDomainCCLocalAssemblerBase(const Assembler& assembler,
97 const Element& element,
98 const SolutionVector& curSol,
99 CouplingManager& couplingManager)
101 element,
102 curSol,
103 localView(assembler.gridGeometry(domainId)),
104 localView(assembler.gridVariables(domainId).curGridVolVars()),
105 localView(assembler.gridVariables(domainId).prevGridVolVars()),
106 localView(assembler.gridVariables(domainId).gridFluxVarsCache()),
108 (element.partitionType() == Dune::GhostEntity))
109 , couplingManager_(couplingManager)
110 {}
111
116 template<class JacobianMatrixRow, class GridVariablesTuple>
117 void assembleJacobianAndResidual(JacobianMatrixRow& jacRow, SubSolutionVector& res, GridVariablesTuple& gridVariables)
118 {
119 this->asImp_().bindLocalViews();
120
121 // for the diagonal jacobian block
122 const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element());
123 res[globalI] = this->asImp_().assembleJacobianAndResidualImplInverse(jacRow[domainId], *std::get<domainId>(gridVariables));
124
125 // for the coupling blocks
126 using namespace Dune::Hybrid;
127 forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
128 {
129 if (i != id)
130 this->assembleJacobianCoupling(i, jacRow, res[globalI], gridVariables);
131 });
132 }
133
138 template<std::size_t otherId, class JacRow, class GridVariables,
139 typename std::enable_if_t<(otherId == id), int> = 0>
140 void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
141 const LocalResidualValues& res, GridVariables& gridVariables)
142 {}
143
147 template<std::size_t otherId, class JacRow, class GridVariables,
148 typename std::enable_if_t<(otherId != id), int> = 0>
149 void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
150 const LocalResidualValues& res, GridVariables& gridVariables)
151 {
152 this->asImp_().assembleJacobianCoupling(domainJ, jacRow[domainJ], res, *std::get<domainId>(gridVariables));
153 }
154
158 void assembleResidual(SubSolutionVector& res)
159 {
160 this->asImp_().bindLocalViews();
161 const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element());
162 res[globalI] = this->evalLocalResidual()[0]; // forward to the internal implementation
163 }
164
168 ElementResidualVector evalLocalSourceResidual(const Element& element, const ElementVolumeVariables& elemVolVars) const
169 {
170 // initialize the residual vector for all scvs in this element
171 ElementResidualVector residual(this->fvGeometry().numScv());
172
173 // evaluate the volume terms (storage + source terms)
174 // forward to the local residual specialized for the discretization methods
175 for (auto&& scv : scvs(this->fvGeometry()))
176 {
177 const auto& curVolVars = elemVolVars[scv];
178 auto source = this->localResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv);
179 source *= -scv.volume()*curVolVars.extrusionFactor();
180 residual[scv.indexInElement()] = std::move(source);
181 }
182
183 return residual;
184 }
185
189 ElementResidualVector evalLocalSourceResidual(const Element& neighbor) const
190 { return this->evalLocalSourceResidual(neighbor, implicit ? this->curElemVolVars() : this->prevElemVolVars()); }
191
195 LocalResidualValues evalLocalStorageResidual() const
196 {
197 return this->localResidual().evalStorage(this->element(), this->fvGeometry(), this->prevElemVolVars(), this->curElemVolVars())[0];
198 }
199
203 LocalResidualValues evalFluxResidual(const Element& neighbor,
204 const SubControlVolumeFace& scvf) const
205 {
206 const auto& elemVolVars = implicit ? this->curElemVolVars() : this->prevElemVolVars();
207 return this->localResidual().evalFlux(problem(), neighbor, this->fvGeometry(), elemVolVars, this->elemFluxVarsCache(), scvf);
208 }
209
214 {
215 // get some references for convenience
216 const auto& element = this->element();
217 const auto& curSol = this->curSol()[domainId];
218 auto&& fvGeometry = this->fvGeometry();
219 auto&& curElemVolVars = this->curElemVolVars();
220 auto&& elemFluxVarsCache = this->elemFluxVarsCache();
221
222 // bind the caches
223 couplingManager_.bindCouplingContext(domainId, element, this->assembler());
224 fvGeometry.bind(element);
225
226 if (implicit)
227 {
230 if (!this->assembler().isStationaryProblem())
231 this->prevElemVolVars().bindElement(element, fvGeometry, this->assembler().prevSol()[domainId]);
232 }
233 else
234 {
235 auto& prevElemVolVars = this->prevElemVolVars();
236 const auto& prevSol = this->assembler().prevSol()[domainId];
237
239 prevElemVolVars.bind(element, fvGeometry, prevSol);
241 }
242 }
243
245 const Problem& problem() const
246 { return this->assembler().problem(domainId); }
247
249 CouplingManager& couplingManager()
250 { return couplingManager_; }
251
252private:
253 CouplingManager& couplingManager_;
254};
255
267template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
269
276template<std::size_t id, class TypeTag, class Assembler>
277class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
278: public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler,
280{
281 using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
282 using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/true>;
283
285 using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
286
288 using FVElementGeometry = typename GridGeometry::LocalView;
289 using GridView = typename GridGeometry::GridView;
290 using Element = typename GridView::template Codim<0>::Entity;
291
293 enum { dim = GridView::dimension };
294
295 static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
296 static constexpr bool enableGridVolVarsCache = getPropValue<TypeTag, Properties::EnableGridVolumeVariablesCache>();
297 static constexpr int maxElementStencilSize = GridGeometry::maxElementStencilSize;
298 static constexpr auto domainI = Dune::index_constant<id>();
299
300public:
302
309 template<class JacobianMatrixDiagBlock, class GridVariables>
310 LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
311 {
313 // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
314 // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
315 // actual element. In the actual element we evaluate the derivative of the entire residual. //
317
318 // get some aliases for convenience
319 const auto& element = this->element();
320 const auto& fvGeometry = this->fvGeometry();
321 const auto& gridGeometry = fvGeometry.gridGeometry();
322 auto&& curElemVolVars = this->curElemVolVars();
323 auto&& elemFluxVarsCache = this->elemFluxVarsCache();
324
325 // get stencil informations
326 const auto globalI = gridGeometry.elementMapper().index(element);
327 const auto& connectivityMap = gridGeometry.connectivityMap();
328 const auto numNeighbors = connectivityMap[globalI].size();
329
330 // container to store the neighboring elements
331 Dune::ReservedVector<Element, maxElementStencilSize> neighborElements;
332 neighborElements.resize(numNeighbors);
333
334 // assemble the undeflected residual
336 Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
337 origResiduals[0] = this->evalLocalResidual()[0];
338
339 // get the elements in which we need to evaluate the fluxes
340 // and calculate these in the undeflected state
341 unsigned int j = 1;
342 for (const auto& dataJ : connectivityMap[globalI])
343 {
344 neighborElements[j-1] = gridGeometry.element(dataJ.globalJ);
345 for (const auto scvfIdx : dataJ.scvfsJ)
346 origResiduals[j] += this->evalFluxResidual(neighborElements[j-1], fvGeometry.scvf(scvfIdx));
347
348 ++j;
349 }
350
351 // reference to the element's scv (needed later) and corresponding vol vars
352 const auto& scv = fvGeometry.scv(globalI);
353 auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
354
355 // save a copy of the original privars and vol vars in order
356 // to restore the original solution after deflection
357 const auto& curSol = this->curSol()[domainI];
358 const auto origPriVars = curSol[globalI];
359 const auto origVolVars = curVolVars;
360
361 // element solution container to be deflected
362 auto elemSol = elementSolution<FVElementGeometry>(origPriVars);
363
364 // derivatives in the neighbors with repect to the current elements
365 // in index 0 we save the derivative of the element residual with respect to it's own dofs
366 Residuals partialDerivs(numNeighbors + 1);
367
368 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
369 {
370 // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
371 // as we always solve for a delta of the solution with repect to the initial solution this
372 // results in a delta of zero for ghosts, we still need to do the neighbor derivatives though
373 // so we are not done yet here.
374 partialDerivs = 0.0;
375 if (this->elementIsGhost()) partialDerivs[0][pvIdx] = 1.0;
376
377 auto evalResiduals = [&](Scalar priVar)
378 {
379 Residuals partialDerivsTmp(numNeighbors + 1);
380 partialDerivsTmp = 0.0;
381 // update the volume variables and the flux var cache
382 elemSol[0][pvIdx] = priVar;
383 this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx);
384 curVolVars.update(elemSol, this->problem(), element, scv);
385 elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
386 if (enableGridFluxVarsCache)
387 gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
388
389 // calculate the residual with the deflected primary variables
390 if (!this->elementIsGhost()) partialDerivsTmp[0] = this->evalLocalResidual()[0];
391
392 // calculate the fluxes in the neighbors with the deflected primary variables
393 for (std::size_t k = 0; k < connectivityMap[globalI].size(); ++k)
394 for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
395 partialDerivsTmp[k+1] += this->evalFluxResidual(neighborElements[k], fvGeometry.scvf(scvfIdx));
396
397 return partialDerivsTmp;
398 };
399
400 // derive the residuals numerically
401 static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
402 static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
403 NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals,
404 eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
405
406 // add the current partial derivatives to the global jacobian matrix
407 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
408 {
409 // the diagonal entries
410 A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
411
412 // off-diagonal entries
413 j = 1;
414 for (const auto& dataJ : connectivityMap[globalI])
415 A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
416 }
417
418 // restore the original state of the scv's volume variables
419 curVolVars = origVolVars;
420
421 // restore the current element solution
422 elemSol[0][pvIdx] = origPriVars[pvIdx];
423
424 // restore the undeflected state of the coupling context
425 this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx);
426 }
427
428 // restore original state of the flux vars cache in case of global caching.
429 // This has to be done in order to guarantee that everything is in an undeflected
430 // state before the assembly of another element is called. In the case of local caching
431 // this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
432 // We only have to do this for the last primary variable, for all others the flux var cache
433 // is updated with the correct element volume variables before residual evaluations
434 if (enableGridFluxVarsCache)
435 gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
436
437 // compute potential additional derivatives causes by the coupling
438 typename ParentType::LocalResidual::ElementResidualVector origElementResidual({origResiduals[0]});
439 this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origElementResidual, A, gridVariables);
440
441 // return the original residual
442 return origResiduals[0];
443 }
444
449 template<std::size_t otherId, class JacobianBlock, class GridVariables>
450 void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
451 const LocalResidualValues& res, GridVariables& gridVariables)
452 {
454 // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
456
457 // get some aliases for convenience
458 const auto& element = this->element();
459 const auto& fvGeometry = this->fvGeometry();
460 const auto& gridGeometry = fvGeometry.gridGeometry();
461 auto&& curElemVolVars = this->curElemVolVars();
462 auto&& elemFluxVarsCache = this->elemFluxVarsCache();
463
464 // get stencil informations
465 const auto globalI = gridGeometry.elementMapper().index(element);
466 const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
467 const auto& curSolJ = this->curSol()[domainJ];
468
469 // make sure the flux variables cache does not contain any artifacts from prior differentiation
470 elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
471
472 // convenience lambda for call to update self
473 auto updateCoupledVariables = [&] ()
474 {
475 // Update ourself after the context has been modified. Depending on the
476 // type of caching, other objects might have to be updated. All ifs can be optimized away.
477 if (enableGridFluxVarsCache)
478 {
479 if (enableGridVolVarsCache)
480 {
481 this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
482 this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
483 }
484 else
485 {
486 this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache());
487 this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
488 }
489 }
490 else
491 {
492 if (enableGridVolVarsCache)
493 this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
494 else
495 this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
496 }
497 };
498
499 for (const auto& globalJ : stencil)
500 {
501 // undeflected privars and privars to be deflected
502 const auto origPriVarsJ = curSolJ[globalJ];
503 auto priVarsJ = origPriVarsJ;
504
505 const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0];
506
507 for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx)
508 {
509 auto evalCouplingResidual = [&](Scalar priVar)
510 {
511 // update the volume variables and the flux var cache
512 priVarsJ[pvIdx] = priVar;
513 this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
514 updateCoupledVariables();
515 return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0];
516 };
517
518 // derive the residuals numerically
519 LocalResidualValues partialDeriv(0.0);
520 const auto& paramGroup = this->assembler().problem(domainJ).paramGroup();
521 static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod");
522 static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup);
523 NumericDifferentiation::partialDerivative(evalCouplingResidual, priVarsJ[pvIdx], partialDeriv, origResidual,
524 epsCoupl(priVarsJ[pvIdx], pvIdx), numDiffMethod);
525
526 // add the current partial derivatives to the global jacobian matrix
527 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
528 A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
529
530 // restore the current element solution
531 priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
532
533 // restore the undeflected state of the coupling context
534 this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
535 }
536
537 // restore original state of the flux vars cache and/or vol vars in case of global caching.
538 // This has to be done in order to guarantee that the undeflected residual computation done
539 // above is correct when jumping to the next coupled dof of domainJ.
540 updateCoupledVariables();
541 }
542 }
543};
544
551template<std::size_t id, class TypeTag, class Assembler>
552class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
553: public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler,
554 SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, false>, false>
555{
556 using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>;
557 using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/false>;
558
560 using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
561
562 static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
563 static constexpr auto domainI = Dune::index_constant<id>();
564
565public:
567
578 template<class JacobianMatrixDiagBlock, class GridVariables>
579 LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
580 {
581 if (this->assembler().isStationaryProblem())
582 DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
583
584 // assemble the undeflected residual
585 const auto residual = this->evalLocalResidual()[0];
586 const auto storageResidual = this->evalLocalStorageResidual();
587
589 // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
590 // neighboring elements all derivatives are zero. For the assembled element only the storage //
591 // derivatives are non-zero. //
593
594 // get some aliases for convenience
595 const auto& element = this->element();
596 const auto& fvGeometry = this->fvGeometry();
597 const auto& gridGeometry = fvGeometry.gridGeometry();
598 auto&& curElemVolVars = this->curElemVolVars();
599
600 // reference to the element's scv (needed later) and corresponding vol vars
601 const auto globalI = gridGeometry.elementMapper().index(element);
602 const auto& scv = fvGeometry.scv(globalI);
603 auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
604
605 // save a copy of the original privars and vol vars in order
606 // to restore the original solution after deflection
607 const auto& curSol = this->curSol()[domainI];
608 const auto origPriVars = curSol[globalI];
609 const auto origVolVars = curVolVars;
610
611 // element solution container to be deflected
612 auto elemSol = elementSolution(element, curSol, gridGeometry);
613
614 // derivatives in the neighbors with repect to the current elements
615 LocalResidualValues partialDeriv;
616 for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
617 {
618 // reset derivatives of element dof with respect to itself
619 partialDeriv = 0.0;
620
621 auto evalStorage = [&](Scalar priVar)
622 {
623 // update the volume variables and calculate
624 // the residual with the deflected primary variables
625 elemSol[0][pvIdx] = priVar;
626 curVolVars.update(elemSol, this->problem(), element, scv);
627 return this->evalLocalStorageResidual();
628 };
629
630 // for non-ghosts compute the derivative numerically
631 if (!this->elementIsGhost())
632 {
633 static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
634 static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
635 NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, storageResidual,
636 eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
637 }
638
639 // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
640 // as we always solve for a delta of the solution with repect to the initial solution this
641 // results in a delta of zero for ghosts
642 else partialDeriv[pvIdx] = 1.0;
643
644 // add the current partial derivatives to the global jacobian matrix
645 // only diagonal entries for explicit jacobians
646 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
647 A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
648
649 // restore the original state of the scv's volume variables
650 curVolVars = origVolVars;
651
652 // restore the current element solution
653 elemSol[0][pvIdx] = origPriVars[pvIdx];
654 }
655
656 // return the original residual
657 return residual;
658 }
659
666 template<std::size_t otherId, class JacobianBlock, class GridVariables>
667 void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
668 const LocalResidualValues& res, GridVariables& gridVariables)
669 {}
670};
671
677template<std::size_t id, class TypeTag, class Assembler>
678class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
679: public SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler,
680 SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, true>, true>
681{
682 using ThisType = SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>;
683 using ParentType = SubDomainCCLocalAssemblerBase<id, TypeTag, Assembler, ThisType, /*implicit=*/true>;
685 using LocalResidualValues = GetPropType<TypeTag, Properties::NumEqVector>;
687 using Element = typename GridView::template Codim<0>::Entity;
688
690 enum { dim = GridView::dimension };
691
692 static constexpr auto domainI = Dune::index_constant<id>();
693
694public:
696
703 template<class JacobianMatrixDiagBlock, class GridVariables>
704 LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
705 {
706 // get some aliases for convenience
707 const auto& problem = this->problem();
708 const auto& element = this->element();
709 const auto& fvGeometry = this->fvGeometry();
710 const auto& curElemVolVars = this->curElemVolVars();
711 const auto& elemFluxVarsCache = this->elemFluxVarsCache();
712
713 // get reference to the element's current vol vars
714 const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(element);
715 const auto& scv = fvGeometry.scv(globalI);
716 const auto& volVars = curElemVolVars[scv];
717
718 // if the problem is instationary, add derivative of storage term
719 if (!this->assembler().isStationaryProblem())
720 this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
721
722 // add source term derivatives
723 this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
724
725 // add flux derivatives for each scvf
726 for (const auto& scvf : scvfs(fvGeometry))
727 {
728 // inner faces
729 if (!scvf.boundary())
730 this->localResidual().addFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
731
732 // boundary faces
733 else
734 {
735 const auto& bcTypes = problem.boundaryTypes(element, scvf);
736
737 // add Dirichlet boundary flux derivatives
738 if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
739 this->localResidual().addCCDirichletFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
740
741 // add Robin ("solution dependent Neumann") boundary flux derivatives
742 else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
743 this->localResidual().addRobinFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
744
745 else
746 DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
747 }
748 }
749
750 // return element residual
751 return this->evalLocalResidual()[0];
752 }
753
760 template<std::size_t otherId, class JacobianBlock, class GridVariables>
761 void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
762 const LocalResidualValues& res, GridVariables& gridVariables)
763 {
765 // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
767
768 // get some aliases for convenience
769 const auto& element = this->element();
770 const auto& fvGeometry = this->fvGeometry();
771 const auto& gridGeometry = fvGeometry.gridGeometry();
772 auto&& curElemVolVars = this->curElemVolVars();
773 // auto&& elemFluxVarsCache = this->elemFluxVarsCache();
774
775 // get stencil informations
776 const auto globalI = gridGeometry.elementMapper().index(element);
777 const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
778
779 for (const auto globalJ : stencil)
780 {
781 const auto& elementJ = this->assembler().gridGeometry(domainJ).element(globalJ);
782 this->couplingManager().addCouplingDerivatives(A[globalI][globalJ], domainI, element, fvGeometry, curElemVolVars, domainJ, elementJ);
783 }
784 }
785};
786
787} // end namespace Dumux
788
789#endif
An enum class to define various differentiation methods available in order to compute the derivatives...
An adapter class for local assemblers using numeric differentiation.
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.
Element solution classes and factory functions.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
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
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
Definition: common/pdesolver.hh:35
A base class for all local assemblers.
Definition: fvlocalassemblerbase.hh:48
ElementVolumeVariables & curElemVolVars()
The current element volume variables.
Definition: fvlocalassemblerbase.hh:263
typename LocalResidual::ElementResidualVector ElementResidualVector
Definition: fvlocalassemblerbase.hh:68
Implementation & asImp_()
Definition: fvlocalassemblerbase.hh:307
ElementVolumeVariables & prevElemVolVars()
The element volume variables of the provious time step.
Definition: fvlocalassemblerbase.hh:267
ElementResidualVector evalLocalResidual() const
Convenience function to evaluate the complete local residual for the current element....
Definition: fvlocalassemblerbase.hh:120
FVElementGeometry & fvGeometry()
The global finite volume geometry.
Definition: fvlocalassemblerbase.hh:259
const Assembler & assembler() const
The assembler.
Definition: fvlocalassemblerbase.hh:243
ElementFluxVariablesCache & elemFluxVarsCache()
The element flux variables cache.
Definition: fvlocalassemblerbase.hh:271
LocalResidual & localResidual()
The local residual for the current element.
Definition: fvlocalassemblerbase.hh:275
const Element & element() const
The current element.
Definition: fvlocalassemblerbase.hh:247
const SolutionVector & curSol() const
The current solution.
Definition: fvlocalassemblerbase.hh:255
A helper class for local assemblers using numeric differentiation to determine the epsilon.
Definition: numericepsilon.hh:41
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
Definition: multidomain/couplingmanager.hh:46
A base class for all multidomain local assemblers.
Definition: subdomaincclocalassembler.hh:59
void bindLocalViews()
Prepares all local views necessary for local assembly.
Definition: subdomaincclocalassembler.hh:213
SubDomainCCLocalAssemblerBase(const Assembler &assembler, const Element &element, const SolutionVector &curSol, CouplingManager &couplingManager)
the constructor
Definition: subdomaincclocalassembler.hh:96
ElementResidualVector evalLocalSourceResidual(const Element &neighbor) const
Evaluates the local source term depending on time discretization scheme.
Definition: subdomaincclocalassembler.hh:189
LocalResidualValues evalFluxResidual(const Element &neighbor, const SubControlVolumeFace &scvf) const
Evaluates the fluxes depending on the chose time discretization scheme.
Definition: subdomaincclocalassembler.hh:203
void assembleJacobianAndResidual(JacobianMatrixRow &jacRow, SubSolutionVector &res, GridVariablesTuple &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix....
Definition: subdomaincclocalassembler.hh:117
void assembleJacobianCoupling(Dune::index_constant< otherId > domainJ, JacRow &jacRow, const LocalResidualValues &res, GridVariables &gridVariables)
Assemble the entries in a coupling block of the jacobian. There is no coupling block between a domain...
Definition: subdomaincclocalassembler.hh:140
const Problem & problem() const
return reference to the underlying problem
Definition: subdomaincclocalassembler.hh:245
CouplingManager & couplingManager()
return reference to the coupling manager
Definition: subdomaincclocalassembler.hh:249
static constexpr auto domainId
export the domain id of this sub-domain
Definition: subdomaincclocalassembler.hh:89
GetPropType< TypeTag, Properties::LocalResidual > LocalResidual
the local residual type of this domain
Definition: subdomaincclocalassembler.hh:91
ElementResidualVector evalLocalSourceResidual(const Element &element, const ElementVolumeVariables &elemVolVars) const
Evaluates the local source term for an element and given element volume variables.
Definition: subdomaincclocalassembler.hh:168
LocalResidualValues evalLocalStorageResidual() const
Evaluates the storage terms within the element.
Definition: subdomaincclocalassembler.hh:195
void assembleResidual(SubSolutionVector &res)
Assemble the residual only.
Definition: subdomaincclocalassembler.hh:158
The cell-centered scheme multidomain local assembler.
Definition: subdomaincclocalassembler.hh:268
Cell-centered scheme multidomain local assembler using numeric differentiation and implicit time disc...
Definition: subdomaincclocalassembler.hh:280
void assembleJacobianCoupling(Dune::index_constant< otherId > domainJ, JacobianBlock &A, const LocalResidualValues &res, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: subdomaincclocalassembler.hh:450
LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: subdomaincclocalassembler.hh:310
Cell-centered scheme multidomain local assembler using numeric differentiation and explicit time disc...
Definition: subdomaincclocalassembler.hh:555
LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: subdomaincclocalassembler.hh:579
void assembleJacobianCoupling(Dune::index_constant< otherId > domainJ, JacobianBlock &A, const LocalResidualValues &res, GridVariables &gridVariables)
Computes the coupling derivatives with respect to the given element and adds them to the global matri...
Definition: subdomaincclocalassembler.hh:667
Cell-centered scheme local assembler using analytic differentiation and implicit time discretization.
Definition: subdomaincclocalassembler.hh:681
void assembleJacobianCoupling(Dune::index_constant< otherId > domainJ, JacobianBlock &A, const LocalResidualValues &res, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: subdomaincclocalassembler.hh:761
LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
Computes the derivatives with respect to the given element and adds them to the global matrix.
Definition: subdomaincclocalassembler.hh:704
Declares all properties used in Dumux.