My Project
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 
31 #include <opm/models/immiscible/immisciblemodel.hh>
32 #include <opm/models/io/cubegridvanguard.hh>
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 
43 #include <dune/grid/yaspgrid.hh>
44 
45 #include <dune/common/version.hh>
46 #include <dune/common/fvector.hh>
47 #include <dune/common/fmatrix.hh>
48 
49 #include <sstream>
50 #include <string>
51 #include <type_traits>
52 #include <iostream>
53 
54 namespace Opm {
55 template <class TypeTag>
56 class PowerInjectionProblem;
57 }
58 
59 namespace Opm::Properties {
60 
61 namespace TTag {
63 }
64 
65 // Set the grid implementation to be used
66 template<class TypeTag>
67 struct Grid<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Dune::YaspGrid</*dim=*/1>; };
68 
69 // set the Vanguard property
70 template<class TypeTag>
71 struct Vanguard<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Opm::CubeGridVanguard<TypeTag>; };
72 
73 // Set the problem property
74 template<class TypeTag>
75 struct Problem<TypeTag, TTag::PowerInjectionBaseProblem> { using type = Opm::PowerInjectionProblem<TypeTag>; };
76 
77 // Set the wetting phase
78 template<class TypeTag>
79 struct WettingPhase<TypeTag, TTag::PowerInjectionBaseProblem>
80 {
81 private:
82  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
83 
84 public:
85  using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
86 };
87 
88 // Set the non-wetting phase
89 template<class TypeTag>
90 struct NonwettingPhase<TypeTag, TTag::PowerInjectionBaseProblem>
91 {
92 private:
93  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
94 
95 public:
96  using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
97 };
98 
99 // Set the material Law
100 template<class TypeTag>
101 struct MaterialLaw<TypeTag, TTag::PowerInjectionBaseProblem>
102 {
103 private:
104  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
105  enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
106  enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
107 
108  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
109  using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
110  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
111  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
112 
113  // define the material law which is parameterized by effective
114  // saturations
115  using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
116 
117 public:
118  // define the material law parameterized by absolute saturations
119  using type = Opm::EffToAbsLaw<EffectiveLaw>;
120 };
121 
122 // Write out the filter velocities for this problem
123 template<class TypeTag>
124 struct VtkWriteFilterVelocities<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr bool value = true; };
125 
126 // Disable gravity
127 template<class TypeTag>
128 struct EnableGravity<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr bool value = false; };
129 
130 // define the properties specific for the power injection problem
131 template<class TypeTag>
132 struct DomainSizeX<TypeTag, TTag::PowerInjectionBaseProblem>
133 {
134  using type = GetPropType<TypeTag, Scalar>;
135  static constexpr type value = 100.0;
136 };
137 template<class TypeTag>
138 struct DomainSizeY<TypeTag, TTag::PowerInjectionBaseProblem>
139 {
140  using type = GetPropType<TypeTag, Scalar>;
141  static constexpr type value = 1.0;
142 };
143 template<class TypeTag>
144 struct DomainSizeZ<TypeTag, TTag::PowerInjectionBaseProblem>
145 {
146  using type = GetPropType<TypeTag, Scalar>;
147  static constexpr type value = 1.0;
148 };
149 
150 template<class TypeTag>
151 struct CellsX<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 250; };
152 template<class TypeTag>
153 struct CellsY<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 1; };
154 template<class TypeTag>
155 struct CellsZ<TypeTag, TTag::PowerInjectionBaseProblem> { static constexpr int value = 1; };
156 
157 // The default for the end time of the simulation
158 template<class TypeTag>
159 struct EndTime<TypeTag, TTag::PowerInjectionBaseProblem>
160 {
161  using type = GetPropType<TypeTag, Scalar>;
162  static constexpr type value = 100;
163 };
164 
165 // The default for the initial time step size of the simulation
166 template<class TypeTag>
167 struct InitialTimeStepSize<TypeTag, TTag::PowerInjectionBaseProblem>
168 {
169  using type = GetPropType<TypeTag, Scalar>;
170  static constexpr type value = 1e-3;
171 };
172 
173 } // namespace Opm::Properties
174 
175 namespace Opm {
188 template <class TypeTag>
189 class PowerInjectionProblem : public GetPropType<TypeTag, Properties::BaseProblem>
190 {
191  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
192 
193  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
194  using GridView = GetPropType<TypeTag, Properties::GridView>;
195  using Indices = GetPropType<TypeTag, Properties::Indices>;
196  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
197  using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
198  using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
199  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
200  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
201  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
202  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
203  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
204 
205  enum {
206  // number of phases
207  numPhases = FluidSystem::numPhases,
208 
209  // phase indices
210  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
211  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
212 
213  // equation indices
214  contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
215 
216  // Grid and world dimension
217  dim = GridView::dimension,
218  dimWorld = GridView::dimensionworld
219  };
220 
221  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
222  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
223 
224  using CoordScalar = typename GridView::ctype;
225  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
226 
227  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
228 
229 public:
233  PowerInjectionProblem(Simulator& simulator)
234  : ParentType(simulator)
235  { }
236 
240  void finishInit()
241  {
242  ParentType::finishInit();
243 
244  eps_ = 3e-6;
245  FluidSystem::init();
246 
247  temperature_ = 273.15 + 26.6;
248 
249  // parameters for the Van Genuchten law
250  // alpha and n
251  materialParams_.setVgAlpha(0.00045);
252  materialParams_.setVgN(7.3);
253  materialParams_.finalize();
254 
255  K_ = this->toDimMatrix_(5.73e-08); // [m^2]
256 
257  setupInitialFluidState_();
258  }
259 
264 
268  std::string name() const
269  {
270  std::ostringstream oss;
271  oss << "powerinjection_";
272  if (std::is_same<GetPropType<TypeTag, Properties::FluxModule>,
273  Opm::DarcyFluxModule<TypeTag> >::value)
274  oss << "darcy";
275  else
276  oss << "forchheimer";
277 
278  if (std::is_same<GetPropType<TypeTag, Properties::LocalLinearizerSplice>,
279  Properties::TTag::AutoDiffLocalLinearizer>::value)
280  oss << "_" << "ad";
281  else
282  oss << "_" << "fd";
283 
284  return oss.str();
285  }
286 
290  void endTimeStep()
291  {
292 #ifndef NDEBUG
293  this->model().checkConservativeness();
294 
295  // Calculate storage terms
296  EqVector storage;
297  this->model().globalStorage(storage);
298 
299  // Write mass balance information for rank 0
300  if (this->gridView().comm().rank() == 0) {
301  std::cout << "Storage: " << storage << std::endl << std::flush;
302  }
303 #endif // NDEBUG
304  }
306 
311 
315  template <class Context>
316  const DimMatrix& intrinsicPermeability(const Context& /*context*/,
317  unsigned /*spaceIdx*/,
318  unsigned /*timeIdx*/) const
319  { return K_; }
320 
324  template <class Context>
325  Scalar ergunCoefficient(const Context& /*context*/,
326  unsigned /*spaceIdx*/,
327  unsigned /*timeIdx*/) const
328  { return 0.3866; }
329 
333  template <class Context>
334  Scalar porosity(const Context& /*context*/,
335  unsigned /*spaceIdx*/,
336  unsigned /*timeIdx*/) const
337  { return 0.558; }
338 
342  template <class Context>
343  const MaterialLawParams&
344  materialLawParams(const Context& /*context*/,
345  unsigned /*spaceIdx*/,
346  unsigned /*timeIdx*/) const
347  { return materialParams_; }
348 
352  template <class Context>
353  Scalar temperature(const Context& /*context*/,
354  unsigned /*spaceIdx*/,
355  unsigned /*timeIdx*/) const
356  { return temperature_; }
357 
359 
364 
371  template <class Context>
372  void boundary(BoundaryRateVector& values,
373  const Context& context,
374  unsigned spaceIdx,
375  unsigned timeIdx) const
376  {
377  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
378 
379  if (onLeftBoundary_(pos)) {
380  RateVector massRate(0.0);
381  massRate = 0.0;
382  massRate[contiNEqIdx] = -1.00; // kg / (m^2 * s)
383 
384  // impose a forced flow boundary
385  values.setMassRate(massRate);
386  }
387  else if (onRightBoundary_(pos))
388  // free flow boundary with initial condition on the right
389  values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
390  else
391  values.setNoFlow();
392  }
393 
395 
400 
404  template <class Context>
405  void initial(PrimaryVariables& values,
406  const Context& /*context*/,
407  unsigned /*spaceIdx*/,
408  unsigned /*timeIdx*/) const
409  {
410  // assign the primary variables
411  values.assignNaive(initialFluidState_);
412  }
413 
420  template <class Context>
421  void source(RateVector& rate,
422  const Context& /*context*/,
423  unsigned /*spaceIdx*/,
424  unsigned /*timeIdx*/) const
425  { rate = Scalar(0.0); }
426 
428 
429 private:
430  bool onLeftBoundary_(const GlobalPosition& pos) const
431  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
432 
433  bool onRightBoundary_(const GlobalPosition& pos) const
434  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
435 
436  void setupInitialFluidState_()
437  {
438  initialFluidState_.setTemperature(temperature_);
439 
440  Scalar Sw = 1.0;
441  initialFluidState_.setSaturation(wettingPhaseIdx, Sw);
442  initialFluidState_.setSaturation(nonWettingPhaseIdx, 1 - Sw);
443 
444  Scalar p = 1e5;
445  initialFluidState_.setPressure(wettingPhaseIdx, p);
446  initialFluidState_.setPressure(nonWettingPhaseIdx, p);
447 
448  typename FluidSystem::template ParameterCache<Scalar> paramCache;
449  paramCache.updateAll(initialFluidState_);
450  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
451  initialFluidState_.setDensity(phaseIdx,
452  FluidSystem::density(initialFluidState_, paramCache, phaseIdx));
453  initialFluidState_.setViscosity(phaseIdx,
454  FluidSystem::viscosity(initialFluidState_, paramCache, phaseIdx));
455  }
456  }
457 
458  DimMatrix K_;
459  MaterialLawParams materialParams_;
460 
461  Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
462  Scalar temperature_;
463  Scalar eps_;
464 };
465 
466 } // namespace Opm
467 
468 #endif
1D Problem with very fast injection of gas on the left.
Definition: powerinjectionproblem.hh:190
void finishInit()
Definition: powerinjectionproblem.hh:240
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:316
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:334
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:353
void endTimeStep()
Definition: powerinjectionproblem.hh:290
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:405
Scalar ergunCoefficient(const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:325
std::string name() const
Definition: powerinjectionproblem.hh:268
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:421
PowerInjectionProblem(Simulator &simulator)
Definition: powerinjectionproblem.hh:233
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: powerinjectionproblem.hh:372
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition: powerinjectionproblem.hh:344
Definition: powerinjectionproblem.hh:62