My Project
Loading...
Searching...
No Matches
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/common/Exceptions.hpp>
32
33#include <opm/material/densead/Math.hpp>
34#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
35#include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
36
40
44
54
55#include <iostream>
56#include <sstream>
57#include <string>
58#include <vector>
59
60namespace Opm {
61template <class TypeTag>
62class PvsModel;
63}
64
65namespace Opm::Properties {
66
67namespace TTag {
68
71{
72 using InheritsFrom = std::tuple<MultiPhaseBaseModel>; };
73} // namespace TTag
74
76template<class TypeTag>
77struct LocalResidual<TypeTag, TTag::PvsModel>
79
81template<class TypeTag>
82struct NewtonMethod<TypeTag, TTag::PvsModel>
84
86template<class TypeTag>
87struct Model<TypeTag, TTag::PvsModel>
88{ using type = Opm::PvsModel<TypeTag>; };
89
91template<class TypeTag>
92struct PrimaryVariables<TypeTag, TTag::PvsModel>
94
96template<class TypeTag>
97struct RateVector<TypeTag, TTag::PvsModel>
99
101template<class TypeTag>
104
106template<class TypeTag>
109
111template<class TypeTag>
114
116template<class TypeTag>
117struct Indices<TypeTag, TTag::PvsModel>
118{ using type = Opm::PvsIndices<TypeTag, /*PVIdx=*/0>; };
119
121template<class TypeTag>
122struct EnableEnergy<TypeTag, TTag::PvsModel>
123{ static constexpr bool value = false; };
124
125// disable molecular diffusion by default
126template<class TypeTag>
127struct EnableDiffusion<TypeTag, TTag::PvsModel>
128{ static constexpr bool value = false; };
129
131template<class TypeTag>
132struct PvsPressureBaseWeight<TypeTag, TTag::PvsModel>
133{
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 1.0;
136};
137
139template<class TypeTag>
140struct PvsSaturationsBaseWeight<TypeTag, TTag::PvsModel>
141{
142 using type = GetPropType<TypeTag, Scalar>;
143 static constexpr type value = 1.0;
144};
145
147template<class TypeTag>
149{
150 using type = GetPropType<TypeTag, Scalar>;
151 static constexpr type value = 1.0;
152};
153
154} // namespace Opm::Properties
155
156namespace Opm::Parameters {
157
159struct PvsVerbosity { static constexpr int value = 1; };
160
161} // namespace Opm::Parameters
162
163namespace Opm {
164
261template <class TypeTag>
263 : public MultiPhaseBaseModel<TypeTag>
264{
266
271
276
278 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
279 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
280 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
281
282 using Element = typename GridView::template Codim<0>::Entity;
283 using ElementIterator = typename GridView::template Codim<0>::Iterator;
284
286
287public:
288 PvsModel(Simulator& simulator)
289 : ParentType(simulator)
290 {
291 verbosity_ = Parameters::Get<Parameters::PvsVerbosity>();
292 numSwitched_ = 0;
293 }
294
298 static void registerParameters()
299 {
301
302 // register runtime parameters of the VTK output modules
305
306 if (enableDiffusion)
308
309 if (enableEnergy)
311
312 Parameters::Register<Parameters::PvsVerbosity>
313 ("The verbosity level of the primary variable "
314 "switching model");
315 }
316
320 static std::string name()
321 { return "pvs"; }
322
326 std::string primaryVarName(unsigned pvIdx) const
327 {
328 std::string s;
329 if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
330 return s;
331
332 std::ostringstream oss;
333 if (pvIdx == Indices::pressure0Idx)
334 oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
335 else if (Indices::switch0Idx <= pvIdx
336 && pvIdx < Indices::switch0Idx + numPhases - 1)
337 oss << "switch_" << pvIdx - Indices::switch0Idx;
338 else if (Indices::switch0Idx + numPhases - 1 <= pvIdx
339 && pvIdx < Indices::switch0Idx + numComponents - 1)
340 oss << "auxMoleFrac^" << FluidSystem::componentName(pvIdx);
341 else
342 assert(false);
343
344 return oss.str();
345 }
346
350 std::string eqName(unsigned eqIdx) const
351 {
352 std::string s;
353 if (!(s = EnergyModule::eqName(eqIdx)).empty())
354 return s;
355
356 std::ostringstream oss;
357 if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx
358 + numComponents) {
359 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
360 oss << "continuity^" << FluidSystem::componentName(compIdx);
361 }
362 else
363 assert(false);
364
365 return oss.str();
366 }
367
372 {
373 ParentType::updateFailed();
374 numSwitched_ = 0;
375 }
376
381 {
382 ParentType::updateBegin();
383
384 // find the a reference pressure. The first degree of freedom
385 // might correspond to non-interior entities which would lead
386 // to an undefined value, so we have to iterate...
387 size_t nDof = this->numTotalDof();
388 for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
389 if (this->dofTotalVolume(dofIdx) > 0.0) {
390 referencePressure_ =
391 this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
392 if (referencePressure_ > 0.0)
393 break;
394 }
395 }
396 }
397
401 Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
402 {
403 Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
404 if (tmp > 0)
405 // energy related quantity
406 return tmp;
407
408 if (Indices::pressure0Idx == pvIdx) {
409 return 10 / referencePressure_;
410 }
411
412 if (Indices::switch0Idx <= pvIdx && pvIdx < Indices::switch0Idx
413 + numPhases - 1) {
414 unsigned phaseIdx = pvIdx - Indices::switch0Idx;
415
416 if (!this->solution(/*timeIdx=*/0)[globalDofIdx].phaseIsPresent(phaseIdx))
417 // for saturations, the weight is always 1
418 return 1;
419
420 // for saturations, the PvsMoleSaturationsBaseWeight
421 // property determines the weight
423 }
424
425 // for mole fractions, the PvsMoleFractionsBaseWeight
426 // property determines the weight
428 }
429
433 Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
434 {
435 Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
436 if (tmp > 0)
437 // energy related equation
438 return tmp;
439
440 unsigned compIdx = eqIdx - Indices::conti0EqIdx;
441 assert(compIdx <= numComponents);
442
443 // make all kg equal
444 return FluidSystem::molarMass(compIdx);
445 }
446
451 {
452 ParentType::advanceTimeLevel();
453 numSwitched_ = 0;
454 }
455
460 bool switched() const
461 { return numSwitched_ > 0; }
462
466 template <class DofEntity>
467 void serializeEntity(std::ostream& outstream, const DofEntity& dofEntity)
468 {
469 // write primary variables
470 ParentType::serializeEntity(outstream, dofEntity);
471
472 unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
473 if (!outstream.good())
474 throw std::runtime_error("Could not serialize DOF "+std::to_string(dofIdx));
475
476 outstream << this->solution(/*timeIdx=*/0)[dofIdx].phasePresence() << " ";
477 }
478
482 template <class DofEntity>
483 void deserializeEntity(std::istream& instream, const DofEntity& dofEntity)
484 {
485 // read primary variables
486 ParentType::deserializeEntity(instream, dofEntity);
487
488 // read phase presence
489 unsigned dofIdx = static_cast<unsigned>(this->dofMapper().index(dofEntity));
490 if (!instream.good())
491 throw std::runtime_error("Could not deserialize DOF "+std::to_string(dofIdx));
492
493 short tmp;
494 instream >> tmp;
495 this->solution(/*timeIdx=*/0)[dofIdx].setPhasePresence(tmp);
496 this->solution(/*timeIdx=*/1)[dofIdx].setPhasePresence(tmp);
497 }
498
506 void switchPrimaryVars_()
507 {
508 numSwitched_ = 0;
509
510 int succeeded;
511 try {
512 std::vector<bool> visited(this->numGridDof(), false);
513 ElementContext elemCtx(this->simulator_);
514
515 for (const auto& elem : elements(this->gridView_, Dune::Partitions::interior)) {
516 elemCtx.updateStencil(elem);
517
518 size_t numLocalDof = elemCtx.stencil(/*timeIdx=*/0).numPrimaryDof();
519 for (unsigned dofIdx = 0; dofIdx < numLocalDof; ++dofIdx) {
520 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
521
522 if (visited[globalIdx])
523 continue;
524 visited[globalIdx] = true;
525
526 // compute the intensive quantities of the current degree of freedom
527 auto& priVars = this->solution(/*timeIdx=*/0)[globalIdx];
528 elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
529 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
530
531 // evaluate primary variable switch
532 short oldPhasePresence = priVars.phasePresence();
533
534 // set the primary variables and the new phase state
535 // from the current fluid state
536 priVars.assignNaive(intQuants.fluidState());
537
538 if (oldPhasePresence != priVars.phasePresence()) {
539 if (verbosity_ > 1)
540 printSwitchedPhases_(elemCtx,
541 dofIdx,
542 intQuants.fluidState(),
544 priVars);
545 ++numSwitched_;
546 }
547 }
548 }
549
550 succeeded = 1;
551 }
552 catch (...)
553 {
554 std::cout << "rank " << this->simulator_.gridView().comm().rank()
555 << " caught an exception during primary variable switching"
556 << "\n" << std::flush;
557 succeeded = 0;
558 }
559 succeeded = this->simulator_.gridView().comm().min(succeeded);
560
561 if (!succeeded)
562 throw NumericalProblem("A process did not succeed in adapting the primary variables");
563
564 // make sure that if there was a variable switch in an
565 // other partition we will also set the switch flag
566 // for our partition.
567 numSwitched_ = this->gridView_.comm().sum(numSwitched_);
568
569 if (verbosity_ > 0)
570 this->simulator_.model().newtonMethod().endIterMsg()
571 << ", num switched=" << numSwitched_;
572 }
573
574 template <class FluidState>
575 void printSwitchedPhases_(const ElementContext& elemCtx,
576 unsigned dofIdx,
577 const FluidState& fs,
578 short oldPhasePresence,
579 const PrimaryVariables& newPv) const
580 {
582
583 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
584 bool oldPhasePresent = (oldPhasePresence& (1 << phaseIdx)) > 0;
585 bool newPhasePresent = newPv.phaseIsPresent(phaseIdx);
587 continue;
588
589 const auto& pos = elemCtx.pos(dofIdx, /*timeIdx=*/0);
591 std::cout << "'" << FluidSystem::phaseName(phaseIdx)
592 << "' phase disappears at position " << pos
593 << ". saturation=" << fs.saturation(phaseIdx)
594 << std::flush;
595 }
596 else {
597 Scalar sumx = 0;
598 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
599 sumx += FsToolbox::value(fs.moleFraction(phaseIdx, compIdx));
600
601 std::cout << "'" << FluidSystem::phaseName(phaseIdx)
602 << "' phase appears at position " << pos
603 << " sum x = " << sumx << std::flush;
604 }
605 }
606
607 std::cout << ", new primary variables: ";
608 newPv.print();
609 std::cout << "\n" << std::flush;
610 }
611
612 void registerOutputModules_()
613 {
614 ParentType::registerOutputModules_();
615
616 // add the VTK output modules which are meaningful for the model
617 this->addOutputModule(new Opm::VtkPhasePresenceModule<TypeTag>(this->simulator_));
618 this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
619 if (enableDiffusion)
620 this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
621 if (enableEnergy)
622 this->addOutputModule(new Opm::VtkEnergyModule<TypeTag>(this->simulator_));
623 }
624
625 mutable Scalar referencePressure_;
626
627 // number of switches of the phase state in the last Newton
628 // iteration
629 unsigned numSwitched_;
630
631 // verbosity of the model
632 int verbosity_;
633};
634}
635
636#endif
Provides the auxiliary methods required for consideration of the energy equation.
Definition energymodule.hh:50
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
Definition multiphasebasemodel.hh:153
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition multiphasebasemodel.hh:179
Implements a rate vector on the boundary for the fully implicit compositional multi-phase primary var...
Definition pvsboundaryratevector.hh:47
Contains all data which is required to calculate all fluxes at a flux integration point for the prima...
Definition pvsextensivequantities.hh:54
The indices for the compositional multi-phase primary variable switching model.
Definition pvsindices.hh:48
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition pvsintensivequantities.hh:61
Element-wise calculation of the local residual for the compositional multi-phase primary variable swi...
Definition pvslocalresidual.hh:48
A generic compositional multi-phase model using primary-variable switching.
Definition pvsmodel.hh:264
void advanceTimeLevel()
Called by the problem if a time integration was successful, post processing of the solution is done a...
Definition pvsmodel.hh:450
void serializeEntity(std::ostream &outstream, const DofEntity &dofEntity)
Write the current solution for a degree of freedom to a restart file.
Definition pvsmodel.hh:467
std::string eqName(unsigned eqIdx) const
Given an equation index, return a human readable name.
Definition pvsmodel.hh:350
std::string primaryVarName(unsigned pvIdx) const
Given an primary variable index, return a human readable name.
Definition pvsmodel.hh:326
void updateBegin()
Called by the update() method before it tries to apply the newton method.
Definition pvsmodel.hh:380
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
Returns the relative weight of an equation.
Definition pvsmodel.hh:433
static std::string name()
Definition pvsmodel.hh:320
bool switched() const
Return true if the primary variables were switched for at least one vertex after the last timestep.
Definition pvsmodel.hh:460
static void registerParameters()
Register all run-time parameters for the PVS compositional model.
Definition pvsmodel.hh:298
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
Returns the relative weight of a primary variable for calculating relative errors.
Definition pvsmodel.hh:401
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:483
void updateFailed()
Called by the update() method if it was unsuccessful.
Definition pvsmodel.hh:371
A newton solver which is specific to the compositional multi-phase PVS model.
Definition pvsnewtonmethod.hh:52
Represents the primary variables used in the primary variable switching compositional model.
Definition pvsprimaryvariables.hh:60
Implements a vector representing molar, mass or volumetric rates.
Definition pvsratevector.hh:53
VTK output module for the fluid composition.
Definition vtkcompositionmodule.hpp:57
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkcompositionmodule.hpp:86
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
Definition vtkdiffusionmodule.hpp:58
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkdiffusionmodule.hpp:87
VTK output module for quantities which make sense for models which assume thermal equilibrium.
Definition vtkenergymodule.hpp:58
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkenergymodule.hpp:86
VTK output module for the fluid composition.
Definition vtkphasepresencemodule.hpp:48
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition vtkphasepresencemodule.hpp:72
Classes required for molecular diffusion.
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
A base class for fully-implicit multi-phase porous-media flow models which assume multiple fluid phas...
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
Implements a rate vector on the boundary for the fully implicit compositional multi-phase primary var...
Contains all data which is required to calculate all fluxes at a flux integration point for the prima...
The indices for the compositional multi-phase primary variable switching model.
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Element-wise calculation of the local residual for the compositional multi-phase primary variable swi...
A newton solver which is specific to the compositional multi-phase PVS model.
Represents the primary variables used in the primary variable switching compositional model.
Declares the properties required for the compositional multi-phase primary variable switching model.
Implements a vector representing molar, mass or volumetric rates.
The verbosity of the model (0 -> do not print anything, 2 -> spam stdout a lot)
Definition pvsmodel.hh:159
Type of object for specifying boundary conditions.
Definition fvbaseproperties.hh:119
Enable diffusive fluxes?
Definition multiphasebaseproperties.hh:79
Specify whether energy should be considered as a conservation quantity or not.
Definition multiphasebaseproperties.hh:76
Data required to calculate a flux over a face.
Definition fvbaseproperties.hh:149
Enumerations used by the model.
Definition multiphasebaseproperties.hh:48
The secondary variables within a sub-control volume.
Definition fvbaseproperties.hh:133
The type of the local residual function.
Definition fvbaseproperties.hh:94
The type of the model.
Definition basicproperties.hh:88
Specifies the type of the actual Newton method.
Definition newtonmethodproperties.hh:32
A vector of primary variables within a sub-control volume.
Definition fvbaseproperties.hh:130
The basis value for the weight of the mole fraction primary variables.
Definition pvsproperties.hh:51
The basis value for the weight of the pressure primary variable.
Definition pvsproperties.hh:45
The basis value for the weight of the saturation primary variables.
Definition pvsproperties.hh:48
Vector containing volumetric or areal rates of quantities.
Definition fvbaseproperties.hh:116
The type tag for the isothermal single phase problems.
Definition pvsmodel.hh:71
VTK output module for the fluid composition.
VTK output module for quantities which make sense for models which incorperate molecular diffusion.
VTK output module for quantities which make sense for models which assume thermal equilibrium.