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