3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
fvlocalresidual.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 *****************************************************************************/
24#ifndef DUMUX_FV_LOCAL_RESIDUAL_HH
25#define DUMUX_FV_LOCAL_RESIDUAL_HH
26
27#include <dune/common/exceptions.hh>
28#include <dune/istl/bvector.hh>
29
34
35namespace Dumux {
36
43template<class TypeTag>
45{
50 using Element = typename GridView::template Codim<0>::Entity;
51 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
54 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
55 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
58 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
60 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
63
64public:
67
69 FVLocalResidual(const Problem* problem,
70 const TimeLoop* timeLoop = nullptr)
71 : problem_(problem)
72 , timeLoop_(timeLoop)
73 {}
74
80 // \{
81
96 const Element &element,
97 const GridGeometry& gridGeometry,
98 const GridVariables& gridVariables,
99 const SolutionVector& sol) const
100 {
101 // make sure FVElementGeometry and volume variables are bound to the element
102 auto fvGeometry = localView(gridGeometry);
103 fvGeometry.bind(element);
104
105 auto elemVolVars = localView(gridVariables.curGridVolVars());
106 elemVolVars.bind(element, fvGeometry, sol);
107
108 ElementResidualVector storage(fvGeometry.numScv());
109
110 // calculate the amount of conservation each quantity inside
111 // all sub control volumes
112 for (auto&& scv : scvs(fvGeometry))
113 {
114 auto localScvIdx = scv.localDofIndex();
115 const auto& volVars = elemVolVars[scv];
116 storage[localScvIdx] = asImp().computeStorage(problem, scv, volVars);
117 storage[localScvIdx] *= scv.volume() * volVars.extrusionFactor();
118 }
119
120 return storage;
121 }
122
123 // \}
124
129 // \{
130
143 ElementResidualVector evalStorage(const Element& element,
144 const FVElementGeometry& fvGeometry,
145 const ElementVolumeVariables& prevElemVolVars,
146 const ElementVolumeVariables& curElemVolVars) const
147 {
148 assert(timeLoop_ && "no time loop set for storage term evaluation");
149
150 // initialize the residual vector for all scvs in this element
151 ElementResidualVector residual(fvGeometry.numScv());
152
153 // evaluate the volume terms (storage + source terms)
154 // forward to the local residual specialized for the discretization methods
155 for (auto&& scv : scvs(fvGeometry))
156 asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
157
158 return residual;
159 }
160
162 const FVElementGeometry& fvGeometry,
163 const ElementVolumeVariables& elemVolVars,
164 const ElementFluxVariablesCache& elemFluxVarsCache,
165 const ElementBoundaryTypes &bcTypes) const
166 {
167 // initialize the residual vector for all scvs in this element
168 ElementResidualVector residual(fvGeometry.numScv());
169
170 // evaluate the volume terms (storage + source terms)
171 // forward to the local residual specialized for the discretization methods
172 for (auto&& scv : scvs(fvGeometry))
173 asImp().evalSource(residual, this->problem(), element, fvGeometry, elemVolVars, scv);
174
175 // forward to the local residual specialized for the discretization methods
176 for (auto&& scvf : scvfs(fvGeometry))
177 asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVolVars, bcTypes, elemFluxVarsCache, scvf);
178
179 return residual;
180 }
181
182 // \}
183
184
189 // \{
190
200 NumEqVector computeStorage(const Problem& problem,
201 const SubControlVolume& scv,
202 const VolumeVariables& volVars) const
203 {
204 DUNE_THROW(Dune::NotImplemented, "This model does not implement a storage method!");
205 }
206
220 NumEqVector computeSource(const Problem& problem,
221 const Element& element,
222 const FVElementGeometry& fvGeometry,
223 const ElementVolumeVariables& elemVolVars,
224 const SubControlVolume &scv) const
225 {
226 NumEqVector source(0.0);
227
228 // add contributions from volume flux sources
229 source += problem.source(element, fvGeometry, elemVolVars, scv);
230
231 // add contribution from possible point sources
232 source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
233
234 return source;
235 }
236
251 NumEqVector computeFlux(const Problem& problem,
252 const Element& element,
253 const FVElementGeometry& fvGeometry,
254 const ElementVolumeVariables& elemVolVars,
255 const SubControlVolumeFace& scvf,
256 const ElementFluxVariablesCache& elemFluxVarsCache) const
257 {
258 DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
259 }
260
261 // \}
262
267 // \{
268
285 const Problem& problem,
286 const Element& element,
287 const FVElementGeometry& fvGeometry,
288 const ElementVolumeVariables& prevElemVolVars,
289 const ElementVolumeVariables& curElemVolVars,
290 const SubControlVolume& scv) const
291 {
292 const auto& curVolVars = curElemVolVars[scv];
293 const auto& prevVolVars = prevElemVolVars[scv];
294
295 // mass balance within the element. this is the
296 // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
297 // euler as time discretization.
298 //
299 // TODO: We might need a more explicit way for
300 // doing the time discretization...
301
303 NumEqVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
304 NumEqVector storage = asImp().computeStorage(problem, scv, curVolVars);
305
306 prevStorage *= prevVolVars.extrusionFactor();
307 storage *= curVolVars.extrusionFactor();
308
309 storage -= prevStorage;
310 storage *= scv.volume();
311 storage /= timeLoop_->timeStepSize();
312
313 residual[scv.localDofIndex()] += storage;
314 }
315
330 const Problem& problem,
331 const Element& element,
332 const FVElementGeometry& fvGeometry,
333 const ElementVolumeVariables& curElemVolVars,
334 const SubControlVolume& scv) const
335 {
337 const auto& curVolVars = curElemVolVars[scv];
338 NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
339 source *= scv.volume()*curVolVars.extrusionFactor();
340
342 residual[scv.localDofIndex()] -= source;
343 }
344
362 const Problem& problem,
363 const Element& element,
364 const FVElementGeometry& fvGeometry,
365 const ElementVolumeVariables& elemVolVars,
366 const ElementBoundaryTypes& elemBcTypes,
367 const ElementFluxVariablesCache& elemFluxVarsCache,
368 const SubControlVolumeFace& scvf) const {}
369
383 NumEqVector evalFlux(const Problem& problem,
384 const Element& element,
385 const FVElementGeometry& fvGeometry,
386 const ElementVolumeVariables& elemVolVars,
387 const ElementFluxVariablesCache& elemFluxVarsCache,
388 const SubControlVolumeFace& scvf) const
389 {
390 return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
391 }
392
393 //\}
394
398 // \{
399
401 template<class PartialDerivativeMatrix>
402 void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
403 const Problem& problem,
404 const Element& element,
405 const FVElementGeometry& fvGeometry,
406 const VolumeVariables& curVolVars,
407 const SubControlVolume& scv) const
408 {
409 DUNE_THROW(Dune::NotImplemented, "analytic storage derivative");
410 }
411
413 template<class PartialDerivativeMatrix>
414 void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
415 const Problem& problem,
416 const Element& element,
417 const FVElementGeometry& fvGeometry,
418 const VolumeVariables& curVolVars,
419 const SubControlVolume& scv) const
420 {
421 DUNE_THROW(Dune::NotImplemented, "analytic source derivative");
422 }
423
425 template<class PartialDerivativeMatrices, class T = TypeTag>
426 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod != DiscretizationMethod::box, void>
427 addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
428 const Problem& problem,
429 const Element& element,
430 const FVElementGeometry& fvGeometry,
431 const ElementVolumeVariables& curElemVolVars,
432 const ElementFluxVariablesCache& elemFluxVarsCache,
433 const SubControlVolumeFace& scvf) const
434 {
435 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for cell-centered models");
436 }
437
439 template<class JacobianMatrix, class T = TypeTag>
440 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethod::box, void>
441 addFluxDerivatives(JacobianMatrix& A,
442 const Problem& problem,
443 const Element& element,
444 const FVElementGeometry& fvGeometry,
445 const ElementVolumeVariables& curElemVolVars,
446 const ElementFluxVariablesCache& elemFluxVarsCache,
447 const SubControlVolumeFace& scvf) const
448 {
449 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for box models");
450 }
451
453 template<class PartialDerivativeMatrices>
454 void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
455 const Problem& problem,
456 const Element& element,
457 const FVElementGeometry& fvGeometry,
458 const ElementVolumeVariables& curElemVolVars,
459 const ElementFluxVariablesCache& elemFluxVarsCache,
460 const SubControlVolumeFace& scvf) const
461 {
462 DUNE_THROW(Dune::NotImplemented, "analytic Dirichlet flux derivative");
463 }
464
466 template<class PartialDerivativeMatrices>
467 void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
468 const Problem& problem,
469 const Element& element,
470 const FVElementGeometry& fvGeometry,
471 const ElementVolumeVariables& curElemVolVars,
472 const ElementFluxVariablesCache& elemFluxVarsCache,
473 const SubControlVolumeFace& scvf) const
474 {
475 DUNE_THROW(Dune::NotImplemented, "analytic Robin flux derivative");
476 }
477
478 //\}
479
483 // \{
484
486 const Problem& problem() const
487 { return *problem_; }
488
491 const TimeLoop& timeLoop() const
492 { return *timeLoop_; }
493
495 bool isStationary() const
496 { return !timeLoop_; }
497
498 // \}
499protected:
500
501 Implementation &asImp()
502 { return *static_cast<Implementation*>(this); }
503
504 const Implementation &asImp() const
505 { return *static_cast<const Implementation*>(this); }
506
507private:
508 const Problem* problem_;
509 const TimeLoop* timeLoop_;
510};
511
512} // end namespace Dumux
513
514#endif
A arithmetic block vector type based on DUNE's reserved vector.
Manages the handling of time dependent problems.
The available discretization methods in Dumux.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
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
The element-wise residual for finite volume schemes.
Definition: fvlocalresidual.hh:45
Implementation & asImp()
Definition: fvlocalresidual.hh:501
NumEqVector computeSource(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolume &scv) const
Calculate the source term of the equation.
Definition: fvlocalresidual.hh:220
void addRobinFluxDerivatives(PartialDerivativeMatrices &derivativeMatrices, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the derivative of Robin type boundary conditions ("solution dependent Neumann")
Definition: fvlocalresidual.hh:467
NumEqVector computeStorage(const Problem &problem, const SubControlVolume &scv, const VolumeVariables &volVars) const
Calculate the source term of the equation.
Definition: fvlocalresidual.hh:200
void evalFlux(ElementResidualVector &residual, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementBoundaryTypes &elemBcTypes, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the source local residual, i.e. the deviation of the source term from zero and add to the res...
Definition: fvlocalresidual.hh:361
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::discMethod !=DiscretizationMethod::box, void > addFluxDerivatives(PartialDerivativeMatrices &derivativeMatrices, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the derivative of the flux residual.
Definition: fvlocalresidual.hh:427
const Problem & problem() const
the problem
Definition: fvlocalresidual.hh:486
void addSourceDerivatives(PartialDerivativeMatrix &partialDerivatives, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const VolumeVariables &curVolVars, const SubControlVolume &scv) const
Compute the derivative of the source residual.
Definition: fvlocalresidual.hh:414
const TimeLoop & timeLoop() const
Definition: fvlocalresidual.hh:491
void addStorageDerivatives(PartialDerivativeMatrix &partialDerivatives, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const VolumeVariables &curVolVars, const SubControlVolume &scv) const
Compute the derivative of the storage residual.
Definition: fvlocalresidual.hh:402
bool isStationary() const
returns true if the residual is stationary
Definition: fvlocalresidual.hh:495
NumEqVector evalFlux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the flux local residual, i.e. the deviation of the flux term from zero.
Definition: fvlocalresidual.hh:383
const Implementation & asImp() const
Definition: fvlocalresidual.hh:504
ElementResidualVector evalStorage(const Problem &problem, const Element &element, const GridGeometry &gridGeometry, const GridVariables &gridVariables, const SolutionVector &sol) const
Compute the storage term for the current solution.
Definition: fvlocalresidual.hh:95
FVLocalResidual(const Problem *problem, const TimeLoop *timeLoop=nullptr)
the constructor
Definition: fvlocalresidual.hh:69
ElementResidualVector evalStorage(const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &prevElemVolVars, const ElementVolumeVariables &curElemVolVars) const
Compute the storage local residual, i.e. the deviation of the storage term from zero for instationary...
Definition: fvlocalresidual.hh:143
void evalSource(ElementResidualVector &residual, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const SubControlVolume &scv) const
Compute the source local residual, i.e. the deviation of the source term from zero.
Definition: fvlocalresidual.hh:329
ElementResidualVector evalFluxAndSource(const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const ElementBoundaryTypes &bcTypes) const
Definition: fvlocalresidual.hh:161
void evalStorage(ElementResidualVector &residual, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &prevElemVolVars, const ElementVolumeVariables &curElemVolVars, const SubControlVolume &scv) const
Compute the storage local residual, i.e. the deviation of the storage term from zero for instationary...
Definition: fvlocalresidual.hh:284
std::enable_if_t< GetPropType< T, Properties::GridGeometry >::discMethod==DiscretizationMethod::box, void > addFluxDerivatives(JacobianMatrix &A, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the derivative of the flux residual for the box method.
Definition: fvlocalresidual.hh:441
void addCCDirichletFluxDerivatives(PartialDerivativeMatrices &derivativeMatrices, const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &curElemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const SubControlVolumeFace &scvf) const
Compute the derivative of the Dirichlet flux residual for cell-centered schemes.
Definition: fvlocalresidual.hh:454
NumEqVector computeFlux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const ElementFluxVariablesCache &elemFluxVarsCache) const
Calculate the source term of the equation.
Definition: fvlocalresidual.hh:251
A arithmetic block vector type based on DUNE's reserved vector.
Definition: reservedblockvector.hh:38
Manages the handling of time dependent problems.
Definition: timeloop.hh:72
virtual Scalar timeStepSize() const =0
Returns the suggested time step length .
The default time loop for instationary simulations.
Definition: timeloop.hh:117
Declares all properties used in Dumux.