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/material/common/Valgrind.hpp>
37
38#include <dune/istl/bvector.hh>
39#include <dune/grid/common/geometry.hh>
40
41#include <dune/common/fvector.hh>
42
43#include <dune/common/classname.hh>
44
45#include <cmath>
46
47namespace Opm {
56template<class TypeTag>
58{
59private:
61
63 using Element = typename GridView::template Codim<0>::Entity;
64
74
75 static constexpr bool useVolumetricResidual = getPropValue<TypeTag, Properties::UseVolumetricResidual>();
76
77 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
78 enum { extensiveStorageTerm = getPropValue<TypeTag, Properties::ExtensiveStorageTerm>() };
79
80 using Toolbox = MathToolbox<Evaluation>;
81 using EvalVector = Dune::FieldVector<Evaluation, numEq>;
82
83 // copying the local residual class is not a good idea
85 {}
86
87public:
88 using LocalEvalBlockVector = Dune::BlockVector<EvalVector, aligned_allocator<EvalVector, alignof(EvalVector)> >;
89
91 { }
92
94 { }
95
99 static void registerParameters()
100 { }
101
107 { return internalResidual_; }
108
115 const EvalVector& residual(unsigned dofIdx) const
116 { return internalResidual_[dofIdx]; }
117
128 void eval(const Problem& problem, const Element& element)
129 {
130 ElementContext elemCtx(problem);
131 elemCtx.updateAll(element);
132 eval(elemCtx);
133 }
134
144 void eval(ElementContext& elemCtx)
145 {
146 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
147 internalResidual_.resize(numDof);
148 asImp_().eval(internalResidual_, elemCtx);
149 }
150
159 ElementContext& elemCtx) const
160 {
161 assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
162
163 residual = 0.0;
164
165 // evaluate the flux terms
166 asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
167
168 // evaluate the storage and the source terms
169 asImp_().evalVolumeTerms_(residual, elemCtx);
170
171 // evaluate the boundary conditions
172 asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
173
174 if (useVolumetricResidual) {
175 // make the residual volume specific (i.e., make it incorrect mass per cubic
176 // meter instead of total mass)
177 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
178 for (unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) {
179 if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0.0) {
180 // interior DOF
181 Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);
182
183 assert(std::isfinite(dofVolume));
184 Valgrind::CheckDefined(dofVolume);
185
186 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
187 residual[dofIdx][eqIdx] /= dofVolume;
188 }
189 }
190 }
191 }
192
205 const ElementContext& elemCtx,
206 unsigned timeIdx) const
207 {
208 // the derivative of the storage term depends on the current primary variables;
209 // for time indices != 0, the storage term is constant (because these solutions
210 // are not changed by the Newton method!)
211 if (timeIdx == 0) {
212 // calculate the amount of conservation each quantity inside
213 // all primary sub control volumes
214 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
215 for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
216 storage[dofIdx] = 0.0;
217
218 // the volume of the associated DOF
219 Scalar alpha =
220 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
221 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
222
223 // If the degree of freedom which we currently look at is the one at the
224 // center of attention, we need to consider the derivatives for the
225 // storage term, else the storage term is constant w.r.t. the primary
226 // variables of the focused DOF.
227 if (dofIdx == elemCtx.focusDofIndex()) {
228 asImp_().computeStorage(storage[dofIdx],
229 elemCtx,
230 dofIdx,
231 timeIdx);
232
233 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
234 storage[dofIdx][eqIdx] *= alpha;
235 }
236 else {
237 Dune::FieldVector<Scalar, numEq> tmp;
238 asImp_().computeStorage(tmp,
239 elemCtx,
240 dofIdx,
241 timeIdx);
242
243 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
244 storage[dofIdx][eqIdx] = tmp[eqIdx]*alpha;
245 }
246 }
247 }
248 else {
249 // for all previous solutions, the storage term does _not_ depend on the
250 // current primary variables, so we use scalars to store it.
251 if (elemCtx.enableStorageCache()) {
252 size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
253 for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
254 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
255 const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
256 for (unsigned eqIdx=0; eqIdx < numEq; eqIdx++)
257 storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
258 }
259 }
260 else {
261 // calculate the amount of conservation each quantity inside
262 // all primary sub control volumes
263 Dune::FieldVector<Scalar, numEq> tmp;
264 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
265 for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
266 tmp = 0.0;
267 asImp_().computeStorage(tmp,
268 elemCtx,
269 dofIdx,
270 timeIdx);
271 tmp *=
272 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
273 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
274
275 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
276 storage[dofIdx][eqIdx] = tmp[eqIdx];
277 }
278 }
279 }
280
281#ifndef NDEBUG
282 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
283 for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
284 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
285 Valgrind::CheckDefined(storage[dofIdx][eqIdx]);
286 assert(isfinite(storage[dofIdx][eqIdx]));
287 }
288 }
289#endif
290 }
291
300 const ElementContext& elemCtx,
301 unsigned timeIdx) const
302 {
303 RateVector flux;
304
305 const auto& stencil = elemCtx.stencil(timeIdx);
306 // calculate the mass flux over the sub-control volume faces
307 size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
308 for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
309 const auto& face = stencil.interiorFace(scvfIdx);
310 unsigned i = face.interiorIndex();
311 unsigned j = face.exteriorIndex();
312
313 Valgrind::SetUndefined(flux);
314 asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
315 Valgrind::CheckDefined(flux);
316#ifndef NDEBUG
317 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
318 assert(isfinite(flux[eqIdx]));
319#endif
320
321 Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
322 alpha *= face.area();
323 Valgrind::CheckDefined(alpha);
324 assert(alpha > 0.0);
325 assert(isfinite(alpha));
326
327 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
328 flux[eqIdx] *= alpha;
329
330 // The balance equation for a finite volume is given by
331 //
332 // dStorage/dt + Flux = Source
333 //
334 // where the 'Flux' and the 'Source' terms represent the
335 // mass per second which leaves the finite
336 // volume. Re-arranging this, we get
337 //
338 // dStorage/dt + Flux - Source = 0
339 //
340 // Since the mass flux as calculated by computeFlux() goes out of sub-control
341 // volume i and into sub-control volume j, we need to add the flux to finite
342 // volume i and subtract it from finite volume j
343 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
344 assert(isfinite(flux[eqIdx]));
345 residual[i][eqIdx] += flux[eqIdx];
346 residual[j][eqIdx] -= flux[eqIdx];
347 }
348 }
349
350#if !defined NDEBUG
351 // in debug mode, ensure that the residual is well-defined
352 size_t numDof = elemCtx.numDof(timeIdx);
353 for (unsigned i=0; i < numDof; i++) {
354 for (unsigned j = 0; j < numEq; ++ j) {
355 assert(isfinite(residual[i][j]));
356 Valgrind::CheckDefined(residual[i][j]);
357 }
358 }
359#endif
360
361 }
362
364 // The following methods _must_ be overloaded by the actual flow
365 // models!
367
375 void computeStorage(EqVector&,
376 const ElementContext&,
377 unsigned,
378 unsigned) const
379 {
380 throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
381 +" does not implement the required method 'computeStorage()'");
382 }
383
391 void computeFlux(RateVector&,
392 const ElementContext&,
393 unsigned,
394 unsigned) const
395 {
396 throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
397 +" does not implement the required method 'computeFlux()'");
398 }
399
406 void computeSource(RateVector&,
407 const ElementContext&,
408 unsigned,
409 unsigned) const
410 {
411 throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
412 +" does not implement the required method 'computeSource()'");
413 }
414
415protected:
420 const ElementContext& elemCtx,
421 unsigned timeIdx) const
422 {
423 if (!elemCtx.onBoundary())
424 return;
425
426 BoundaryContext boundaryCtx(elemCtx);
427 // move the iterator to the first boundary
428 if(boundaryCtx.intersection(0).neighbor())
429 boundaryCtx.increment();
430
431 // evaluate the boundary for all boundary faces of the current context
432 size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0);
433 for (unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx, boundaryCtx.increment()) {
434 // add the residual of all vertices of the boundary
435 // segment
437 boundaryCtx,
438 faceIdx,
439 timeIdx);
440 }
441
442#if !defined NDEBUG
443 // in debug mode, ensure that the residual and the storage terms are well-defined
444 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
445 for (unsigned i=0; i < numDof; i++) {
446 for (unsigned j = 0; j < numEq; ++ j) {
447 assert(isfinite(residual[i][j]));
448 Valgrind::CheckDefined(residual[i][j]);
449 }
450 }
451#endif
452
453 }
454
460 const BoundaryContext& boundaryCtx,
461 unsigned boundaryFaceIdx,
462 unsigned timeIdx) const
463 {
464 BoundaryRateVector values;
465
466 Valgrind::SetUndefined(values);
467 boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
468 Valgrind::CheckDefined(values);
469
470 const auto& stencil = boundaryCtx.stencil(timeIdx);
471 unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
472 const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
473 for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx) {
474 values[eqIdx] *=
475 stencil.boundaryFace(boundaryFaceIdx).area()
476 * insideIntQuants.extrusionFactor();
477
478 Valgrind::CheckDefined(values[eqIdx]);
479 assert(isfinite(values[eqIdx]));
480 }
481
482 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
483 residual[dofIdx][eqIdx] += values[eqIdx];
484 }
485
492 ElementContext& elemCtx) const
493 {
494 EvalVector tmp;
495 EqVector tmp2;
496 RateVector sourceRate;
497
498 tmp = 0.0;
499 tmp2 = 0.0;
500
501 // evaluate the volumetric terms (storage + source terms)
502 size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
503 for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
504 Scalar extrusionFactor =
505 elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
506 Valgrind::CheckDefined(extrusionFactor);
507 assert(isfinite(extrusionFactor));
508 assert(extrusionFactor > 0.0);
509 Scalar scvVolume =
510 elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
511 Valgrind::CheckDefined(scvVolume);
512 assert(isfinite(scvVolume));
513 assert(scvVolume > 0.0);
514
515 // if the model uses extensive quantities in its storage term, and we use
516 // automatic differention and current DOF is also not the one we currently
517 // focus on, the storage term does not need any derivatives!
518 if (!extensiveStorageTerm &&
519 !std::is_same<Scalar, Evaluation>::value &&
520 dofIdx != elemCtx.focusDofIndex())
521 {
522 asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0);
523 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
524 tmp[eqIdx] = tmp2[eqIdx];
525 }
526 else
527 asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0);
528
529#ifndef NDEBUG
530 Valgrind::CheckDefined(tmp);
531 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
532 assert(isfinite(tmp[eqIdx]));
533#endif
534
535 if (elemCtx.enableStorageCache()) {
536 const auto& model = elemCtx.model();
537 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
538 if (model.newtonMethod().numIterations() == 0 &&
539 !elemCtx.haveStashedIntensiveQuantities())
540 {
541 if (!elemCtx.problem().recycleFirstIterationStorage()) {
542 // we re-calculate the storage term for the solution of the
543 // previous time step from scratch instead of using the one of
544 // the first iteration of the current time step.
545 tmp2 = 0.0;
546 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1);
547 asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
548 }
549 else {
550 // if the storage term is cached and we're in the first iteration
551 // of the time step, use the storage term of the first iteration
552 // as the one as the solution of the last time step (this assumes
553 // that the initial guess for the solution at the end of the time
554 // step is the same as the solution at the beginning of the time
555 // step. This is usually true, but some fancy preprocessing
556 // scheme might invalidate that assumption.)
557 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
558 tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
559 }
560
561 Valgrind::CheckDefined(tmp2);
562
563 model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2);
564 }
565 else {
566 // if the mass storage at the beginning of the time step is not cached,
567 // if the storage term is cached and we're not looking at the first
568 // iteration of the time step, we take the cached data.
569 tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1);
570 Valgrind::CheckDefined(tmp2);
571 }
572 }
573 else {
574 // if the mass storage at the beginning of the time step is not cached,
575 // we re-calculate it from scratch.
576 tmp2 = 0.0;
577 asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
578 Valgrind::CheckDefined(tmp2);
579 }
580
581 // Use the implicit Euler time discretization
582 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
583 double dt = elemCtx.simulator().timeStepSize();
584 assert(dt > 0);
585 tmp[eqIdx] -= tmp2[eqIdx];
586 tmp[eqIdx] *= scvVolume / dt;
587
588 residual[dofIdx][eqIdx] += tmp[eqIdx];
589 }
590
591 Valgrind::CheckDefined(residual[dofIdx]);
592
593 // deal with the source term
594 asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
595
596 // if the model uses extensive quantities in its storage term, and we use
597 // automatic differention and current DOF is also not the one we currently
598 // focus on, the storage term does not need any derivatives!
599 if (!extensiveStorageTerm &&
600 !std::is_same<Scalar, Evaluation>::value &&
601 dofIdx != elemCtx.focusDofIndex())
602 {
603 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
604 residual[dofIdx][eqIdx] -= scalarValue(sourceRate[eqIdx])*scvVolume;
605 }
606 else {
607 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
608 sourceRate[eqIdx] *= scvVolume;
609 residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
610 }
611 }
612
613 Valgrind::CheckDefined(residual[dofIdx]);
614 }
615
616#if !defined NDEBUG
617 // in debug mode, ensure that the residual is well-defined
618 size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
619 for (unsigned i=0; i < numDof; i++) {
620 for (unsigned j = 0; j < numEq; ++ j) {
621 assert(isfinite(residual[i][j]));
622 Valgrind::CheckDefined(residual[i][j]);
623 }
624 }
625#endif
626 }
627
628
629private:
630 Implementation& asImp_()
631 { return *static_cast<Implementation*>(this); }
632
633 const Implementation& asImp_() const
634 { return *static_cast<const Implementation*>(this); }
635
636 LocalEvalBlockVector internalResidual_;
637};
638
639} // namespace Opm
640
641#endif
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1....
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:58
Dune::BlockVector< EvalVector, aligned_allocator< EvalVector, alignof(EvalVector)> > LocalEvalBlockVector
Definition: fvbaselocalresidual.hh:88
FvBaseLocalResidual()
Definition: fvbaselocalresidual.hh:90
void evalVolumeTerms_(LocalEvalBlockVector &residual, 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:491
static void registerParameters()
Register all run-time parameters for the local residual.
Definition: fvbaselocalresidual.hh:99
void computeSource(RateVector &, const ElementContext &, unsigned, unsigned) const
Calculate the source term of the equation.
Definition: fvbaselocalresidual.hh:406
void evalFluxes(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Add the flux term to a local residual.
Definition: fvbaselocalresidual.hh:299
void computeStorage(EqVector &, const ElementContext &, unsigned, unsigned) const
Evaluate the amount all conservation quantities (e.g. phase mass) within a finite sub-control volume.
Definition: fvbaselocalresidual.hh:375
void evalStorage(LocalEvalBlockVector &storage, const ElementContext &elemCtx, unsigned timeIdx) const
Calculate the amount of all conservation quantities stored in all element's sub-control volumes for a...
Definition: fvbaselocalresidual.hh:204
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e. the deviation of the conservation equations from zero and store the ...
Definition: fvbaselocalresidual.hh:128
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:459
void eval(LocalEvalBlockVector &residual, ElementContext &elemCtx) const
Compute the local residual, i.e. the deviation of the conservation equations from zero.
Definition: fvbaselocalresidual.hh:158
void computeFlux(RateVector &, const ElementContext &, unsigned, unsigned) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume.
Definition: fvbaselocalresidual.hh:391
~FvBaseLocalResidual()
Definition: fvbaselocalresidual.hh:93
void evalBoundary_(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Evaluate the boundary conditions of an element.
Definition: fvbaselocalresidual.hh:419
const LocalEvalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:106
void eval(ElementContext &elemCtx)
Compute the local residual, i.e. the deviation of the conservation equations from zero and store the ...
Definition: fvbaselocalresidual.hh:144
const EvalVector & residual(unsigned dofIdx) const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:115
Definition: alignedallocator.hh:114
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilboundaryratevector.hh:37
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:235
This file provides the infrastructure to retrieve run-time parameters.