fvbaselocalresidual.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  This file is part of the Open Porous Media project (OPM).
5 
6  OPM 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 2 of the License, or
9  (at your option) any later version.
10 
11  OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
29 #define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
30 
31 #include "fvbaseproperties.hh"
32 
35 
36 #include <opm/common/Valgrind.hpp>
37 #include <opm/common/Unused.hpp>
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
40 
41 #include <dune/istl/bvector.hh>
42 #include <dune/grid/common/geometry.hh>
43 
44 #include <dune/common/fvector.hh>
45 
46 #include <dune/common/classname.hh>
47 
48 #include <cmath>
49 
50 namespace Ewoms {
59 template<class TypeTag>
61 {
62 private:
63  typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
64 
65  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
66  typedef typename GridView::template Codim<0>::Entity Element;
67 
68  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
69  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
70  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
71  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
72  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
73  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
74  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
75  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
76  typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
77 
78  static constexpr bool useVolumetricResidual = GET_PROP_VALUE(TypeTag, UseVolumetricResidual);
79 
80  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
81  enum { extensiveStorageTerm = GET_PROP_VALUE(TypeTag, ExtensiveStorageTerm) };
82 
83  typedef Opm::MathToolbox<Evaluation> Toolbox;
84  typedef Dune::FieldVector<Evaluation, numEq> EvalVector;
85 
86  // copying the local residual class is not a good idea
88  {}
89 
90 public:
91  typedef Dune::BlockVector<EvalVector, Ewoms::aligned_allocator<EvalVector, alignof(EvalVector)> > LocalEvalBlockVector;
92 
94  { }
95 
97  { }
98 
102  static void registerParameters()
103  { }
104 
109  const LocalEvalBlockVector& residual() const
110  { return internalResidual_; }
111 
118  const EvalVector& residual(unsigned dofIdx) const
119  { return internalResidual_[dofIdx]; }
120 
131  void eval(const Problem& problem, const Element& element)
132  {
133  ElementContext elemCtx(problem);
134  elemCtx.updateAll(element);
135  eval(elemCtx);
136  }
137 
147  void eval(const ElementContext& elemCtx)
148  {
149  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
150  internalResidual_.resize(numDof);
151  asImp_().eval(internalResidual_, elemCtx);
152  }
153 
161  void eval(LocalEvalBlockVector& residual,
162  const ElementContext& elemCtx) const
163  {
164  assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
165 
166  residual = 0.0;
167 
168  // evaluate the flux terms
169  asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
170 
171  // evaluate the storage and the source terms
172  asImp_().evalVolumeTerms_(residual, elemCtx);
173 
174  // evaluate the boundary conditions
175  asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
176 
177  if (useVolumetricResidual) {
178  // make the residual volume specific (i.e., make it incorrect mass per cubic
179  // meter instead of total mass)
180  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
181  for (unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) {
182  if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0.0) {
183  // interior DOF
184  Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);
185 
186  assert(std::isfinite(dofVolume));
187  Opm::Valgrind::CheckDefined(dofVolume);
188 
189  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
190  residual[dofIdx][eqIdx] /= dofVolume;
191  }
192  }
193  }
194  }
195 
207  void evalStorage(LocalEvalBlockVector& storage,
208  const ElementContext& elemCtx,
209  unsigned timeIdx) const
210  {
211  // the derivative of the storage term depends on the current primary variables;
212  // for time indices != 0, the storage term is constant (because these solutions
213  // are not changed by the Newton method!)
214  if (timeIdx == 0) {
215  // calculate the amount of conservation each quantity inside
216  // all primary sub control volumes
217  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
218  for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
219  storage[dofIdx] = 0.0;
220 
221  // the volume of the associated DOF
222  Scalar alpha =
223  elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
224  * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
225 
226  // If the degree of freedom which we currently look at is the one at the
227  // center of attention, we need to consider the derivatives for the
228  // storage term, else the storage term is constant w.r.t. the primary
229  // variables of the focused DOF.
230  if (dofIdx == elemCtx.focusDofIndex()) {
231  asImp_().computeStorage(storage[dofIdx],
232  elemCtx,
233  dofIdx,
234  timeIdx);
235 
236  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
237  storage[dofIdx][eqIdx] *= alpha;
238  }
239  else {
240  Dune::FieldVector<Scalar, numEq> tmp;
241  asImp_().computeStorage(tmp,
242  elemCtx,
243  dofIdx,
244  timeIdx);
245 
246  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
247  storage[dofIdx][eqIdx] = tmp[eqIdx]*alpha;
248  }
249  }
250  }
251  else {
252  // for all previous solutions, the storage term does _not_ depend on the
253  // current primary variables, so we use scalars to store it.
254  if (elemCtx.enableStorageCache()) {
255  size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
256  for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
257  unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
258  const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
259  for (unsigned eqIdx=0; eqIdx < numEq; eqIdx++)
260  storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
261  }
262  }
263  else {
264  // calculate the amount of conservation each quantity inside
265  // all primary sub control volumes
266  Dune::FieldVector<Scalar, numEq> tmp;
267  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
268  for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
269  tmp = 0.0;
270  asImp_().computeStorage(tmp,
271  elemCtx,
272  dofIdx,
273  timeIdx);
274  tmp *=
275  elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
276  * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
277 
278  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
279  storage[dofIdx][eqIdx] = tmp[eqIdx];
280  }
281  }
282  }
283  }
284 
292  void evalFluxes(LocalEvalBlockVector& residual,
293  const ElementContext& elemCtx,
294  unsigned timeIdx) const
295  {
296  RateVector flux;
297 
298  const auto& stencil = elemCtx.stencil(timeIdx);
299  // calculate the mass flux over the sub-control volume faces
300  size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
301  for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
302  const auto& face = stencil.interiorFace(scvfIdx);
303  unsigned i = face.interiorIndex();
304  unsigned j = face.exteriorIndex();
305 
306  Opm::Valgrind::SetUndefined(flux);
307  asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
308  Opm::Valgrind::CheckDefined(flux);
309 
310  Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
311  alpha *= face.area();
312  Opm::Valgrind::CheckDefined(alpha);
313  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
314  flux[eqIdx] *= alpha;
315 
316  // The balance equation for a finite volume is given by
317  //
318  // dStorage/dt + Flux = Source
319  //
320  // where the 'Flux' and the 'Source' terms represent the
321  // mass per second which leaves the finite
322  // volume. Re-arranging this, we get
323  //
324  // dStorage/dt + Flux - Source = 0
325  //
326  // Since the mass flux as calculated by computeFlux() goes out of sub-control
327  // volume i and into sub-control volume j, we need to add the flux to finite
328  // volume i and subtract it from finite volume j
329  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
330  residual[i][eqIdx] += flux[eqIdx];
331  residual[j][eqIdx] -= flux[eqIdx];
332  }
333  }
334 
335 #if !defined NDEBUG
336  // in debug mode, ensure that the residual is well-defined
337  size_t numDof = elemCtx.numDof(timeIdx);
338  for (unsigned i=0; i < numDof; i++) {
339  for (unsigned j = 0; j < numEq; ++ j) {
340  assert(std::isfinite(Toolbox::value(residual[i][j])));
341  Opm::Valgrind::CheckDefined(residual[i][j]);
342  }
343  }
344 #endif
345 
346  }
347 
349  // The following methods _must_ be overloaded by the actual flow
350  // models!
352 
360  void computeStorage(EqVector& storage OPM_UNUSED,
361  const ElementContext& elemCtx OPM_UNUSED,
362  unsigned dofIdx OPM_UNUSED,
363  unsigned timeIdx OPM_UNUSED) const
364  {
365  OPM_THROW(std::logic_error,
366  "Not implemented: The local residual " << Dune::className<Implementation>()
367  << " does not implement the required method 'computeStorage()'");
368  }
369 
377  void computeFlux(RateVector& flux OPM_UNUSED,
378  const ElementContext& elemCtx OPM_UNUSED,
379  unsigned scvfIdx OPM_UNUSED,
380  unsigned timeIdx OPM_UNUSED) const
381  {
382  OPM_THROW(std::logic_error,
383  "Not implemented: The local residual " << Dune::className<Implementation>()
384  << " does not implement the required method 'computeFlux()'");
385  }
386 
393  void computeSource(RateVector& source OPM_UNUSED,
394  const ElementContext& elemCtx OPM_UNUSED,
395  unsigned dofIdx OPM_UNUSED,
396  unsigned timeIdx OPM_UNUSED) const
397  {
398  OPM_THROW(std::logic_error,
399  "Not implemented: The local residual " << Dune::className<Implementation>()
400  << " does not implement the required method 'computeSource()'");
401  }
402 
403 protected:
407  void evalBoundary_(LocalEvalBlockVector& residual,
408  const ElementContext& elemCtx,
409  unsigned timeIdx) const
410  {
411  if (!elemCtx.onBoundary())
412  return;
413 
414  BoundaryContext boundaryCtx(elemCtx);
415 
416  // evaluate the boundary for all boundary faces of the current context
417  size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0);
418  for (unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx) {
419  // add the residual of all vertices of the boundary
420  // segment
422  boundaryCtx,
423  faceIdx,
424  timeIdx);
425  }
426 
427 #if !defined NDEBUG
428  // in debug mode, ensure that the residual and the storage terms are well-defined
429  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
430  for (unsigned i=0; i < numDof; i++) {
431  for (unsigned j = 0; j < numEq; ++ j) {
432  assert(std::isfinite(Toolbox::value(residual[i][j])));
433  Opm::Valgrind::CheckDefined(residual[i][j]);
434  }
435  }
436 #endif
437 
438  }
439 
444  void evalBoundarySegment_(LocalEvalBlockVector& residual,
445  const BoundaryContext& boundaryCtx,
446  unsigned boundaryFaceIdx,
447  unsigned timeIdx) const
448  {
449  BoundaryRateVector values;
450 
451  Opm::Valgrind::SetUndefined(values);
452  boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
453  Opm::Valgrind::CheckDefined(values);
454 
455  const auto& stencil = boundaryCtx.stencil(timeIdx);
456  unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
457  const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
458  for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx)
459  values[eqIdx] *=
460  stencil.boundaryFace(boundaryFaceIdx).area()
461  * insideIntQuants.extrusionFactor();
462 
463  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
464  residual[dofIdx][eqIdx] += values[eqIdx];
465  }
466 
472  void evalVolumeTerms_(LocalEvalBlockVector& residual,
473  const ElementContext& elemCtx) const
474  {
475  EvalVector tmp;
476  EqVector tmp2;
477  RateVector sourceRate;
478 
479  tmp = 0.0;
480  tmp2 = 0.0;
481 
482  // evaluate the volumetric terms (storage + source terms)
483  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
484  for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
485  Scalar extrusionFactor =
486  elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
487  Scalar scvVolume =
488  elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
489  Opm::Valgrind::CheckDefined(scvVolume);
490 
491  // if the model uses extensive quantities in its storage term, and we use
492  // automatic differention and current DOF is also not the one we currently
493  // focus on, the storage term does not need any derivatives!
494  if (!extensiveStorageTerm &&
495  !std::is_same<Scalar, Evaluation>::value &&
496  dofIdx != elemCtx.focusDofIndex())
497  {
498  asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0);
499  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
500  tmp[eqIdx] = tmp2[eqIdx];
501  }
502  else
503  asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0);
504  Opm::Valgrind::CheckDefined(tmp);
505 
506  if (elemCtx.enableStorageCache()) {
507  const auto& model = elemCtx.model();
508  unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
509  if (model.newtonMethod().numIterations() == 0 &&
510  !elemCtx.haveStashedIntensiveQuantities())
511  {
512  // if the storage term is cached and we're in the first iteration of
513  // the time step, update the cache of the storage term (this assumes
514  // that the initial guess for the solution at the end of the time
515  // step is the same as the solution at the beginning of the time
516  // step. This is usually true, but some fancy preprocessing scheme
517  // might invalidate that assumption.)
518  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
519  tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
520  Opm::Valgrind::CheckDefined(tmp2);
521 
522  model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2);
523  }
524  else {
525  // if the storage term is cached and we're not looking at the first
526  // iteration of the time step, we take the cached data.
527  tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1);
528  Opm::Valgrind::CheckDefined(tmp2);
529  }
530  }
531  else {
532  // if the mass storage at the beginning of the time step is not cached,
533  // we re-calculate it from scratch.
534  tmp2 = 0.0;
535  asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
536  Opm::Valgrind::CheckDefined(tmp2);
537  }
538 
539  // Use the implicit Euler time discretization
540  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
541  tmp[eqIdx] -= tmp2[eqIdx];
542  tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
543 
544  residual[dofIdx][eqIdx] += tmp[eqIdx];
545  }
546 
547  Opm::Valgrind::CheckDefined(residual[dofIdx]);
548 
549  // deal with the source term
550  asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
551 
552  // if the model uses extensive quantities in its storage term, and we use
553  // automatic differention and current DOF is also not the one we currently
554  // focus on, the storage term does not need any derivatives!
555  if (!extensiveStorageTerm &&
556  !std::is_same<Scalar, Evaluation>::value &&
557  dofIdx != elemCtx.focusDofIndex())
558  {
559  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
560  residual[dofIdx][eqIdx] -= Opm::scalarValue(sourceRate[eqIdx])*scvVolume;
561  }
562  else {
563  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
564  sourceRate[eqIdx] *= scvVolume;
565  residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
566  }
567  }
568 
569  Opm::Valgrind::CheckDefined(residual[dofIdx]);
570  }
571 
572 #if !defined NDEBUG
573  // in debug mode, ensure that the residual is well-defined
574  size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
575  for (unsigned i=0; i < numDof; i++) {
576  for (unsigned j = 0; j < numEq; ++ j) {
577  assert(Opm::isfinite(residual[i][j]));
578  Opm::Valgrind::CheckDefined(residual[i][j]);
579  }
580  }
581 #endif
582  }
583 
584 
585 private:
586  Implementation& asImp_()
587  { return *static_cast<Implementation*>(this); }
588 
589  const Implementation& asImp_() const
590  { return *static_cast<const Implementation*>(this); }
591 
592  LocalEvalBlockVector internalResidual_;
593 };
594 
595 } // namespace Ewoms
596 
597 #endif
void eval(LocalEvalBlockVector &residual, const ElementContext &elemCtx) const
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:161
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1...
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:131
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
This file provides the infrastructure to retrieve run-time parameters.
Declare the properties used by the infrastructure code of the finite volume discretizations.
void evalFluxes(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Add the flux term to a local residual.
Definition: fvbaselocalresidual.hh:292
const LocalEvalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:109
void computeFlux(RateVector &flux OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume...
Definition: fvbaselocalresidual.hh:377
void computeStorage(EqVector &storage OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the amount all conservation quantities (e.g.
Definition: fvbaselocalresidual.hh:360
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:60
const EvalVector & residual(unsigned dofIdx) const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:118
void evalBoundary_(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Evaluate the boundary conditions of an element.
Definition: fvbaselocalresidual.hh:407
void evalStorage(LocalEvalBlockVector &storage, const ElementContext &elemCtx, unsigned timeIdx) const
Calculate the amount of all conservation quantities stored in all element&#39;s sub-control volumes for a...
Definition: fvbaselocalresidual.hh:207
void evalVolumeTerms_(LocalEvalBlockVector &residual, const ElementContext &elemCtx) const
Add the change in the storage terms and the source term to the local residual of all sub-control volu...
Definition: fvbaselocalresidual.hh:472
void computeSource(RateVector &source OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Calculate the source term of the equation.
Definition: fvbaselocalresidual.hh:393
void evalBoundarySegment_(LocalEvalBlockVector &residual, const BoundaryContext &boundaryCtx, unsigned boundaryFaceIdx, unsigned timeIdx) const
Evaluate all boundary conditions for a single sub-control volume face to the local residual...
Definition: fvbaselocalresidual.hh:444
void eval(const ElementContext &elemCtx)
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:147
static void registerParameters()
Register all run-time parameters for the local residual.
Definition: fvbaselocalresidual.hh:102