All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
powerinjectionproblem.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_POWER_INJECTION_PROBLEM_HH
29 #define EWOMS_POWER_INJECTION_PROBLEM_HH
30 
33 
34 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
35 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
36 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
37 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
38 #include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
39 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
40 #include <opm/material/components/SimpleH2O.hpp>
41 #include <opm/material/components/Air.hpp>
42 #include <opm/common/Unused.hpp>
43 
44 #include <dune/grid/yaspgrid.hh>
45 
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
49 
50 #include <sstream>
51 #include <string>
52 #include <type_traits>
53 #include <iostream>
54 
55 namespace Ewoms {
56 template <class TypeTag>
58 }
59 
60 namespace Ewoms {
61 namespace Properties {
62 NEW_TYPE_TAG(PowerInjectionBaseProblem);
63 
64 // Set the grid implementation to be used
65 SET_TYPE_PROP(PowerInjectionBaseProblem, Grid, Dune::YaspGrid</*dim=*/1>);
66 
67 // set the GridManager property
68 SET_TYPE_PROP(PowerInjectionBaseProblem, GridManager,
70 
71 // Set the problem property
72 SET_TYPE_PROP(PowerInjectionBaseProblem, Problem,
74 
75 // Set the wetting phase
76 SET_PROP(PowerInjectionBaseProblem, WettingPhase)
77 {
78 private:
79  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
80 
81 public:
82  typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
83 };
84 
85 // Set the non-wetting phase
86 SET_PROP(PowerInjectionBaseProblem, NonwettingPhase)
87 {
88 private:
89  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
90 
91 public:
92  typedef Opm::GasPhase<Scalar, Opm::Air<Scalar> > type;
93 };
94 
95 // Set the material Law
96 SET_PROP(PowerInjectionBaseProblem, MaterialLaw)
97 {
98 private:
99  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
100  enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
101  enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
102 
103  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
104  typedef Opm::TwoPhaseMaterialTraits<Scalar,
105  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
106  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>
107  Traits;
108 
109  // define the material law which is parameterized by effective
110  // saturations
111  typedef Opm::RegularizedVanGenuchten<Traits> EffectiveLaw;
112 
113 public:
114  // define the material law parameterized by absolute saturations
115  typedef Opm::EffToAbsLaw<EffectiveLaw> type;
116 };
117 
118 // Write out the filter velocities for this problem
119 SET_BOOL_PROP(PowerInjectionBaseProblem, VtkWriteFilterVelocities, true);
120 
121 // Disable gravity
122 SET_BOOL_PROP(PowerInjectionBaseProblem, EnableGravity, false);
123 
124 // define the properties specific for the power injection problem
125 SET_SCALAR_PROP(PowerInjectionBaseProblem, DomainSizeX, 100.0);
126 SET_SCALAR_PROP(PowerInjectionBaseProblem, DomainSizeY, 1.0);
127 SET_SCALAR_PROP(PowerInjectionBaseProblem, DomainSizeZ, 1.0);
128 
129 SET_INT_PROP(PowerInjectionBaseProblem, CellsX, 250);
130 SET_INT_PROP(PowerInjectionBaseProblem, CellsY, 1);
131 SET_INT_PROP(PowerInjectionBaseProblem, CellsZ, 1);
132 
133 // The default for the end time of the simulation
134 SET_SCALAR_PROP(PowerInjectionBaseProblem, EndTime, 100);
135 
136 // The default for the initial time step size of the simulation
137 SET_SCALAR_PROP(PowerInjectionBaseProblem, InitialTimeStepSize, 1e-3);
138 } // namespace Properties
139 } // namespace Ewoms
140 
141 namespace Ewoms {
154 template <class TypeTag>
155 class PowerInjectionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
156 {
157  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
158 
159  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
160  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
161  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
162  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
163  typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
164  typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
165  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
166  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
167  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
168  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
169  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
170 
171  enum {
172  // number of phases
173 
174  // phase indices
175  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
176  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
177 
178  // equation indices
179  contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
180 
181  // Grid and world dimension
182  dim = GridView::dimension,
183  dimWorld = GridView::dimensionworld
184  };
185 
186  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
187  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
188 
189  typedef typename GridView::ctype CoordScalar;
190  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
191 
192  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
193 
194 public:
199  : ParentType(simulator)
200  { }
201 
205  void finishInit()
206  {
207  ParentType::finishInit();
208 
209  eps_ = 3e-6;
210  FluidSystem::init();
211 
212  temperature_ = 273.15 + 26.6;
213 
214  // parameters for the Van Genuchten law
215  // alpha and n
216  materialParams_.setVgAlpha(0.00045);
217  materialParams_.setVgN(7.3);
218  materialParams_.finalize();
219 
220  K_ = this->toDimMatrix_(5.73e-08); // [m^2]
221 
222  setupInitialFluidState_();
223  }
224 
228 
233  std::string name() const
234  {
235  std::ostringstream oss;
236  oss << "powerinjection_";
237  if (std::is_same<typename GET_PROP_TYPE(TypeTag, FluxModule),
239  oss << "darcy";
240  else
241  oss << "forchheimer";
242 
243  if (std::is_same<typename GET_PROP_TYPE(TypeTag, LocalLinearizerSplice),
244  TTAG(AutoDiffLocalLinearizer)>::value)
245  oss << "_" << "ad";
246  else
247  oss << "_" << "fd";
248 
249  return oss.str();
250  }
251 
255  void endTimeStep()
256  {
257 #ifndef NDEBUG
258  this->model().checkConservativeness();
259 
260  // Calculate storage terms
261  EqVector storage;
262  this->model().globalStorage(storage);
263 
264  // Write mass balance information for rank 0
265  if (this->gridView().comm().rank() == 0) {
266  std::cout << "Storage: " << storage << std::endl << std::flush;
267  }
268 #endif // NDEBUG
269  }
271 
275 
280  template <class Context>
281  const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED,
282  unsigned spaceIdx OPM_UNUSED,
283  unsigned timeIdx OPM_UNUSED) const
284  { return K_; }
285 
289  template <class Context>
290  Scalar ergunCoefficient(const Context& context OPM_UNUSED,
291  unsigned spaceIdx OPM_UNUSED,
292  unsigned timeIdx OPM_UNUSED) const
293  { return 0.3866; }
294 
298  template <class Context>
299  Scalar porosity(const Context& context OPM_UNUSED,
300  unsigned spaceIdx OPM_UNUSED,
301  unsigned timeIdx OPM_UNUSED) const
302  { return 0.558; }
303 
307  template <class Context>
308  const MaterialLawParams&
309  materialLawParams(const Context& context OPM_UNUSED,
310  unsigned spaceIdx OPM_UNUSED,
311  unsigned timeIdx OPM_UNUSED) const
312  { return materialParams_; }
313 
317  template <class Context>
318  Scalar temperature(const Context& context OPM_UNUSED,
319  unsigned spaceIdx OPM_UNUSED,
320  unsigned timeIdx OPM_UNUSED) const
321  { return temperature_; }
322 
324 
328 
336  template <class Context>
337  void boundary(BoundaryRateVector& values,
338  const Context& context,
339  unsigned spaceIdx,
340  unsigned timeIdx) const
341  {
342  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
343 
344  if (onLeftBoundary_(pos)) {
345  RateVector massRate(0.0);
346  massRate = 0.0;
347  massRate[contiNEqIdx] = -1.00; // kg / (m^2 * s)
348 
349  // impose a forced flow boundary
350  values.setMassRate(massRate);
351  }
352  else if (onRightBoundary_(pos))
353  // free flow boundary with initial condition on the right
354  values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
355  else
356  values.setNoFlow();
357  }
358 
360 
364 
369  template <class Context>
370  void initial(PrimaryVariables& values,
371  const Context& context OPM_UNUSED,
372  unsigned spaceIdx OPM_UNUSED,
373  unsigned timeIdx OPM_UNUSED) const
374  {
375  // assign the primary variables
376  values.assignNaive(initialFluidState_);
377  }
378 
385  template <class Context>
386  void source(RateVector& rate,
387  const Context& context OPM_UNUSED,
388  unsigned spaceIdx OPM_UNUSED,
389  unsigned timeIdx OPM_UNUSED) const
390  { rate = Scalar(0.0); }
391 
393 
394 private:
395  bool onLeftBoundary_(const GlobalPosition& pos) const
396  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
397 
398  bool onRightBoundary_(const GlobalPosition& pos) const
399  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
400 
401  void setupInitialFluidState_()
402  {
403  initialFluidState_.setTemperature(temperature_);
404 
405  Scalar Sw = 1.0;
406  initialFluidState_.setSaturation(wettingPhaseIdx, Sw);
407  initialFluidState_.setSaturation(nonWettingPhaseIdx, 1 - Sw);
408 
409  Scalar p = 1e5;
410  initialFluidState_.setPressure(wettingPhaseIdx, p);
411  initialFluidState_.setPressure(nonWettingPhaseIdx, p);
412  }
413 
414  DimMatrix K_;
415  MaterialLawParams materialParams_;
416 
417  Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
418  Scalar temperature_;
419  Scalar eps_;
420 };
421 
422 } // namespace Ewoms
423 
424 #endif
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Specifies a flux module which uses the Darcy relation.
Definition: darcyfluxmodule.hh:65
std::string name() const
The problem name.
Definition: powerinjectionproblem.hh:233
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: powerinjectionproblem.hh:337
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: powerinjectionproblem.hh:386
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the initial value for a control volume.
Definition: powerinjectionproblem.hh:370
Provides a grid manager which a regular grid made of quadrilaterals.
void endTimeStep()
Called by the simulator after each time integration.
Definition: powerinjectionproblem.hh:255
Scalar ergunCoefficient(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Returns the Ergun coefficient.
Definition: powerinjectionproblem.hh:290
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: powerinjectionproblem.hh:318
#define TTAG(TypeTagName)
Convert a type tag name to a type.
Definition: propertysystem.hh:138
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: powerinjectionproblem.hh:299
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: powerinjectionproblem.hh:205
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: powerinjectionproblem.hh:309
Provides a grid manager which a regular grid made of quadrilaterals.
Definition: cubegridmanager.hh:65
A fully-implicit multi-phase flow model which assumes immiscibility of the phases.
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: powerinjectionproblem.hh:281
1D Problem with very fast injection of gas on the left.
Definition: powerinjectionproblem.hh:57
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394