28 #ifndef EWOMS_PVS_MODEL_HH
29 #define EWOMS_PVS_MODEL_HH
31 #include <opm/material/densead/Math.hpp>
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>
61 template <
class TypeTag>
66 namespace Properties {
217 template <
class TypeTag>
219 :
public MultiPhaseBaseModel<TypeTag>
221 typedef MultiPhaseBaseModel<TypeTag> ParentType;
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;
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;
235 enum { enableDiffusion =
GET_PROP_VALUE(TypeTag, EnableDiffusion) };
238 typedef typename GridView::template Codim<0>::Entity Element;
239 typedef typename GridView::template Codim<0>::Iterator ElementIterator;
244 PvsModel(Simulator& simulator)
245 : ParentType(simulator)
269 "The verbosity level of the primary variable "
285 if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
288 std::ostringstream oss;
289 if (pvIdx == Indices::pressure0Idx)
290 oss <<
"pressure_" << FluidSystem::phaseName(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);
309 if (!(s = EnergyModule::eqName(eqIdx)).empty())
312 std::ostringstream oss;
313 if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx
315 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
316 oss <<
"continuity^" << FluidSystem::componentName(compIdx);
329 ParentType::updateFailed();
338 ParentType::updateBegin();
343 size_t nDof = this->numTotalDof();
344 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
345 if (this->dofTotalVolume(dofIdx) > 0.0) {
347 this->solution(0)[dofIdx][Indices::pressure0Idx];
348 if (referencePressure_ > 0.0)
359 Scalar tmp = EnergyModule::primaryVarWeight(*
this, globalDofIdx, pvIdx);
364 if (Indices::pressure0Idx == pvIdx) {
365 return 10 / referencePressure_;
368 if (Indices::switch0Idx <= pvIdx && pvIdx < Indices::switch0Idx
370 unsigned phaseIdx = pvIdx - Indices::switch0Idx;
372 if (!this->solution(0)[globalDofIdx].phaseIsPresent(phaseIdx))
389 Scalar
eqWeight(
unsigned globalDofIdx,
unsigned eqIdx)
const
391 Scalar tmp = EnergyModule::eqWeight(*
this, globalDofIdx, eqIdx);
396 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
397 assert(0 <= compIdx && compIdx <= numComponents);
400 return FluidSystem::molarMass(compIdx);
408 ParentType::advanceTimeLevel();
417 {
return numSwitched_ > 0; }
422 template <
class DofEntity>
426 ParentType::serializeEntity(outstream, dofEntity);
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);
432 outstream << this->solution(0)[dofIdx].phasePresence() <<
" ";
438 template <
class DofEntity>
442 ParentType::deserializeEntity(instream, dofEntity);
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);
452 this->solution(0)[dofIdx].setPhasePresence(tmp);
453 this->solution(1)[dofIdx].setPhasePresence(tmp);
463 void switchPrimaryVars_()
469 std::vector<bool> visited(this->numGridDof(),
false);
470 ElementContext elemCtx(this->simulator_);
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)
478 elemCtx.updateStencil(elem);
480 size_t numLocalDof = elemCtx.stencil(0).numPrimaryDof();
481 for (
unsigned dofIdx = 0; dofIdx < numLocalDof; ++dofIdx) {
482 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
484 if (visited[globalIdx])
486 visited[globalIdx] =
true;
489 auto& priVars = this->solution(0)[globalIdx];
490 elemCtx.updateIntensiveQuantities(priVars, dofIdx, 0);
491 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, 0);
494 short oldPhasePresence = priVars.phasePresence();
498 priVars.assignNaive(intQuants.fluidState());
500 if (oldPhasePresence != priVars.phasePresence()) {
502 printSwitchedPhases_(elemCtx,
504 intQuants.fluidState(),
516 std::cout <<
"rank " << this->simulator_.gridView().comm().rank()
517 <<
" caught an exception during primary variable switching"
518 <<
"\n" << std::flush;
521 succeeded = this->simulator_.gridView().comm().min(succeeded);
524 OPM_THROW(Opm::NumericalProblem,
525 "A process did not succeed in adapting the primary variables");
531 numSwitched_ = this->gridView_.comm().sum(numSwitched_);
534 this->simulator_.model().newtonMethod().endIterMsg()
535 <<
", num switched=" << numSwitched_;
538 template <
class Flu
idState>
539 void printSwitchedPhases_(
const ElementContext& elemCtx,
541 const FluidState& fs,
542 short oldPhasePresence,
543 const PrimaryVariables& newPv)
const
545 typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
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)
553 const auto& pos = elemCtx.pos(dofIdx, 0);
554 if (oldPhasePresent && !newPhasePresent) {
555 std::cout <<
"'" << FluidSystem::phaseName(phaseIdx)
556 <<
"' phase disappears at position " << pos
557 <<
". saturation=" << fs.saturation(phaseIdx)
562 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
563 sumx += FsToolbox::value(fs.moleFraction(phaseIdx, compIdx));
565 std::cout <<
"'" << FluidSystem::phaseName(phaseIdx)
566 <<
"' phase appears at position " << pos
567 <<
" sum x = " << sumx << std::flush;
571 std::cout <<
", new primary variables: ";
573 std::cout <<
"\n" << std::flush;
576 void registerOutputModules_()
578 ParentType::registerOutputModules_();
589 mutable Scalar referencePressure_;
593 unsigned numSwitched_;
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.
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: vtkenergymodule.hh:101
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
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...
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition: pvsmodel.hh:282
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
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...
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition: pvsmodel.hh:306
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...
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
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition: pvsmodel.hh:357
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
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
bool switched() const
Return true if the primary variables were switched for at least one vertex after the last timestep...
Definition: pvsmodel.hh:416
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