All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
ncplocalresidual.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_NCP_LOCAL_RESIDUAL_HH
29 #define EWOMS_NCP_LOCAL_RESIDUAL_HH
30 
31 #include "ncpproperties.hh"
32 
35 
36 #include <opm/common/Valgrind.hpp>
37 
38 namespace Ewoms {
45 template <class TypeTag>
46 class NcpLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
47 {
48  typedef typename GET_PROP_TYPE(TypeTag, DiscLocalResidual) ParentType;
49  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
50  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
51  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
52  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
53  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
54  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
55  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
56 
57  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
58  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
59  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
60  enum { ncp0EqIdx = Indices::ncp0EqIdx };
61  enum { conti0EqIdx = Indices::conti0EqIdx };
62 
63  enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
65 
66  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
68 
69  typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
70  typedef Dune::BlockVector<EvalEqVector> ElemEvalEqVector;
71  typedef Opm::MathToolbox<Evaluation> Toolbox;
72 
73 public:
77  template <class LhsEval>
78  void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
79  const ElementContext& elemCtx,
80  unsigned dofIdx,
81  unsigned timeIdx,
82  unsigned phaseIdx) const
83  {
84  const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
85  const auto& fluidState = intQuants.fluidState();
86 
87  // compute storage term of all components within all phases
88  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
89  unsigned eqIdx = conti0EqIdx + compIdx;
90  storage[eqIdx] +=
91  Toolbox::template decay<LhsEval>(fluidState.molarity(phaseIdx, compIdx))
92  * Toolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx))
93  * Toolbox::template decay<LhsEval>(intQuants.porosity());
94  }
95 
96  EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
97  }
98 
102  template <class LhsEval>
103  void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
104  const ElementContext& elemCtx,
105  unsigned dofIdx,
106  unsigned timeIdx) const
107  {
108  storage = 0;
109  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
110  addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx);
111 
112  EnergyModule::addSolidHeatStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx));
113  }
114 
118  void computeFlux(RateVector& flux,
119  const ElementContext& elemCtx,
120  unsigned scvfIdx,
121  unsigned timeIdx) const
122  {
123  flux = 0.0;
124  addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
125  Opm::Valgrind::CheckDefined(flux);
126 
127  addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
128  Opm::Valgrind::CheckDefined(flux);
129  }
130 
134  void addAdvectiveFlux(RateVector& flux,
135  const ElementContext& elemCtx,
136  unsigned scvfIdx,
137  unsigned timeIdx) const
138  {
139  const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
140 
141  unsigned focusDofIdx = elemCtx.focusDofIndex();
142  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
143  // data attached to upstream and the downstream DOFs
144  // of the current phase
145  unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
146  const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
147 
148  // this is a bit hacky because it is specific to the element-centered
149  // finite volume scheme. (N.B. that if finite differences are used to
150  // linearize the system of equations, it does not matter.)
151  if (upIdx == focusDofIdx) {
152  Evaluation tmp =
153  up.fluidState().molarDensity(phaseIdx)
154  * extQuants.volumeFlux(phaseIdx);
155 
156  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
157  flux[conti0EqIdx + compIdx] +=
158  tmp*up.fluidState().moleFraction(phaseIdx, compIdx);
159  }
160  }
161  else {
162  Evaluation tmp =
163  Toolbox::value(up.fluidState().molarDensity(phaseIdx))
164  * extQuants.volumeFlux(phaseIdx);
165 
166  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
167  flux[conti0EqIdx + compIdx] +=
168  tmp*Toolbox::value(up.fluidState().moleFraction(phaseIdx, compIdx));
169  }
170  }
171  }
172 
173  EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
174  }
175 
179  void addDiffusiveFlux(RateVector& flux,
180  const ElementContext& elemCtx,
181  unsigned scvfIdx,
182  unsigned timeIdx) const
183  {
184  DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
185  EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
186  }
187 
194  void computeSource(RateVector& source,
195  const ElementContext& elemCtx,
196  unsigned dofIdx,
197  unsigned timeIdx) const
198  {
199  Opm::Valgrind::SetUndefined(source);
200  elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
201  Opm::Valgrind::CheckDefined(source);
202 
203  // evaluate the NCPs (i.e., the "phase presence" equations)
204  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
205  source[ncp0EqIdx + phaseIdx] =
206  phaseNcp(elemCtx, dofIdx, timeIdx, phaseIdx);
207  }
208  }
209 
213  template <class LhsEval = Evaluation>
214  LhsEval phaseNcp(const ElementContext& elemCtx,
215  unsigned dofIdx,
216  unsigned timeIdx,
217  unsigned phaseIdx) const
218  {
219  const auto& fluidState = elemCtx.intensiveQuantities(dofIdx, timeIdx).fluidState();
220  typedef typename std::remove_const<typename std::remove_reference<decltype(fluidState)>::type>::type FluidState;
221 
222  typedef Opm::MathToolbox<LhsEval> LhsToolbox;
223 
224  const LhsEval& a = phaseNotPresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
225  const LhsEval& b = phasePresentIneq_<FluidState, LhsEval>(fluidState, phaseIdx);
226  return LhsToolbox::min(a, b);
227  }
228 
229 private:
234  template <class FluidState, class LhsEval>
235  LhsEval phasePresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
236  {
237  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
238 
239  return FsToolbox::template decay<LhsEval>(fluidState.saturation(phaseIdx));
240  }
241 
246  template <class FluidState, class LhsEval>
247  LhsEval phaseNotPresentIneq_(const FluidState& fluidState, unsigned phaseIdx) const
248  {
249  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
250 
251  // difference of sum of mole fractions in the phase from 100%
252  LhsEval a = 1.0;
253  for (unsigned i = 0; i < numComponents; ++i)
254  a -= FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
255  return a;
256  }
257 };
258 
259 } // namespace Ewoms
260 
261 #endif
void addPhaseStorage(Dune::FieldVector< LhsEval, numEq > &storage, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned phaseIdx) const
Adds the amount all conservation quantities (e.g.
Definition: ncplocalresidual.hh:78
void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume...
Definition: ncplocalresidual.hh:118
Details needed to calculate the local residual in the compositional multi-phase NCP-model ...
Definition: ncplocalresidual.hh:46
void computeSource(RateVector &source, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Calculate the source term of the equation.
Definition: ncplocalresidual.hh:194
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:59
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Declares the properties required for the NCP compositional multi-phase model.
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
void computeStorage(Dune::FieldVector< LhsEval, numEq > &storage, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Evaluate the amount all conservation quantities (e.g.
Definition: ncplocalresidual.hh:103
void addAdvectiveFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Add the advective mass flux at a given flux integration point.
Definition: ncplocalresidual.hh:134
LhsEval phaseNcp(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx, unsigned phaseIdx) const
Returns the value of the NCP-function for a phase.
Definition: ncplocalresidual.hh:214
Classes required for molecular diffusion.
void addDiffusiveFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Adds the diffusive flux at a given flux integration point.
Definition: ncplocalresidual.hh:179