pvsmodel.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_PVS_MODEL_HH
29 #define EWOMS_PVS_MODEL_HH
30 
31 #include <opm/material/densead/Math.hpp>
32 
33 #include "pvsproperties.hh"
34 #include "pvslocalresidual.hh"
35 #include "pvsnewtonmethod.hh"
36 #include "pvsprimaryvariables.hh"
37 #include "pvsratevector.hh"
38 #include "pvsboundaryratevector.hh"
41 #include "pvsindices.hh"
42 
49 
50 #include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
51 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
52 #include <opm/common/ErrorMacros.hpp>
53 #include <opm/common/Exceptions.hpp>
54 
55 #include <iostream>
56 #include <sstream>
57 #include <string>
58 #include <vector>
59 
60 namespace Ewoms {
61 template <class TypeTag>
62 class PvsModel;
63 }
64 
65 namespace Ewoms {
66 namespace Properties {
69  VtkPhasePresence,
70  VtkComposition,
71  VtkEnergy,
72  VtkDiffusion));
73 
76  LocalResidual,
78 
81 
84 
87 
90 
93 
96 
99 
102 
103 // set the model to a medium verbosity
104 SET_INT_PROP(PvsModel, PvsVerbosity, 1);
105 
107 SET_BOOL_PROP(PvsModel, EnableEnergy, false);
108 
109 // disable molecular diffusion by default
110 SET_BOOL_PROP(PvsModel, EnableDiffusion, false);
111 
113 SET_SCALAR_PROP(PvsModel, PvsPressureBaseWeight, 1.0);
114 
116 SET_SCALAR_PROP(PvsModel, PvsSaturationsBaseWeight, 1.0);
117 
119 SET_SCALAR_PROP(PvsModel, PvsMoleFractionsBaseWeight, 1.0);
120 } // namespace Properties
121 
217 template <class TypeTag>
218 class PvsModel
219  : public MultiPhaseBaseModel<TypeTag>
220 {
221  typedef MultiPhaseBaseModel<TypeTag> ParentType;
222 
223  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
224  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
225  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
226  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
227 
228  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
229  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
230  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
231  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
232 
233  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
234  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
235  enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
236  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
237 
238  typedef typename GridView::template Codim<0>::Entity Element;
239  typedef typename GridView::template Codim<0>::Iterator ElementIterator;
240 
241  typedef Ewoms::EnergyModule<TypeTag, enableEnergy> EnergyModule;
242 
243 public:
244  PvsModel(Simulator& simulator)
245  : ParentType(simulator)
246  {
247  verbosity_ = EWOMS_GET_PARAM(TypeTag, int, PvsVerbosity);
248  numSwitched_ = 0;
249  }
250 
254  static void registerParameters()
255  {
257 
258  // register runtime parameters of the VTK output modules
261 
262  if (enableDiffusion)
264 
265  if (enableEnergy)
267 
268  EWOMS_REGISTER_PARAM(TypeTag, int, PvsVerbosity,
269  "The verbosity level of the primary variable "
270  "switching model");
271  }
272 
276  static std::string name()
277  { return "pvs"; }
278 
282  std::string primaryVarName(unsigned pvIdx) const
283  {
284  std::string s;
285  if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
286  return s;
287 
288  std::ostringstream oss;
289  if (pvIdx == Indices::pressure0Idx)
290  oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
291  else if (Indices::switch0Idx <= pvIdx
292  && pvIdx < Indices::switch0Idx + numPhases - 1)
293  oss << "switch_" << pvIdx - Indices::switch0Idx;
294  else if (Indices::switch0Idx + numPhases - 1 <= pvIdx
295  && pvIdx < Indices::switch0Idx + numComponents - 1)
296  oss << "auxMoleFrac^" << FluidSystem::componentName(pvIdx);
297  else
298  assert(false);
299 
300  return oss.str();
301  }
302 
306  std::string eqName(unsigned eqIdx) const
307  {
308  std::string s;
309  if (!(s = EnergyModule::eqName(eqIdx)).empty())
310  return s;
311 
312  std::ostringstream oss;
313  if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx
314  + numComponents) {
315  unsigned compIdx = eqIdx - Indices::conti0EqIdx;
316  oss << "continuity^" << FluidSystem::componentName(compIdx);
317  }
318  else
319  assert(false);
320 
321  return oss.str();
322  }
323 
328  {
329  ParentType::updateFailed();
330  numSwitched_ = 0;
331  }
332 
336  void updateBegin()
337  {
338  ParentType::updateBegin();
339 
340  // find the a reference pressure. The first degree of freedom
341  // might correspond to non-interior entities which would lead
342  // to an undefined value, so we have to iterate...
343  size_t nDof = this->numTotalDof();
344  for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
345  if (this->dofTotalVolume(dofIdx) > 0.0) {
346  referencePressure_ =
347  this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
348  if (referencePressure_ > 0.0)
349  break;
350  }
351  }
352  }
353 
357  Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
358  {
359  Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
360  if (tmp > 0)
361  // energy related quantity
362  return tmp;
363 
364  if (Indices::pressure0Idx == pvIdx) {
365  return 10 / referencePressure_;
366  }
367 
368  if (Indices::switch0Idx <= pvIdx && pvIdx < Indices::switch0Idx
369  + numPhases - 1) {
370  unsigned phaseIdx = pvIdx - Indices::switch0Idx;
371 
372  if (!this->solution(/*timeIdx=*/0)[globalDofIdx].phaseIsPresent(phaseIdx))
373  // for saturations, the weight is always 1
374  return 1;
375 
376  // for saturations, the PvsMoleSaturationsBaseWeight
377  // property determines the weight
378  return GET_PROP_VALUE(TypeTag, PvsSaturationsBaseWeight);
379  }
380 
381  // for mole fractions, the PvsMoleFractionsBaseWeight
382  // property determines the weight
383  return GET_PROP_VALUE(TypeTag, PvsMoleFractionsBaseWeight);
384  }
385 
389  Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
390  {
391  Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
392  if (tmp > 0)
393  // energy related equation
394  return tmp;
395 
396  unsigned compIdx = eqIdx - Indices::conti0EqIdx;
397  assert(0 <= compIdx && compIdx <= numComponents);
398 
399  // make all kg equal
400  return FluidSystem::molarMass(compIdx);
401  }
402 
407  {
408  ParentType::advanceTimeLevel();
409  numSwitched_ = 0;
410  }
411 
416  bool switched() const
417  { return numSwitched_ > 0; }
418 
422  template <class DofEntity>
423  void serializeEntity(std::ostream& outstream, const DofEntity& dofEntity)
424  {
425  // write primary variables
426  ParentType::serializeEntity(outstream, dofEntity);
427 
428  unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
429  if (!outstream.good())
430  OPM_THROW(std::runtime_error, "Could not serialize DOF " << dofIdx);
431 
432  outstream << this->solution(/*timeIdx=*/0)[dofIdx].phasePresence() << " ";
433  }
434 
438  template <class DofEntity>
439  void deserializeEntity(std::istream& instream, const DofEntity& dofEntity)
440  {
441  // read primary variables
442  ParentType::deserializeEntity(instream, dofEntity);
443 
444  // read phase presence
445  unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
446  if (!instream.good())
447  OPM_THROW(std::runtime_error,
448  "Could not deserialize DOF " << dofIdx);
449 
450  short tmp;
451  instream >> tmp;
452  this->solution(/*timeIdx=*/0)[dofIdx].setPhasePresence(tmp);
453  this->solution(/*timeIdx=*/1)[dofIdx].setPhasePresence(tmp);
454  }
455 
463  void switchPrimaryVars_()
464  {
465  numSwitched_ = 0;
466 
467  int succeeded;
468  try {
469  std::vector<bool> visited(this->numGridDof(), false);
470  ElementContext elemCtx(this->simulator_);
471 
472  ElementIterator elemIt = this->gridView_.template begin<0>();
473  ElementIterator elemEndIt = this->gridView_.template end<0>();
474  for (; elemIt != elemEndIt; ++elemIt) {
475  const Element& elem = *elemIt;
476  if (elem.partitionType() != Dune::InteriorEntity)
477  continue;
478  elemCtx.updateStencil(elem);
479 
480  size_t numLocalDof = elemCtx.stencil(/*timeIdx=*/0).numPrimaryDof();
481  for (unsigned dofIdx = 0; dofIdx < numLocalDof; ++dofIdx) {
482  unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
483 
484  if (visited[globalIdx])
485  continue;
486  visited[globalIdx] = true;
487 
488  // compute the intensive quantities of the current degree of freedom
489  auto& priVars = this->solution(/*timeIdx=*/0)[globalIdx];
490  elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
491  const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
492 
493  // evaluate primary variable switch
494  short oldPhasePresence = priVars.phasePresence();
495 
496  // set the primary variables and the new phase state
497  // from the current fluid state
498  priVars.assignNaive(intQuants.fluidState());
499 
500  if (oldPhasePresence != priVars.phasePresence()) {
501  if (verbosity_ > 1)
502  printSwitchedPhases_(elemCtx,
503  dofIdx,
504  intQuants.fluidState(),
505  oldPhasePresence,
506  priVars);
507  ++numSwitched_;
508  }
509  }
510  }
511 
512  succeeded = 1;
513  }
514  catch (...)
515  {
516  std::cout << "rank " << this->simulator_.gridView().comm().rank()
517  << " caught an exception during primary variable switching"
518  << "\n" << std::flush;
519  succeeded = 0;
520  }
521  succeeded = this->simulator_.gridView().comm().min(succeeded);
522 
523  if (!succeeded) {
524  OPM_THROW(Opm::NumericalProblem,
525  "A process did not succeed in adapting the primary variables");
526  }
527 
528  // make sure that if there was a variable switch in an
529  // other partition we will also set the switch flag
530  // for our partition.
531  numSwitched_ = this->gridView_.comm().sum(numSwitched_);
532 
533  if (verbosity_ > 0)
534  this->simulator_.model().newtonMethod().endIterMsg()
535  << ", num switched=" << numSwitched_;
536  }
537 
538  template <class FluidState>
539  void printSwitchedPhases_(const ElementContext& elemCtx,
540  unsigned dofIdx,
541  const FluidState& fs,
542  short oldPhasePresence,
543  const PrimaryVariables& newPv) const
544  {
545  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
546 
547  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
548  bool oldPhasePresent = (oldPhasePresence& (1 << phaseIdx)) > 0;
549  bool newPhasePresent = newPv.phaseIsPresent(phaseIdx);
550  if (oldPhasePresent == newPhasePresent)
551  continue;
552 
553  const auto& pos = elemCtx.pos(dofIdx, /*timeIdx=*/0);
554  if (oldPhasePresent && !newPhasePresent) {
555  std::cout << "'" << FluidSystem::phaseName(phaseIdx)
556  << "' phase disappears at position " << pos
557  << ". saturation=" << fs.saturation(phaseIdx)
558  << std::flush;
559  }
560  else {
561  Scalar sumx = 0;
562  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
563  sumx += FsToolbox::value(fs.moleFraction(phaseIdx, compIdx));
564 
565  std::cout << "'" << FluidSystem::phaseName(phaseIdx)
566  << "' phase appears at position " << pos
567  << " sum x = " << sumx << std::flush;
568  }
569  }
570 
571  std::cout << ", new primary variables: ";
572  newPv.print();
573  std::cout << "\n" << std::flush;
574  }
575 
576  void registerOutputModules_()
577  {
578  ParentType::registerOutputModules_();
579 
580  // add the VTK output modules which are meaningful for the model
581  this->addOutputModule(new Ewoms::VtkPhasePresenceModule<TypeTag>(this->simulator_));
582  this->addOutputModule(new Ewoms::VtkCompositionModule<TypeTag>(this->simulator_));
583  if (enableDiffusion)
584  this->addOutputModule(new Ewoms::VtkDiffusionModule<TypeTag>(this->simulator_));
585  if (enableEnergy)
586  this->addOutputModule(new Ewoms::VtkEnergyModule<TypeTag>(this->simulator_));
587  }
588 
589  mutable Scalar referencePressure_;
590 
591  // number of switches of the phase state in the last Newton
592  // iteration
593  unsigned numSwitched_;
594 
595  // verbosity of the model
596  int verbosity_;
597 };
598 }
599 
600 #endif
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
The indices for the compositional multi-phase primary variable switching model.
Definition: pvsindices.hh:45
Element-wise calculation of the local residual for the compositional multi-phase primary variable swi...
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Represents the primary variables used in the primary variable switching compositional model...
Definition: pvsprimaryvariables.hh:59
Implements a vector representing molar, mass or volumetric rates.
Definition: baseauxiliarymodule.hh:37
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkenergymodule.hh:101
bool switched() const
Return true if the primary variables were switched for at least one vertex after the last timestep...
Definition: pvsmodel.hh:416
static std::string name()
Definition: pvsmodel.hh:276
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkcompositionmodule.hh:105
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition: pvsmodel.hh:327
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:59
Contains all data which is required to calculate all fluxes at a flux integration point for the prima...
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void serializeEntity(std::ostream &outstream, const DofEntity &dofEntity)
Write the current solution for a degree of freedom to a restart file.
Definition: pvsmodel.hh:423
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: pvsmodel.hh:357
Implements a vector representing molar, mass or volumetric rates.
Definition: pvsratevector.hh:50
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
VTK output module for quantities which make sense for models which incorperate molecular diffusion...
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
static void registerParameters()
Register all run-time parameters for the PVS compositional model.
Definition: pvsmodel.hh:254
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkphasepresencemodule.hh:80
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
VTK output module for the fluid composition.
Definition: vtkphasepresencemodule.hh:57
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Element-wise calculation of the local residual for the compositional multi-phase primary variable swi...
Definition: pvslocalresidual.hh:47
Contains all data which is required to calculate all fluxes at a flux integration point for the prima...
Definition: pvsextensivequantities.hh:50
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition: pvsmodel.hh:389
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition: vtkdiffusionmodule.hh:101
VTK output module for quantities which make sense for models which assume thermal equilibrium...
VTK output module for the fluid composition.
A generic compositional multi-phase model using primary-variable switching.
Definition: pvsmodel.hh:62
The indices for the compositional multi-phase primary variable switching model.
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
Implements a rate vector on the boundary for the fully implicit compositional multi-phase primary var...
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: pvsmodel.hh:282
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
Represents the primary variables used in the primary variable switching compositional model...
VTK output module for quantities which make sense for models which incorperate molecular diffusion...
Definition: vtkdiffusionmodule.hh:73
A newton solver which is specific to the compositional multi-phase PVS model.
Definition: pvsnewtonmethod.hh:42
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition: pvsintensivequantities.hh:56
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: multiphasebasemodel.hh:147
Declares the properties required for the compositional multi-phase primary variable switching model...
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Definition: multiphasebasemodel.hh:47
Implements a rate vector on the boundary for the fully implicit compositional multi-phase primary var...
Definition: pvsboundaryratevector.hh:46
VTK output module for quantities which make sense for models which assume thermal equilibrium...
Definition: vtkenergymodule.hh:73
VTK output module for the fluid composition.
Definition: vtkcompositionmodule.hh:77
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: pvsmodel.hh:306
Classes required for molecular diffusion.
void deserializeEntity(std::istream &instream, const DofEntity &dofEntity)
Reads the current solution variables for a degree of freedom from a restart file. ...
Definition: pvsmodel.hh:439
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition: pvsmodel.hh:406
A newton solver which is specific to the compositional multi-phase PVS model.
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition: pvsmodel.hh:336