All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
fvbasediscretization.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_DISCRETIZATION_HH
29 #define EWOMS_FV_BASE_DISCRETIZATION_HH
30 
31 #include <opm/material/densead/Math.hpp>
32 
33 #include "fvbaseproperties.hh"
34 #include "fvbaselinearizer.hh"
37 #include "fvbaselocalresidual.hh"
38 #include "fvbaseelementcontext.hh"
39 #include "fvbaseboundarycontext.hh"
41 #include "fvbaseconstraints.hh"
42 #include "fvbasediscretization.hh"
44 #include "fvbasenewtonmethod.hh"
48 
55 #include <ewoms/common/timer.hh>
57 
58 #include <opm/material/common/MathToolbox.hpp>
59 #include <opm/common/Valgrind.hpp>
60 #include <opm/common/Unused.hpp>
61 #include <opm/common/ErrorMacros.hpp>
62 #include <opm/common/Exceptions.hpp>
63 
64 #include <dune/common/fvector.hh>
65 #include <dune/common/fmatrix.hh>
66 #include <dune/istl/bvector.hh>
67 
68 #if HAVE_DUNE_FEM
69 #include <dune/fem/space/common/adaptmanager.hh>
70 #include <dune/fem/space/common/restrictprolongtuple.hh>
71 #include <dune/fem/function/blockvectorfunction.hh>
72 #include <dune/fem/misc/capabilities.hh>
73 #endif
74 
75 #include <limits>
76 #include <list>
77 #include <sstream>
78 #include <string>
79 #include <vector>
80 
81 namespace Ewoms {
82 template<class TypeTag>
84 
85 namespace Properties {
88 
91  Dune::MultipleCodimMultipleGeomTypeMapper<typename GET_PROP_TYPE(TypeTag, GridView),
92  Dune::MCMGVertexLayout>);
93 
95 SET_TYPE_PROP(FvBaseDiscretization, ElementMapper,
96  Dune::MultipleCodimMultipleGeomTypeMapper<typename GET_PROP_TYPE(TypeTag, GridView),
97  Dune::MCMGElementLayout>);
98 
100 SET_PROP(FvBaseDiscretization, BorderListCreator)
101 {
102  typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
103  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
104 public:
105  typedef Ewoms::Linear::NullBorderListCreator<GridView, DofMapper> type;
106 };
107 
108 SET_TYPE_PROP(FvBaseDiscretization, DiscLocalResidual, Ewoms::FvBaseLocalResidual<TypeTag>);
109 
110 SET_TYPE_PROP(FvBaseDiscretization, DiscIntensiveQuantities, Ewoms::FvBaseIntensiveQuantities<TypeTag>);
111 SET_TYPE_PROP(FvBaseDiscretization, DiscExtensiveQuantities, Ewoms::FvBaseExtensiveQuantities<TypeTag>);
112 
114 SET_TYPE_PROP(FvBaseDiscretization, GradientCalculator, Ewoms::FvBaseGradientCalculator<TypeTag>);
115 
117 SET_PROP(FvBaseDiscretization, JacobianMatrix)
118 {
119 private:
120  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
121  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
122  typedef typename Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
123 public:
124  typedef typename Dune::BCRSMatrix<MatrixBlock> type;
125 };
126 
129 SET_INT_PROP(FvBaseDiscretization, MaxTimeStepDivisions, 10);
130 
134 SET_TYPE_PROP(FvBaseDiscretization, EqVector,
135  Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
136  GET_PROP_VALUE(TypeTag, NumEq)>);
137 
143 SET_TYPE_PROP(FvBaseDiscretization, RateVector,
144  typename GET_PROP_TYPE(TypeTag, EqVector));
145 
149 SET_TYPE_PROP(FvBaseDiscretization, BoundaryRateVector,
150  typename GET_PROP_TYPE(TypeTag, RateVector));
151 
155 SET_TYPE_PROP(FvBaseDiscretization, Constraints, Ewoms::FvBaseConstraints<TypeTag>);
156 
160 SET_TYPE_PROP(FvBaseDiscretization, ElementEqVector,
161  Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, EqVector)>);
162 
166 SET_TYPE_PROP(FvBaseDiscretization, GlobalEqVector,
167  Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, EqVector)>);
168 
172 SET_TYPE_PROP(FvBaseDiscretization, PrimaryVariables, Ewoms::FvBasePrimaryVariables<TypeTag>);
173 
177 SET_TYPE_PROP(FvBaseDiscretization, SolutionVector,
178  Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, PrimaryVariables)>);
179 
185 SET_TYPE_PROP(FvBaseDiscretization, IntensiveQuantities, Ewoms::FvBaseIntensiveQuantities<TypeTag>);
186 
190 SET_TYPE_PROP(FvBaseDiscretization, ElementContext, Ewoms::FvBaseElementContext<TypeTag>);
191 SET_TYPE_PROP(FvBaseDiscretization, BoundaryContext, Ewoms::FvBaseBoundaryContext<TypeTag>);
192 SET_TYPE_PROP(FvBaseDiscretization, ConstraintsContext, Ewoms::FvBaseConstraintsContext<TypeTag>);
193 
197 SET_TYPE_PROP(FvBaseDiscretization, ThreadManager, Ewoms::ThreadManager<TypeTag>);
198 SET_INT_PROP(FvBaseDiscretization, ThreadsPerProcess, 1);
199 SET_BOOL_PROP(FvBaseDiscretization, UseLinearizationLock, true);
200 
204 SET_TYPE_PROP(FvBaseDiscretization, Linearizer, Ewoms::FvBaseLinearizer<TypeTag>);
205 
207 #if 0
208 // requires GCC 4.6 or later to be able call the constexpr function here
209 SET_SCALAR_PROP(FvBaseDiscretization, MaxTimeStepSize, std::numeric_limits<Scalar>::infinity());
210 #else
211 SET_SCALAR_PROP(FvBaseDiscretization, MaxTimeStepSize, 1e100);
212 #endif
213 SET_SCALAR_PROP(FvBaseDiscretization, MinTimeStepSize, 0.0);
215 
217 SET_BOOL_PROP(FvBaseDiscretization, EnableGridAdaptation, false);
218 
220 SET_BOOL_PROP(FvBaseDiscretization, EnableVtkOutput, true);
221 
223 SET_INT_PROP(FvBaseDiscretization, VtkOutputFormat, Dune::VTK::ascii);
224 
225 // disable caching the storage term by default
226 SET_BOOL_PROP(FvBaseDiscretization, EnableStorageCache, false);
227 
228 // disable constraints by default
229 SET_BOOL_PROP(FvBaseDiscretization, EnableConstraints, false);
230 
231 // by default, disable the intensive quantity cache. If the intensive quantities are
232 // relatively cheap to calculate, the cache basically does not yield any performance
233 // impact because of the intensive quantity cache will cause additional pressure on the
234 // CPU caches...
235 SET_BOOL_PROP(FvBaseDiscretization, EnableIntensiveQuantityCache, false);
236 
237 // do not use thermodynamic hints by default. If you enable this, make sure to also
238 // enable the intensive quantity cache above to avoid getting an exception...
239 SET_BOOL_PROP(FvBaseDiscretization, EnableThermodynamicHints, false);
240 
241 // if the deflection of the newton method is large, we do not need to solve the linear
242 // approximation accurately. Assuming that the value for the current solution is quite
243 // close to the final value, a reduction of 3 orders of magnitude in the defect should be
244 // sufficient...
245 SET_SCALAR_PROP(FvBaseDiscretization, LinearSolverTolerance, 1e-3);
246 
248 SET_INT_PROP(FvBaseDiscretization, TimeDiscHistorySize, 2);
249 
252 SET_BOOL_PROP(FvBaseDiscretization, ExtensiveStorageTerm, false);
253 
254 // use volumetric residuals is default
255 SET_BOOL_PROP(FvBaseDiscretization, UseVolumetricResidual, true);
256 
257 } // namespace Properties
258 
264 template<class TypeTag>
265 class FvBaseDiscretization
266 {
267  typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
268  typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
269  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
270  typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
271  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
272  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
273  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
274  typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
275  typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
276  typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
277  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
278  typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
279  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
280  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
281  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
282  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
283  typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
284  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
285  typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
286  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
287  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
288  typedef typename GET_PROP_TYPE(TypeTag, GradientCalculator) GradientCalculator;
289  typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
290  typedef typename GET_PROP_TYPE(TypeTag, DiscBaseOutputModule) DiscBaseOutputModule;
291  typedef typename GET_PROP_TYPE(TypeTag, GridCommHandleFactory) GridCommHandleFactory;
292  typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
293  typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
294 
295  typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) LocalLinearizer;
296  typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
297 
298  enum {
299  numEq = GET_PROP_VALUE(TypeTag, NumEq),
300  historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize),
301  };
302 
303  typedef std::vector<IntensiveQuantities, Ewoms::aligned_allocator<IntensiveQuantities, alignof(IntensiveQuantities)> > IntensiveQuantitiesVector;
304 
305  typedef typename GridView::template Codim<0>::Entity Element;
306  typedef typename GridView::template Codim<0>::Iterator ElementIterator;
307 
308  typedef Opm::MathToolbox<Evaluation> Toolbox;
309  typedef Dune::FieldVector<Evaluation, numEq> VectorBlock;
310  typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
311 
312  typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector;
313 
314  class BlockVectorWrapper
315  {
316  protected:
317  SolutionVector blockVector_;
318  public:
319  BlockVectorWrapper(const std::string& name OPM_UNUSED, const size_t size)
320  : blockVector_(size)
321  {}
322 
323  SolutionVector& blockVector()
324  { return blockVector_; }
325  const SolutionVector& blockVector() const
326  { return blockVector_; }
327  };
328 
329 #if HAVE_DUNE_FEM
330  typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace;
331 
332  // discrete function storing solution data
333  typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables> DiscreteFunction;
334 
335  // problem restriction and prolongation operator for adaptation
336  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
337  typedef typename Problem :: RestrictProlongOperator ProblemRestrictProlongOperator;
338 
339  // discrete function restriction and prolongation operator for adaptation
340  typedef Dune::Fem::RestrictProlongDefault< DiscreteFunction > DiscreteFunctionRestrictProlong;
341  typedef Dune::Fem::RestrictProlongTuple< DiscreteFunctionRestrictProlong, ProblemRestrictProlongOperator > RestrictProlong;
342  // adaptation classes
343  typedef Dune::Fem::AdaptationManager<Grid, RestrictProlong > AdaptationManager;
344 #else
345  typedef BlockVectorWrapper DiscreteFunction;
346  typedef size_t DiscreteFunctionSpace;
347 #endif
348 
349  // copying a discretization object is not a good idea
350  FvBaseDiscretization(const FvBaseDiscretization& );
351 
352 public:
353  // this constructor required to be explicitly specified because
354  // we've defined a constructor above which deletes all implicitly
355  // generated constructors in C++.
356  FvBaseDiscretization(Simulator& simulator)
357  : simulator_(simulator)
358  , gridView_(simulator.gridView())
359  , elementMapper_(gridView_)
360  , vertexMapper_(gridView_)
361  , newtonMethod_(simulator)
362  , localLinearizer_(ThreadManager::maxThreads())
363  , linearizer_(new Linearizer())
364 #if HAVE_DUNE_FEM
365  , space_( simulator.gridManager().gridPart() )
366 #else
367  , space_( asImp_().numGridDof() )
368 #endif
369  , enableGridAdaptation_( EWOMS_GET_PARAM(TypeTag, bool, EnableGridAdaptation) )
370  {
371 #if HAVE_DUNE_FEM
372  if (enableGridAdaptation_ && !Dune::Fem::Capabilities::isLocallyAdaptive<Grid>::v)
373  OPM_THROW(Opm::NotImplemented,
374  "Grid adaptation enabled, but chosen Grid is not capable of adaptivity");
375 #else
376  if (enableGridAdaptation_)
377  OPM_THROW(Opm::NotImplemented,
378  "Grid adaptation currently requires the presence of the dune-fem module");
379 #endif
380  bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
381  if (enableGridAdaptation_ && !isEcfv)
382  OPM_THROW(Opm::NotImplemented,
383  "Grid adaptation currently only works for the element-centered finite "
384  "volume discretization (is: " << Dune::className<Discretization>() << ")");
385 
386  enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
387 
388  size_t numDof = asImp_().numGridDof();
389  for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
390  solution_[timeIdx].reset(new DiscreteFunction("solution", space_));
391 
392  if (storeIntensiveQuantities()) {
393  intensiveQuantityCache_[timeIdx].resize(numDof);
394  intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
395  }
396 
397  if (enableStorageCache_)
398  storageCache_[timeIdx].resize(numDof);
399  }
400 
401  resizeAndResetIntensiveQuantitiesCache_();
402  asImp_().registerOutputModules_();
403  }
404 
405  ~FvBaseDiscretization()
406  {
407  // delete all output modules
408  auto modIt = outputModules_.begin();
409  const auto& modEndIt = outputModules_.end();
410  for (; modIt != modEndIt; ++modIt)
411  delete *modIt;
412 
413  delete linearizer_;
414  }
415 
419  static void registerParameters()
420  {
421  Linearizer::registerParameters();
422  LocalLinearizer::registerParameters();
423  LocalResidual::registerParameters();
424  GradientCalculator::registerParameters();
425  IntensiveQuantities::registerParameters();
426  ExtensiveQuantities::registerParameters();
428 
429  // register runtime parameters of the output modules
431 
432  EWOMS_REGISTER_PARAM(TypeTag, bool, EnableGridAdaptation, "Enable adaptive grid refinement/coarsening");
433  EWOMS_REGISTER_PARAM(TypeTag, bool, EnableVtkOutput, "Global switch for turing on writing VTK files");
434  EWOMS_REGISTER_PARAM(TypeTag, bool, EnableThermodynamicHints, "Enable thermodynamic hints");
435  EWOMS_REGISTER_PARAM(TypeTag, bool, EnableIntensiveQuantityCache, "Turn on caching of intensive quantities");
436  EWOMS_REGISTER_PARAM(TypeTag, bool, EnableStorageCache, "Store previous storage terms and avoid re-calculating them.");
437  }
438 
442  void finishInit()
443  {
444  // initialize the volume of the finite volumes to zero
445  size_t numDof = asImp_().numGridDof();
446  dofTotalVolume_.resize(numDof);
447  std::fill(dofTotalVolume_.begin(), dofTotalVolume_.end(), 0.0);
448 
449  ElementContext elemCtx(simulator_);
450  gridTotalVolume_ = 0.0;
451 
452  // iterate through the grid and evaluate the initial condition
453  ElementIterator elemIt = gridView_.template begin</*codim=*/0>();
454  const ElementIterator& elemEndIt = gridView_.template end</*codim=*/0>();
455  for (; elemIt != elemEndIt; ++elemIt) {
456  const Element& elem = *elemIt;
457  const bool isInteriorElement = elem.partitionType() == Dune::InteriorEntity;
458  // ignore everything which is not in the interior if the
459  // current process' piece of the grid
460  if (!isInteriorElement)
461  continue;
462 
463  // deal with the current element
464  elemCtx.updateStencil(elem);
465  const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
466 
467  // loop over all element vertices, i.e. sub control volumes
468  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++) {
469  // map the local degree of freedom index to the global one
470  unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
471 
472  Scalar dofVolume = stencil.subControlVolume(dofIdx).volume();
473  dofTotalVolume_[globalIdx] += dofVolume;
474  if (isInteriorElement)
475  gridTotalVolume_ += dofVolume;
476  }
477  }
478 
479  // determine which DOFs should be considered to lie fully in the interior of the
480  // local process grid partition: those which do not have a non-zero volume
481  // before taking the peer processes into account...
482  isLocalDof_.resize(numDof);
483  for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx)
484  isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
485 
486  // add the volumes of the DOFs on the process boundaries
487  const auto sumHandle =
488  GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
489  asImp_().dofMapper());
490  gridView_.communicate(*sumHandle,
491  Dune::InteriorBorder_All_Interface,
492  Dune::ForwardCommunication);
493 
494  // sum up the volumes of the grid partitions
495  gridTotalVolume_ = gridView_.comm().sum(gridTotalVolume_);
496 
497  linearizer_->init(simulator_);
498  for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId)
499  localLinearizer_[threadId].init(simulator_);
500 
501  resizeAndResetIntensiveQuantitiesCache_();
502  if (storeIntensiveQuantities()) {
503  // invalidate all cached intensive quantities
504  for (unsigned timeIdx = 0; timeIdx < historySize; ++ timeIdx)
505  invalidateIntensiveQuantitiesCache(timeIdx);
506  }
507  }
508 
512  bool enableGridAdaptation() const
513  { return enableGridAdaptation_; }
514 
520  {
521  // first set the whole domain to zero
522  SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
523  uCur = Scalar(0.0);
524 
525  ElementContext elemCtx(simulator_);
526 
527  // iterate through the grid and evaluate the initial condition
528  ElementIterator elemIt = gridView_.template begin</*codim=*/0>();
529  const ElementIterator& elemEndIt = gridView_.template end</*codim=*/0>();
530  for (; elemIt != elemEndIt; ++elemIt) {
531  const Element& elem = *elemIt;
532  // ignore everything which is not in the interior if the
533  // current process' piece of the grid
534  if (elem.partitionType() != Dune::InteriorEntity)
535  continue;
536 
537  // deal with the current element
538  elemCtx.updateStencil(elem);
539 
540  // loop over all element vertices, i.e. sub control volumes
541  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++)
542  {
543  // map the local degree of freedom index to the global one
544  unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
545 
546  // let the problem do the dirty work of nailing down
547  // the initial solution.
548  simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
549  asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
550  uCur[globalIdx].checkDefined();
551  }
552  }
553 
554  // synchronize the ghost DOFs (if necessary)
555  asImp_().syncOverlap();
556 
557  // also set the solution of the "previous" time steps to the
558  // initial solution.
559  for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx)
560  solution(timeIdx) = solution(/*timeIdx=*/0);
561 
562  simulator_.problem().initialSolutionApplied();
563  }
564 
569  void prefetch(const Element& elem OPM_UNUSED) const
570  {
571  // do nothing by default
572  }
573 
578  { return newtonMethod_; }
579 
583  const NewtonMethod& newtonMethod() const
584  { return newtonMethod_; }
585 
601  const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
602  {
603  if (!enableThermodynamicHints_())
604  return 0;
605 
606  // the intensive quantities cache doubles as thermodynamic hint
607  return cachedIntensiveQuantities(globalIdx, timeIdx);
608  }
609 
621  const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
622  {
623  if (!enableIntensiveQuantitiesCache_() ||
624  !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx])
625  return 0;
626 
627  if (timeIdx > 0 && enableStorageCache_)
628  // with the storage cache enabled, only the intensive quantities for the most
629  // recent time step are cached!
630  return 0;
631 
632  return &intensiveQuantityCache_[timeIdx][globalIdx];
633  }
634 
643  void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
644  unsigned globalIdx,
645  unsigned timeIdx) const
646  {
647  if (!storeIntensiveQuantities())
648  return;
649 
650  intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
651  intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = true;
652  }
653 
662  unsigned timeIdx,
663  bool newValue) const
664  {
665  if (!storeIntensiveQuantities())
666  return;
667 
668  intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue;
669  }
670 
676  void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
677  {
678  if (storeIntensiveQuantities()) {
679  std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
680  intensiveQuantityCacheUpToDate_[timeIdx].end(),
681  /*value=*/false);
682  }
683  }
684 
693  void shiftIntensiveQuantityCache(unsigned numSlots = 1)
694  {
695  if (!storeIntensiveQuantities())
696  return;
697 
698  if (enableStorageCache()) {
699  // if the storage term is cached, the intensive quantities of the previous
700  // time steps do not need to be accessed, and we can thus spare ourselves to
701  // copy the objects for the intensive quantities.
702  return;
703  }
704 
705  assert(numSlots > 0);
706 
707  for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++ timeIdx) {
708  intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
709  intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
710  }
711 
712  // the cache for the most recent time indices do not need to be invalidated
713  // because the solution for them did not change (TODO: that assumes that there is
714  // no post-processing of the solution after a time step! fix it?)
715  }
716 
723  bool enableStorageCache() const
724  { return enableStorageCache_; }
725 
736  const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
737  {
738  assert(enableStorageCache_);
739  return storageCache_[timeIdx][globalIdx];
740  }
741 
753  void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
754  {
755  assert(enableStorageCache_);
756  storageCache_[timeIdx][globalIdx] = value;
757  }
758 
766  Scalar globalResidual(GlobalEqVector& dest,
767  const SolutionVector& u) const
768  {
769  SolutionVector tmp(asImp_().solution(/*timeIdx=*/0));
770  mutableSolution(/*timeIdx=*/0) = u;
771  Scalar res = asImp_().globalResidual(dest);
772  mutableSolution(/*timeIdx=*/0) = tmp;
773  return res;
774  }
775 
782  Scalar globalResidual(GlobalEqVector& dest) const
783  {
784  dest = 0;
785 
786  OmpMutex mutex;
787  ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
788 #ifdef _OPENMP
789 #pragma omp parallel
790 #endif
791  {
792  // Attention: the variables below are thread specific and thus cannot be
793  // moved in front of the #pragma!
794  unsigned threadId = ThreadManager::threadId();
795  ElementContext elemCtx(simulator_);
796  ElementIterator elemIt = threadedElemIt.beginParallel();
797  LocalEvalBlockVector residual, storageTerm;
798 
799  for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
800  const Element& elem = *elemIt;
801  if (elem.partitionType() != Dune::InteriorEntity)
802  continue;
803 
804  elemCtx.updateAll(elem);
805  residual.resize(elemCtx.numDof(/*timeIdx=*/0));
806  storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
807  asImp_().localResidual(threadId).eval(residual, elemCtx);
808 
809  size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
810  ScopedLock addLock(mutex);
811  for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
812  unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
813  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
814  dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
815  }
816  addLock.unlock();
817  }
818  }
819 
820  // add up the residuals on the process borders
821  const auto sumHandle =
822  GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
823  gridView_.communicate(*sumHandle,
824  Dune::InteriorBorder_InteriorBorder_Interface,
825  Dune::ForwardCommunication);
826 
827  // calculate the square norm of the residual. this is not
828  // entirely correct, since the residual for the finite volumes
829  // which are on the boundary are counted once for every
830  // process. As often in life: shit happens (, we don't care)...
831  Scalar result2 = dest.two_norm2();
832  result2 = asImp_().gridView().comm().sum(result2);
833 
834  return std::sqrt(result2);
835  }
836 
843  void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
844  {
845  storage = 0;
846 
847  OmpMutex mutex;
848  ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
849 #ifdef _OPENMP
850 #pragma omp parallel
851 #endif
852  {
853  // Attention: the variables below are thread specific and thus cannot be
854  // moved in front of the #pragma!
855  unsigned threadId = ThreadManager::threadId();
856  ElementContext elemCtx(simulator_);
857  ElementIterator elemIt = threadedElemIt.beginParallel();
858  LocalEvalBlockVector elemStorage;
859 
860  // in this method, we need to disable the storage cache because we want to
861  // evaluate the storage term for other time indices than the most recent one
862  elemCtx.setEnableStorageCache(false);
863 
864  for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
865  const Element& elem = *elemIt;
866  if (elem.partitionType() != Dune::InteriorEntity)
867  continue; // ignore ghost and overlap elements
868 
869  elemCtx.updateStencil(elem);
870  elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
871 
872  size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
873  elemStorage.resize(numPrimaryDof);
874 
875  localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
876 
877  ScopedLock addLock(mutex);
878  for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx)
879  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
880  storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
881  addLock.unlock();
882  }
883  }
884 
885  storage = gridView_.comm().sum(storage);
886  }
887 
895  void checkConservativeness(Scalar OPM_OPTIM_UNUSED tolerance = -1, bool OPM_OPTIM_UNUSED verbose=false) const
896  {
897 #ifndef NDEBUG
898  Scalar totalBoundaryArea(0.0);
899  Scalar totalVolume(0.0);
900  EvalEqVector totalRate(0.0);
901 
902  // take the newton tolerance times the total volume of the grid if we're not
903  // given an explicit tolerance...
904  if (tolerance <= 0) {
905  tolerance =
906  simulator_.model().newtonMethod().tolerance()
907  * simulator_.model().gridTotalVolume()
908  * 1000;
909  }
910 
911  // we assume the implicit Euler time discretization for now...
912  assert(historySize == 2);
913 
914  EqVector storageBeginTimeStep(0.0);
915  globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
916 
917  EqVector storageEndTimeStep(0.0);
918  globalStorage(storageEndTimeStep, /*timeIdx=*/0);
919 
920  // calculate the rate at the boundary and the source rate
921  ElementContext elemCtx(simulator_);
922  elemCtx.setEnableStorageCache(false);
923  auto eIt = simulator_.gridView().template begin</*codim=*/0>();
924  const auto& elemEndIt = simulator_.gridView().template end</*codim=*/0>();
925  for (; eIt != elemEndIt; ++eIt) {
926  if (eIt->partitionType() != Dune::InteriorEntity)
927  continue; // ignore ghost and overlap elements
928 
929  elemCtx.updateAll(*eIt);
930 
931  // handle the boundary terms
932  if (elemCtx.onBoundary()) {
933  BoundaryContext boundaryCtx(elemCtx);
934 
935  for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
936  BoundaryRateVector values;
937  simulator_.problem().boundary(values,
938  boundaryCtx,
939  faceIdx,
940  /*timeIdx=*/0);
941  Opm::Valgrind::CheckDefined(values);
942 
943  unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
944  const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
945 
946  Scalar bfArea =
947  boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0)
948  * insideIntQuants.extrusionFactor();
949 
950  for (unsigned i = 0; i < values.size(); ++i)
951  values[i] *= bfArea;
952 
953  totalBoundaryArea += bfArea;
954  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
955  totalRate[eqIdx] += values[eqIdx];
956  }
957  }
958 
959  // deal with the source terms
960  for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
961  RateVector values;
962  simulator_.problem().source(values,
963  elemCtx,
964  dofIdx,
965  /*timeIdx=*/0);
966  Opm::Valgrind::CheckDefined(values);
967 
968  const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
969  Scalar dofVolume =
970  elemCtx.dofVolume(dofIdx, /*timeIdx=*/0)
971  * intQuants.extrusionFactor();
972  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
973  totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
974  totalVolume += dofVolume;
975  }
976  }
977 
978  // summarize everything over all processes
979  const auto& comm = simulator_.gridView().comm();
980  totalRate = comm.sum(totalRate);
981  totalBoundaryArea = comm.sum(totalBoundaryArea);
982  totalVolume = comm.sum(totalVolume);
983 
984  if (comm.rank() == 0) {
985  EqVector storageRate = storageBeginTimeStep;
986  storageRate -= storageEndTimeStep;
987  storageRate /= simulator_.timeStepSize();
988  if (verbose) {
989  std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
990  std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
991  std::cout << "rate based on storage terms: " << storageRate << "\n";
992  std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
993  std::cout << "difference in rates: ";
994  for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx)
995  std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
996  std::cout << "\n";
997  }
998  for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
999  Scalar eps =
1000  (std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx]))*tolerance;
1001  eps = std::max(tolerance, eps);
1002  assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
1003  }
1004  }
1005 #endif // NDEBUG
1006  }
1007 
1013  Scalar dofTotalVolume(unsigned globalIdx) const
1014  { return dofTotalVolume_[globalIdx]; }
1015 
1021  bool isLocalDof(unsigned globalIdx) const
1022  { return isLocalDof_[globalIdx]; }
1023 
1028  Scalar gridTotalVolume() const
1029  { return gridTotalVolume_; }
1030 
1036  const SolutionVector& solution(unsigned timeIdx) const
1037  {
1038  return solution_[timeIdx]->blockVector();
1039  }
1040 
1044  SolutionVector& solution(unsigned timeIdx)
1045  {
1046  return solution_[timeIdx]->blockVector();
1047  }
1048 
1049  protected:
1053  SolutionVector& mutableSolution(unsigned timeIdx) const
1054  {
1055  return solution_[timeIdx]->blockVector();
1056  }
1057 
1058  public:
1063  const Linearizer& linearizer() const
1064  { return *linearizer_; }
1065 
1070  Linearizer& linearizer()
1071  { return *linearizer_; }
1072 
1081  const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
1082  { return localLinearizer_[openMpThreadId]; }
1086  LocalLinearizer& localLinearizer(unsigned openMpThreadId)
1087  { return localLinearizer_[openMpThreadId]; }
1088 
1092  const LocalResidual& localResidual(unsigned openMpThreadId) const
1093  { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1097  LocalResidual& localResidual(unsigned openMpThreadId)
1098  { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
1099 
1107  Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
1108  {
1109  Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
1110  return 1.0/std::max(absPv, 1.0);
1111  }
1112 
1119  Scalar eqWeight(unsigned globalVertexIdx OPM_UNUSED, unsigned eqIdx OPM_UNUSED) const
1120  { return 1.0; }
1121 
1131  Scalar relativeDofError(unsigned vertexIdx,
1132  const PrimaryVariables& pv1,
1133  const PrimaryVariables& pv2) const
1134  {
1135  Scalar result = 0.0;
1136  for (unsigned j = 0; j < numEq; ++j) {
1137  Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
1138  Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
1139  //Scalar eqErr = std::abs(pv1[j] - pv2[j]);
1140  //eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
1141 
1142  result = std::max(result, eqErr);
1143  }
1144  return result;
1145  }
1146 
1152  bool update()
1153  {
1154  Ewoms::TimerGuard prePostProcessGuard(prePostProcessTimer_);
1155 
1156 #if HAVE_VALGRIND
1157  for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i)
1158  asImp_().solution(/*timeIdx=*/0)[i].checkDefined();
1159 #endif // HAVE_VALGRIND
1160 
1161  // make sure all timers are prestine
1162  prePostProcessTimer_.halt();
1163  linearizeTimer_.halt();
1164  solveTimer_.halt();
1165  updateTimer_.halt();
1166 
1167  prePostProcessTimer_.start();
1168  asImp_().updateBegin();
1169  prePostProcessTimer_.stop();
1170 
1171  bool converged = false;
1172 
1173  try {
1174  converged = newtonMethod_.apply();
1175  }
1176  catch(...) {
1177  prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1178  linearizeTimer_ += newtonMethod_.linearizeTimer();
1179  solveTimer_ += newtonMethod_.solveTimer();
1180  updateTimer_ += newtonMethod_.updateTimer();
1181 
1182  throw;
1183  }
1184 
1185  prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
1186  linearizeTimer_ += newtonMethod_.linearizeTimer();
1187  solveTimer_ += newtonMethod_.solveTimer();
1188  updateTimer_ += newtonMethod_.updateTimer();
1189 
1190  prePostProcessTimer_.start();
1191  if (converged)
1192  asImp_().updateSuccessful();
1193  else
1194  asImp_().updateFailed();
1195  prePostProcessTimer_.stop();
1196 
1197 #if HAVE_VALGRIND
1198  // make sure that the "non-pseudo" primary variables are defined. Note that
1199  // because of the padding, we can't just simply ask valgrind to check the whole
1200  // solution vectors for definedness...
1201  for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
1202  for (size_t eqIdx = 0; eqIdx < numEq; ++eqIdx)
1203  Opm::Valgrind::CheckDefined(asImp_().solution(/*timeIdx=*/0)[i][eqIdx]);
1204  }
1205 #endif // HAVE_VALGRIND
1206 
1207  return converged;
1208  }
1209 
1218  { }
1219 
1226  { }
1227 
1233  { }
1234 
1238  void adaptGrid()
1239  {
1240 #if HAVE_DUNE_FEM
1241  // adapt the grid if enabled and if all dependencies are available
1242  // adaptation is only done if markForGridAdaptation returns true
1243  if (enableGridAdaptation_)
1244  {
1245  // check if problem allows for adaptation and cells were marked
1246  if( simulator_.problem().markForGridAdaptation() )
1247  {
1248  // adapt the grid and load balance if necessary
1249  adaptationManager().adapt();
1250 
1251  // if the grid has potentially changed, we need to re-create the
1252  // supporting data structures.
1253  resetLinearizer();
1254  finishInit();
1255 
1256  // notify the problem that the grid has changed
1257  simulator_.problem().gridChanged();
1258 
1259  // update the entity mappers
1260  elementMapper_.update();
1261  vertexMapper_.update();
1262 
1263  // notify the modules for visualization output
1264  auto outIt = outputModules_.begin();
1265  auto outEndIt = outputModules_.end();
1266  for (; outIt != outEndIt; ++outIt)
1267  (*outIt)->allocBuffers();
1268  }
1269  }
1270 #endif
1271  }
1272 
1279  {
1280  // Reset the current solution to the one of the
1281  // previous time step so that we can start the next
1282  // update at a physically meaningful solution.
1283  solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
1284  invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
1285  }
1286 
1295  {
1296  // at this point we can adapt the grid
1297  asImp_().adaptGrid();
1298 
1299  // make the current solution the previous one.
1300  solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
1301 
1302  // shift the intensive quantities cache by one position in the
1303  // history
1304  asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
1305  }
1306 
1314  template <class Restarter>
1315  void serialize(Restarter& res OPM_UNUSED)
1316  {
1317  OPM_THROW(std::runtime_error,
1318  "Not implemented: The discretization chosen for this problem does not support"
1319  " restart files. (serialize() method unimplemented)");
1320  }
1321 
1329  template <class Restarter>
1330  void deserialize(Restarter& res OPM_UNUSED)
1331  {
1332  OPM_THROW(std::runtime_error,
1333  "Not implemented: The discretization chosen for this problem does not support"
1334  " restart files. (deserialize() method unimplemented)");
1335  }
1336 
1345  template <class DofEntity>
1346  void serializeEntity(std::ostream& outstream,
1347  const DofEntity& dof)
1348  {
1349  unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1350 
1351  // write phase state
1352  if (!outstream.good()) {
1353  OPM_THROW(std::runtime_error, "Could not serialize degree of freedom " << dofIdx);
1354  }
1355 
1356  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1357  outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
1358  }
1359  }
1360 
1369  template <class DofEntity>
1370  void deserializeEntity(std::istream& instream,
1371  const DofEntity& dof)
1372  {
1373  unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
1374 
1375  for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
1376  if (!instream.good())
1377  OPM_THROW(std::runtime_error,
1378  "Could not deserialize degree of freedom " << dofIdx);
1379  instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
1380  }
1381  }
1382 
1386  size_t numGridDof() const
1387  { OPM_THROW(std::logic_error,
1388  "The discretization class must implement the numGridDof() method!"); }
1389 
1393  size_t numAuxiliaryDof() const
1394  {
1395  size_t result = 0;
1396  auto auxModIt = auxEqModules_.begin();
1397  const auto& auxModEndIt = auxEqModules_.end();
1398  for (; auxModIt != auxModEndIt; ++auxModIt)
1399  result += (*auxModIt)->numDofs();
1400 
1401  return result;
1402  }
1403 
1407  size_t numTotalDof() const
1408  { return asImp_().numGridDof() + numAuxiliaryDof(); }
1409 
1414  const DofMapper& dofMapper() const
1415  { OPM_THROW(std::logic_error,
1416  "The discretization class must implement the dofMapper() method!"); }
1417 
1421  const VertexMapper& vertexMapper() const
1422  { return vertexMapper_; }
1423 
1427  const ElementMapper& elementMapper() const
1428  { return elementMapper_; }
1429 
1435  {
1436  delete linearizer_;
1437  linearizer_ = new Linearizer;
1438  linearizer_->init(simulator_);
1439  }
1440 
1444  static std::string discretizationName()
1445  { return ""; }
1446 
1452  std::string primaryVarName(unsigned pvIdx) const
1453  {
1454  std::ostringstream oss;
1455  oss << "primary variable_" << pvIdx;
1456  return oss.str();
1457  }
1458 
1464  std::string eqName(unsigned eqIdx) const
1465  {
1466  std::ostringstream oss;
1467  oss << "equation_" << eqIdx;
1468  return oss.str();
1469  }
1470 
1477  void updatePVWeights(const ElementContext& elemCtx OPM_UNUSED) const
1478  { }
1479 
1484  { outputModules_.push_back(newModule); }
1485 
1494  template <class VtkMultiWriter>
1496  const SolutionVector& u,
1497  const GlobalEqVector& deltaU) const
1498  {
1499  typedef std::vector<double> ScalarBuffer;
1500 
1501  GlobalEqVector globalResid(u.size());
1502  asImp_().globalResidual(globalResid, u);
1503 
1504  // create the required scalar fields
1505  size_t numGridDof = asImp_().numGridDof();
1506 
1507  // global defect of the two auxiliary equations
1508  ScalarBuffer* def[numEq];
1509  ScalarBuffer* delta[numEq];
1510  ScalarBuffer* priVars[numEq];
1511  ScalarBuffer* priVarWeight[numEq];
1512  ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numGridDof);
1513  for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1514  priVars[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1515  priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1516  delta[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1517  def[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
1518  }
1519 
1520  for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx)
1521  {
1522  for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
1523  (*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
1524  (*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
1525  (*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
1526  (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
1527  }
1528 
1529  PrimaryVariables uOld(u[globalIdx]);
1530  PrimaryVariables uNew(uOld);
1531  uNew -= deltaU[globalIdx];
1532  (*relError)[globalIdx] = asImp_().relativeDofError(globalIdx, uOld, uNew);
1533  }
1534 
1535  DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relErr");
1536 
1537  for (unsigned i = 0; i < numEq; ++i) {
1538  std::ostringstream oss;
1539  oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
1540  DiscBaseOutputModule::attachScalarDofData_(writer,
1541  *priVars[i],
1542  oss.str());
1543 
1544  oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
1545  DiscBaseOutputModule::attachScalarDofData_(writer,
1546  *delta[i],
1547  oss.str());
1548 
1549  oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
1550  DiscBaseOutputModule::attachScalarDofData_(writer,
1551  *priVarWeight[i],
1552  oss.str());
1553 
1554  oss.str(""); oss << "defect_" << asImp_().eqName(i);
1555  DiscBaseOutputModule::attachScalarDofData_(writer,
1556  *def[i],
1557  oss.str());
1558  }
1559 
1560  asImp_().prepareOutputFields();
1561  asImp_().appendOutputFields(writer);
1562  }
1563 
1568  void prepareOutputFields() const
1569  {
1570  bool needFullContextUpdate = false;
1571  auto modIt = outputModules_.begin();
1572  const auto& modEndIt = outputModules_.end();
1573  for (; modIt != modEndIt; ++modIt) {
1574  (*modIt)->allocBuffers();
1575  needFullContextUpdate = needFullContextUpdate || (*modIt)->needExtensiveQuantities();
1576  }
1577 
1578  // iterate over grid
1579  ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
1580 #ifdef _OPENMP
1581 #pragma omp parallel
1582 #endif
1583  {
1584  ElementContext elemCtx(simulator_);
1585  ElementIterator elemIt = threadedElemIt.beginParallel();
1586  for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
1587  const auto& elem = *elemIt;
1588  if (elem.partitionType() != Dune::InteriorEntity)
1589  // ignore non-interior entities
1590  continue;
1591 
1592  if (needFullContextUpdate)
1593  elemCtx.updateAll(elem);
1594  else {
1595  elemCtx.updatePrimaryStencil(elem);
1596  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1597  }
1598 
1599  modIt = outputModules_.begin();
1600  for (; modIt != modEndIt; ++modIt)
1601  (*modIt)->processElement(elemCtx);
1602  }
1603  }
1604  }
1605 
1611  {
1612  auto modIt = outputModules_.begin();
1613  const auto& modEndIt = outputModules_.end();
1614  for (; modIt != modEndIt; ++modIt)
1615  (*modIt)->commitBuffers(writer);
1616  }
1617 
1621  const GridView& gridView() const
1622  { return gridView_; }
1623 
1635  void addAuxiliaryModule(std::shared_ptr<BaseAuxiliaryModule<TypeTag> > auxMod)
1636  {
1637  auxMod->setDofOffset(numTotalDof());
1638  auxEqModules_.push_back(auxMod);
1639 
1640  // resize the solutions
1641  if (enableGridAdaptation_
1642  && !std::is_same<DiscreteFunction, BlockVectorWrapper>::value)
1643  {
1644  OPM_THROW(Opm::NotImplemented,
1645  "Problems which require auxiliary modules cannot be used in conjunction "
1646  "with dune-fem");
1647  }
1648 
1649  size_t numDof = numTotalDof();
1650  for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx)
1651  solution(timeIdx).resize(numDof);
1652 
1653  auxMod->applyInitial();
1654  }
1655 
1662  {
1663  auxEqModules_.clear();
1664  linearizer_->eraseMatrix();
1665  newtonMethod_.eraseMatrix();
1666  }
1667 
1671  size_t numAuxiliaryModules() const
1672  { return auxEqModules_.size(); }
1673 
1677  std::shared_ptr<BaseAuxiliaryModule<TypeTag> > auxiliaryModule(unsigned auxEqModIdx)
1678  { return auxEqModules_[auxEqModIdx]; }
1679 
1683  std::shared_ptr<const BaseAuxiliaryModule<TypeTag> > auxiliaryModule(unsigned auxEqModIdx) const
1684  { return auxEqModules_[auxEqModIdx]; }
1685 
1690  { return enableIntensiveQuantitiesCache_() || enableThermodynamicHints_(); }
1691 
1692 #if HAVE_DUNE_FEM
1693  AdaptationManager& adaptationManager()
1694  {
1695  if( ! adaptationManager_ )
1696  {
1697  // create adaptation objects here, because when doing so in constructor
1698  // problem is not yet intialized, aka seg fault
1699  restrictProlong_.reset(
1700  new RestrictProlong( DiscreteFunctionRestrictProlong(*(solution_[/*timeIdx=*/ 0] )),
1701  simulator_.problem().restrictProlongOperator() ) );
1702  adaptationManager_.reset( new AdaptationManager( simulator_.gridManager().grid(), *restrictProlong_ ) );
1703  }
1704  return *adaptationManager_;
1705  }
1706 #endif
1707 
1708  const Ewoms::Timer& prePostProcessTimer() const
1709  { return prePostProcessTimer_; }
1710 
1711  const Ewoms::Timer& linearizeTimer() const
1712  { return linearizeTimer_; }
1713 
1714  const Ewoms::Timer& solveTimer() const
1715  { return solveTimer_; }
1716 
1717  const Ewoms::Timer& updateTimer() const
1718  { return updateTimer_; }
1719 
1720 protected:
1721  void resizeAndResetIntensiveQuantitiesCache_()
1722  {
1723  // allocate the storage cache
1724  if (enableStorageCache()) {
1725  size_t numDof = asImp_().numGridDof();
1726  for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
1727  storageCache_[timeIdx].resize(numDof);
1728  }
1729  }
1730 
1731  // allocate the intensive quantities cache
1732  if (storeIntensiveQuantities()) {
1733  size_t numDof = asImp_().numGridDof();
1734  for(unsigned timeIdx=0; timeIdx<historySize; ++timeIdx) {
1735  intensiveQuantityCache_[timeIdx].resize(numDof);
1736  intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
1737  invalidateIntensiveQuantitiesCache(timeIdx);
1738  }
1739  }
1740  }
1741  template <class Context>
1742  void supplementInitialSolution_(PrimaryVariables& priVars OPM_UNUSED,
1743  const Context& context OPM_UNUSED,
1744  unsigned dofIdx OPM_UNUSED,
1745  unsigned timeIdx OPM_UNUSED)
1746  { }
1747 
1748  static bool enableIntensiveQuantitiesCache_()
1749  { return EWOMS_GET_PARAM(TypeTag, bool, EnableIntensiveQuantityCache); }
1750 
1751  static bool enableThermodynamicHints_()
1752  { return EWOMS_GET_PARAM(TypeTag, bool, EnableThermodynamicHints); }
1753 
1762  {
1763  // add the output modules available on all model
1764  auto *mod = new Ewoms::VtkPrimaryVarsModule<TypeTag>(simulator_);
1765  this->outputModules_.push_back(mod);
1766  }
1767 
1771  LocalResidual& localResidual_()
1772  { return localLinearizer_.localResidual(); }
1773 
1777  bool verbose_() const
1778  { return gridView_.comm().rank() == 0; }
1779 
1780  Implementation& asImp_()
1781  { return *static_cast<Implementation*>(this); }
1782  const Implementation& asImp_() const
1783  { return *static_cast<const Implementation*>(this); }
1784 
1785  // the problem we want to solve. defines the constitutive
1786  // relations, matxerial laws, etc.
1787  Simulator& simulator_;
1788 
1789  // the representation of the spatial domain of the problem
1790  GridView gridView_;
1791 
1792  // the mappers for element and vertex entities to global indices
1793  ElementMapper elementMapper_;
1794  VertexMapper vertexMapper_;
1795 
1796  // a vector with all auxiliary equations to be considered
1797  std::vector<std::shared_ptr<BaseAuxiliaryModule<TypeTag> > > auxEqModules_;
1798 
1799  NewtonMethod newtonMethod_;
1800 
1801  Ewoms::Timer prePostProcessTimer_;
1802  Ewoms::Timer linearizeTimer_;
1803  Ewoms::Timer solveTimer_;
1804  Ewoms::Timer updateTimer_;
1805 
1806  // calculates the local jacobian matrix for a given element
1807  std::vector<LocalLinearizer> localLinearizer_;
1808  // Linearizes the problem at the current time step using the
1809  // local jacobian
1810  Linearizer *linearizer_;
1811 
1812  // cur is the current iterative solution, prev the converged
1813  // solution of the previous time step
1814  mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize];
1815  mutable std::vector<bool> intensiveQuantityCacheUpToDate_[historySize];
1816 
1817  DiscreteFunctionSpace space_;
1818  mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
1819 
1820 #if HAVE_DUNE_FEM
1821  std::unique_ptr< RestrictProlong > restrictProlong_;
1822  std::unique_ptr< AdaptationManager> adaptationManager_;
1823 #endif
1824 
1825 
1826  std::list<BaseOutputModule<TypeTag>*> outputModules_;
1827 
1828  Scalar gridTotalVolume_;
1829  std::vector<Scalar> dofTotalVolume_;
1830  std::vector<bool> isLocalDof_;
1831 
1832  bool enableGridAdaptation_;
1833  mutable GlobalEqVector storageCache_[historySize];
1834  bool enableStorageCache_;
1835 };
1836 } // namespace Ewoms
1837 
1838 #endif
void serializeEntity(std::ostream &outstream, const DofEntity &dof)
Write the current solution for a degree of freedom to a restart file.
Definition: fvbasediscretization.hh:1346
void addOutputModule(BaseOutputModule< TypeTag > *newModule)
Add an module for writing visualization output after a timestep.
Definition: fvbasediscretization.hh:1483
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1.58.
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition: timerguard.hh:40
The base class for all output writers.
Definition: baseoutputwriter.hh:43
Manages the initializing and running of time dependent problems.
const DofMapper & dofMapper() const
Mapper to convert the Dune entities of the discretization&#39;s degrees of freedoms are to indices...
Definition: fvbasediscretization.hh:1414
SolutionVector & mutableSolution(unsigned timeIdx) const
const
Definition: fvbasediscretization.hh:1053
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
The common code for the linearizers of non-linear systems of equations.
Definition: fvbaselinearizer.hh:65
Scalar globalResidual(GlobalEqVector &dest) const
Compute the global residual for the current solution vector.
Definition: fvbasediscretization.hh:782
Simplifies multi-threaded capabilities.
Definition: threadmanager.hh:52
bool enableStorageCache() const
Returns true iff the storage term is cached.
Definition: fvbasediscretization.hh:723
bool verbose_() const
Returns whether messages should be printed.
Definition: fvbasediscretization.hh:1777
Scalar gridTotalVolume() const
Returns the volume of the whole grid which represents the spatial domain.
Definition: fvbasediscretization.hh:1028
ScalarBuffer * allocateManagedScalarBuffer(size_t numEntities)
Allocate a managed buffer for a scalar field.
Definition: vtkmultiwriter.hh:154
void syncOverlap()
Syncronize the values of the primary variables on the degrees of freedom that overlap with the neighb...
Definition: fvbasediscretization.hh:1217
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Linearizer & linearizer()
Returns the object which linearizes the global system of equations at the current solution...
Definition: fvbasediscretization.hh:1070
Base class for specifying auxiliary equations.
Simplifies writing multi-file VTK datasets.
Definition: vtkmultiwriter.hh:63
Definition: locks.hh:54
void finishInit()
Apply the initial conditions to the model.
Definition: fvbasediscretization.hh:442
void prepareOutputFields() const
Prepare the quantities relevant for the current solution to be appended to the output writers...
Definition: fvbasediscretization.hh:1568
void applyInitialSolution()
Applies the initial solution for all degrees of freedom to which the model applies.
Definition: fvbasediscretization.hh:519
Represents all quantities which available on boundary segments.
Definition: fvbaseboundarycontext.hh:45
SolutionVector & solution(unsigned timeIdx)
const
Definition: fvbasediscretization.hh:1044
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition: fvbasediscretization.hh:1225
Scalar globalResidual(GlobalEqVector &dest, const SolutionVector &u) const
Compute the global residual for an arbitrary solution vector.
Definition: fvbasediscretization.hh:766
Provides an STL-iterator like interface to iterate over the enties of a GridView in OpenMP threaded a...
Definition: threadedentityiterator.hh:41
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
size_t numAuxiliaryModules() const
Returns the number of modules for auxiliary equations.
Definition: fvbasediscretization.hh:1671
This is a grid manager which does not create any border list.
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: fvbasediscretization.hh:1107
bool enableGridAdaptation() const
Returns whether the grid ought to be adapted to the solution during the simulation.
Definition: fvbasediscretization.hh:512
Simplifies multi-threaded capabilities.
Provide the properties at a face which make sense indepentently of the conserved quantities.
Definition: fvbaseextensivequantities.hh:44
const EqVector & cachedStorage(unsigned globalIdx, unsigned timeIdx) const
Retrieve an entry of the cache for the storage term.
Definition: fvbasediscretization.hh:736
void appendOutputFields(BaseOutputWriter &writer) const
Append the quantities relevant for the current solution to an output writer.
Definition: fvbasediscretization.hh:1610
bool isLocalDof(unsigned globalIdx) const
Returns if the overlap of the volume ofa degree of freedom is non-zero.
Definition: fvbasediscretization.hh:1021
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Represents all quantities which available for calculating constraints.
void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
Invalidate the whole intensive quantity cache for time index.
Definition: fvbasediscretization.hh:676
void updateSuccessful()
Called by the update() method if it was successful.
Definition: fvbasediscretization.hh:1232
Base class for the model specific class which provides access to all intensive (i.e., volume averaged) quantities.
Definition: fvbaseintensivequantities.hh:45
A Newton method for models using a finite volume discretization.
LocalResidual & localResidual_()
Reference to the local residal object.
Definition: fvbasediscretization.hh:1771
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void registerOutputModules_()
Register all output modules which make sense for the model.
Definition: fvbasediscretization.hh:1761
static unsigned threadId()
Return the index of the current OpenMP thread.
Definition: threadmanager.hh:119
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Base class for the model specific class which provides access to all intensive (i.e., volume averaged) quantities.
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides data handles for parallel communication which operate on DOFs.
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
The base class for writer modules.
Definition: baseoutputmodule.hh:80
This class calculates gradients of arbitrary quantities at flux integration points using the two-poin...
Definition: fvbasegradientcalculator.hh:48
The base class for the finite volume discretization schemes.
Provides an encapsulation to measure the system time.
Definition: timer.hh:48
LocalLinearizer & localLinearizer(unsigned openMpThreadId)
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element...
Definition: fvbasediscretization.hh:1086
void clearAuxiliaryModules()
Causes the list of auxiliary equations to be cleared.
Definition: fvbasediscretization.hh:1661
This class implements an exception-safe scoped lock-keeper.
Definition: locks.hh:65
Represents all quantities which available on boundary segments.
const ElementMapper & elementMapper() const
Returns the mapper for elements to indices.
Definition: fvbasediscretization.hh:1427
void checkConservativeness(Scalar OPM_OPTIM_UNUSED tolerance=-1, bool OPM_OPTIM_UNUSED verbose=false) const
Ensure that the difference between the storage terms of the last and of the current time step is cons...
Definition: fvbasediscretization.hh:895
bool update()
Try to progress the model to the next timestep.
Definition: fvbasediscretization.hh:1152
Represents all quantities which available for calculating constraints.
Definition: fvbaseconstraintscontext.hh:43
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hh:113
size_t numAuxiliaryDof() const
Returns the number of degrees of freedom (DOFs) of the auxiliary equations.
Definition: fvbasediscretization.hh:1393
void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector &value) const
Set an entry of the cache for the storage term.
Definition: fvbasediscretization.hh:753
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
static std::string discretizationName()
Returns a string of discretization&#39;s human-readable name.
Definition: fvbasediscretization.hh:1444
Scalar relativeDofError(unsigned vertexIdx, const PrimaryVariables &pv1, const PrimaryVariables &pv2) const
Returns the relative error between two vectors of primary variables.
Definition: fvbasediscretization.hh:1131
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
Calculates the local residual and its Jacobian for a single element of the grid.
void updatePVWeights(const ElementContext &elemCtx OPM_UNUSED) const
Update the weights of all primary variables within an element given the complete set of intensive qua...
Definition: fvbasediscretization.hh:1477
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:60
Class to specify constraints for a finite volume spatial discretization.
Definition: fvbaseconstraints.hh:48
const GridView & gridView() const
Reference to the grid view of the spatial domain.
Definition: fvbasediscretization.hh:1621
The common code for the linearizers of non-linear systems of equations.
const LocalLinearizer & localLinearizer(unsigned openMpThreadId) const
Returns the local jacobian which calculates the local stiffness matrix for an arbitrary element...
Definition: fvbasediscretization.hh:1081
Scalar eqWeight(unsigned globalVertexIdx OPM_UNUSED, unsigned eqIdx OPM_UNUSED) const
Returns the relative weight of an equation.
Definition: fvbasediscretization.hh:1119
The base class for the finite volume discretization schemes.
Definition: fvbasediscretization.hh:83
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition: newtonmethod.hh:209
void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx, unsigned timeIdx, bool newValue) const
Invalidate the cache for a given intensive quantities object.
Definition: fvbasediscretization.hh:661
std::shared_ptr< BaseAuxiliaryModule< TypeTag > > auxiliaryModule(unsigned auxEqModIdx)
Returns a given module for auxiliary equations.
Definition: fvbasediscretization.hh:1677
void prefetch(const Element &elem OPM_UNUSED) const
Allows to improve the performance by prefetching all data which is associated with a given element...
Definition: fvbasediscretization.hh:569
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:48
const Linearizer & linearizer() const
Returns the operator linearizer for the global jacobian of the problem.
Definition: fvbasediscretization.hh:1063
Provides an encapsulation to measure the system time.
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: fvbasediscretization.hh:1452
NewtonMethod & newtonMethod()
Returns the newton method object.
Definition: fvbasediscretization.hh:577
const LocalResidual & localResidual(unsigned openMpThreadId) const
Returns the object to calculate the local residual function.
Definition: fvbasediscretization.hh:1092
std::shared_ptr< const BaseAuxiliaryModule< TypeTag > > auxiliaryModule(unsigned auxEqModIdx) const
Returns a given module for auxiliary equations.
Definition: fvbasediscretization.hh:1683
This class stores an array of IntensiveQuantities objects, one intensive quantities object for each o...
Definition: fvbaseelementcontext.hh:52
void adaptGrid()
Called by the update() method when the grid should be refined.
Definition: fvbasediscretization.hh:1238
Class to specify constraints for a finite volume spatial discretization.
Scalar dofTotalVolume(unsigned globalIdx) const
Returns the volume of a given control volume.
Definition: fvbasediscretization.hh:1013
const IntensiveQuantities * thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
Return the thermodynamic hint for a entity on the grid at given time.
Definition: fvbasediscretization.hh:601
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: fvbasediscretization.hh:1464
LocalResidual & localResidual(unsigned openMpThreadId)
Returns the object to calculate the local residual function. const
Definition: fvbasediscretization.hh:1097
const SolutionVector & solution(unsigned timeIdx) const
Reference to the solution at a given history index as a block vector.
Definition: fvbasediscretization.hh:1036
Represents the primary variables used by the a model.
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
void addConvergenceVtkFields(VtkMultiWriter &writer, const SolutionVector &u, const GlobalEqVector &deltaU) const
Add the vector fields for analysing the convergence of the newton method to the a VTK writer...
Definition: fvbasediscretization.hh:1495
const IntensiveQuantities * cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
Return the cached intensive quantities for a entity on the grid at given time.
Definition: fvbasediscretization.hh:621
A simple class which makes sure that a timer gets stopped if an exception is thrown.
const VertexMapper & vertexMapper() const
Returns the mapper for vertices to indices.
Definition: fvbasediscretization.hh:1421
Provide the properties at a face which make sense indepentently of the conserved quantities.
VTK output module for the fluid composition.
Definition: vtkprimaryvarsmodule.hh:60
void updateCachedIntensiveQuantities(const IntensiveQuantities &intQuants, unsigned globalIdx, unsigned timeIdx) const
Update the intensive quantity cache for a entity on the grid at given time.
Definition: fvbasediscretization.hh:643
void serialize(Restarter &res OPM_UNUSED)
Serializes the current state of the model.
Definition: fvbasediscretization.hh:1315
size_t numTotalDof() const
Returns the total number of degrees of freedom (i.e., grid plux auxiliary DOFs)
Definition: fvbasediscretization.hh:1407
void addAuxiliaryModule(std::shared_ptr< BaseAuxiliaryModule< TypeTag > > auxMod)
Add a module for an auxiliary equation.
Definition: fvbasediscretization.hh:1635
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkprimaryvarsmodule.hh:84
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition: fvbasediscretization.hh:1294
Base class for specifying auxiliary equations.
Definition: baseauxiliarymodule.hh:59
void deserialize(Restarter &res OPM_UNUSED)
Deserializes the state of the model.
Definition: fvbasediscretization.hh:1330
size_t numGridDof() const
Returns the number of degrees of freedom (DOFs) for the computational grid.
Definition: fvbasediscretization.hh:1386
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
static bool storeIntensiveQuantities()
Returns true if the cache for intensive quantities is enabled.
Definition: fvbasediscretization.hh:1689
void globalStorage(EqVector &storage, unsigned timeIdx=0) const
Compute the integral over the domain of the storage terms of all conservation quantities.
Definition: fvbasediscretization.hh:843
void deserializeEntity(std::istream &instream, const DofEntity &dof)
Reads the current solution variables for a degree of freedom from a restart file. ...
Definition: fvbasediscretization.hh:1370
void resetLinearizer()
Resets the Jacobian matrix linearizer, so that the boundary types can be altered. ...
Definition: fvbasediscretization.hh:1434
const NewtonMethod & newtonMethod() const
Returns the newton method object.
Definition: fvbasediscretization.hh:583
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
static void registerParameters()
Register all run-time parameters for the model.
Definition: fvbasediscretization.hh:419
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition: fvbasediscretization.hh:1278
void shiftIntensiveQuantityCache(unsigned numSlots=1)
Move the intensive quantities for a given time index to the back.
Definition: fvbasediscretization.hh:693