My Project
obstacleproblem.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_OBSTACLE_PROBLEM_HH
29 #define EWOMS_OBSTACLE_PROBLEM_HH
30 
31 #include <opm/models/ncp/ncpproperties.hh>
32 
33 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
34 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
35 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
36 #include <opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <opm/material/thermal/ConstantSolidHeatCapLaw.hpp>
41 #include <opm/material/thermal/SomertonThermalConductionLaw.hpp>
42 #include <opm/material/common/Unused.hpp>
43 
44 #include <dune/grid/yaspgrid.hh>
45 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
46 
47 #include <dune/common/version.hh>
48 #include <dune/common/fvector.hh>
49 #include <dune/common/fmatrix.hh>
50 
51 #include <sstream>
52 #include <string>
53 #include <iostream>
54 
55 namespace Opm {
56 template <class TypeTag>
57 class ObstacleProblem;
58 }
59 
60 namespace Opm::Properties {
61 
62 namespace TTag {
64 }
65 
66 // Set the grid type
67 template<class TypeTag>
68 struct Grid<TypeTag, TTag::ObstacleBaseProblem> { using type = Dune::YaspGrid<2>; };
69 
70 // Set the problem property
71 template<class TypeTag>
72 struct Problem<TypeTag, TTag::ObstacleBaseProblem> { using type = Opm::ObstacleProblem<TypeTag>; };
73 
74 // Set fluid configuration
75 template<class TypeTag>
76 struct FluidSystem<TypeTag, TTag::ObstacleBaseProblem>
77 { using type = Opm::H2ON2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
78 
79 // Set the material Law
80 template<class TypeTag>
81 struct MaterialLaw<TypeTag, TTag::ObstacleBaseProblem>
82 {
83 private:
84  // define the material law
85  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
86  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
87  using MaterialTraits = Opm::TwoPhaseMaterialTraits<Scalar,
88  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
89  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>;
90 
91  using EffMaterialLaw = Opm::LinearMaterial<MaterialTraits>;
92 
93 public:
94  using type = Opm::EffToAbsLaw<EffMaterialLaw>;
95 };
96 
97 // Set the thermal conduction law
98 template<class TypeTag>
99 struct ThermalConductionLaw<TypeTag, TTag::ObstacleBaseProblem>
100 {
101 private:
102  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
103  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
104 
105 public:
106  // define the material law parameterized by absolute saturations
107  using type = Opm::SomertonThermalConductionLaw<FluidSystem, Scalar>;
108 };
109 
110 // set the energy storage law for the solid phase
111 template<class TypeTag>
112 struct SolidEnergyLaw<TypeTag, TTag::ObstacleBaseProblem>
113 { using type = Opm::ConstantSolidHeatCapLaw<GetPropType<TypeTag, Properties::Scalar>>; };
114 
115 // Enable gravity
116 template<class TypeTag>
117 struct EnableGravity<TypeTag, TTag::ObstacleBaseProblem> { static constexpr bool value = true; };
118 
119 // The default for the end time of the simulation
120 template<class TypeTag>
121 struct EndTime<TypeTag, TTag::ObstacleBaseProblem>
122 {
123  using type = GetPropType<TypeTag, Scalar>;
124  static constexpr type value = 1e4;
125 };
126 
127 // The default for the initial time step size of the simulation
128 template<class TypeTag>
129 struct InitialTimeStepSize<TypeTag, TTag::ObstacleBaseProblem>
130 {
131  using type = GetPropType<TypeTag, Scalar>;
132  static constexpr type value = 250;
133 };
134 
135 // The default DGF file to load
136 template<class TypeTag>
137 struct GridFile<TypeTag, TTag::ObstacleBaseProblem> { static constexpr auto value = "./data/obstacle_24x16.dgf"; };
138 
139 } // namespace Opm::Properties
140 
141 namespace Opm {
168 template <class TypeTag>
169 class ObstacleProblem : public GetPropType<TypeTag, Properties::BaseProblem>
170 {
171  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
172 
173  using GridView = GetPropType<TypeTag, Properties::GridView>;
174  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
175  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
176  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
177  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
178  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
179  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
180  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
181  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
182  using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
183  using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
184 
185  enum {
186  // Grid and world dimension
187  dim = GridView::dimension,
188  dimWorld = GridView::dimensionworld,
189  numPhases = getPropValue<TypeTag, Properties::NumPhases>(),
190  gasPhaseIdx = FluidSystem::gasPhaseIdx,
191  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
192  H2OIdx = FluidSystem::H2OIdx,
193  N2Idx = FluidSystem::N2Idx
194  };
195 
196  using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
197  using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
198  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
199  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
200  using Model = GetPropType<TypeTag, Properties::Model>;
201 
202 public:
206  ObstacleProblem(Simulator& simulator)
207  : ParentType(simulator)
208  { }
209 
213  void finishInit()
214  {
215  ParentType::finishInit();
216 
217  eps_ = 1e-6;
218  temperature_ = 273.15 + 25; // -> 25°C
219 
220  // initialize the tables of the fluid system
221  Scalar Tmin = temperature_ - 1.0;
222  Scalar Tmax = temperature_ + 1.0;
223  unsigned nT = 3;
224 
225  Scalar pmin = 1.0e5 * 0.75;
226  Scalar pmax = 2.0e5 * 1.25;
227  unsigned np = 1000;
228 
229  FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
230 
231  // intrinsic permeabilities
232  coarseK_ = this->toDimMatrix_(1e-12);
233  fineK_ = this->toDimMatrix_(1e-15);
234 
235  // the porosity
236  finePorosity_ = 0.3;
237  coarsePorosity_ = 0.3;
238 
239  // residual saturations
240  fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
241  fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
242  coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
243  coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
244 
245  // parameters for the linear law, i.e. minimum and maximum
246  // pressures
247  fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
248  fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
249  coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
250  coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
251 
252  /*
253  // entry pressures for Brooks-Corey
254  fineMaterialParams_.setEntryPressure(5e3);
255  coarseMaterialParams_.setEntryPressure(1e3);
256 
257  // Brooks-Corey shape parameters
258  fineMaterialParams_.setLambda(2);
259  coarseMaterialParams_.setLambda(2);
260  */
261 
262  fineMaterialParams_.finalize();
263  coarseMaterialParams_.finalize();
264 
265  // parameters for the somerton law of thermal conduction
266  computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
267  computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
268 
269  // assume constant volumetric heat capacity and granite
270  solidEnergyLawParams_.setSolidHeatCapacity(790.0 // specific heat capacity of granite [J / (kg K)]
271  * 2700.0); // density of granite [kg/m^3]
272  solidEnergyLawParams_.finalize();
273 
274  initFluidStates_();
275  }
276 
280  void endTimeStep()
281  {
282 #ifndef NDEBUG
283  this->model().checkConservativeness();
284 
285  // Calculate storage terms of the individual phases
286  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
287  PrimaryVariables phaseStorage;
288  this->model().globalPhaseStorage(phaseStorage, phaseIdx);
289 
290  if (this->gridView().comm().rank() == 0) {
291  std::cout << "Storage in " << FluidSystem::phaseName(phaseIdx)
292  << "Phase: [" << phaseStorage << "]"
293  << "\n" << std::flush;
294  }
295  }
296 
297  // Calculate total storage terms
298  EqVector storage;
299  this->model().globalStorage(storage);
300 
301  // Write mass balance information for rank 0
302  if (this->gridView().comm().rank() == 0) {
303  std::cout << "Storage total: [" << storage << "]"
304  << "\n" << std::flush;
305  }
306 #endif // NDEBUG
307  }
308 
313 
317  std::string name() const
318  {
319  std::ostringstream oss;
320  oss << "obstacle"
321  << "_" << Model::name();
322  return oss.str();
323  }
324 
330  template <class Context>
331  Scalar temperature(const Context& context OPM_UNUSED,
332  unsigned spaceIdx OPM_UNUSED,
333  unsigned timeIdx OPM_UNUSED) const
334  { return temperature_; }
335 
339  template <class Context>
340  const DimMatrix&
341  intrinsicPermeability(const Context& context,
342  unsigned spaceIdx,
343  unsigned timeIdx) const
344  {
345  if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
346  return fineK_;
347  return coarseK_;
348  }
349 
353  template <class Context>
354  Scalar porosity(const Context& context,
355  unsigned spaceIdx,
356  unsigned timeIdx) const
357  {
358  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
359  if (isFineMaterial_(pos))
360  return finePorosity_;
361  else
362  return coarsePorosity_;
363  }
364 
368  template <class Context>
369  const MaterialLawParams&
370  materialLawParams(const Context& context,
371  unsigned spaceIdx,
372  unsigned timeIdx) const
373  {
374  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
375  if (isFineMaterial_(pos))
376  return fineMaterialParams_;
377  else
378  return coarseMaterialParams_;
379  }
380 
386  template <class Context>
387  const SolidEnergyLawParams&
388  solidEnergyLawParams(const Context& context OPM_UNUSED,
389  unsigned spaceIdx OPM_UNUSED,
390  unsigned timeIdx OPM_UNUSED) const
391  { return solidEnergyLawParams_; }
392 
396  template <class Context>
397  const ThermalConductionLawParams &
398  thermalConductionParams(const Context& context,
399  unsigned spaceIdx,
400  unsigned timeIdx) const
401  {
402  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
403  if (isFineMaterial_(pos))
404  return fineThermalCondParams_;
405  return coarseThermalCondParams_;
406  }
407 
409 
414 
418  template <class Context>
419  void boundary(BoundaryRateVector& values,
420  const Context& context,
421  unsigned spaceIdx,
422  unsigned timeIdx) const
423  {
424  const auto& pos = context.pos(spaceIdx, timeIdx);
425 
426  if (onInlet_(pos))
427  values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
428  else if (onOutlet_(pos))
429  values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
430  else
431  values.setNoFlow();
432  }
433 
435 
440 
444  template <class Context>
445  void initial(PrimaryVariables& values,
446  const Context& context,
447  unsigned spaceIdx,
448  unsigned timeIdx) const
449  {
450  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
451  values.assignMassConservative(outletFluidState_, matParams);
452  }
453 
460  template <class Context>
461  void source(RateVector& rate,
462  const Context& context OPM_UNUSED,
463  unsigned spaceIdx OPM_UNUSED,
464  unsigned timeIdx OPM_UNUSED) const
465  { rate = 0.0; }
466 
468 
469 private:
474  bool isFineMaterial_(const GlobalPosition& pos) const
475  { return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
476 
477  bool onInlet_(const GlobalPosition& globalPos) const
478  {
479  Scalar x = globalPos[0];
480  Scalar y = globalPos[1];
481  return x >= 60 - eps_ && y <= 10;
482  }
483 
484  bool onOutlet_(const GlobalPosition& globalPos) const
485  {
486  Scalar x = globalPos[0];
487  Scalar y = globalPos[1];
488  return x < eps_ && y <= 10;
489  }
490 
491  void initFluidStates_()
492  {
493  initFluidState_(inletFluidState_, coarseMaterialParams_,
494  /*isInlet=*/true);
495  initFluidState_(outletFluidState_, coarseMaterialParams_,
496  /*isInlet=*/false);
497  }
498 
499  template <class FluidState>
500  void initFluidState_(FluidState& fs, const MaterialLawParams& matParams, bool isInlet)
501  {
502  unsigned refPhaseIdx;
503  unsigned otherPhaseIdx;
504 
505  // set the fluid temperatures
506  fs.setTemperature(temperature_);
507 
508  if (isInlet) {
509  // only liquid on inlet
510  refPhaseIdx = liquidPhaseIdx;
511  otherPhaseIdx = gasPhaseIdx;
512 
513  // set liquid saturation
514  fs.setSaturation(liquidPhaseIdx, 1.0);
515 
516  // set pressure of the liquid phase
517  fs.setPressure(liquidPhaseIdx, 2e5);
518 
519  // set the liquid composition to pure water
520  fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
521  fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
522  }
523  else {
524  // elsewhere, only gas
525  refPhaseIdx = gasPhaseIdx;
526  otherPhaseIdx = liquidPhaseIdx;
527 
528  // set gas saturation
529  fs.setSaturation(gasPhaseIdx, 1.0);
530 
531  // set pressure of the gas phase
532  fs.setPressure(gasPhaseIdx, 1e5);
533 
534  // set the gas composition to 99% nitrogen and 1% steam
535  fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
536  fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
537  }
538 
539  // set the other saturation
540  fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
541 
542  // calulate the capillary pressure
543  PhaseVector pC;
544  MaterialLaw::capillaryPressures(pC, matParams, fs);
545  fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
546  + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
547 
548  // make the fluid state consistent with local thermodynamic
549  // equilibrium
550  using ComputeFromReferencePhase = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
551 
552  typename FluidSystem::template ParameterCache<Scalar> paramCache;
553  ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
554  /*setViscosity=*/true,
555  /*setEnthalpy=*/false);
556  }
557 
558  void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
559  {
560  Scalar lambdaWater = 0.6;
561  Scalar lambdaGranite = 2.8;
562 
563  Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
564  * std::pow(lambdaWater, poro);
565  Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
566 
567  params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
568  params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
569  params.setVacuumLambda(lambdaDry);
570  }
571 
572  DimMatrix coarseK_;
573  DimMatrix fineK_;
574 
575  Scalar coarsePorosity_;
576  Scalar finePorosity_;
577 
578  MaterialLawParams fineMaterialParams_;
579  MaterialLawParams coarseMaterialParams_;
580 
581  ThermalConductionLawParams fineThermalCondParams_;
582  ThermalConductionLawParams coarseThermalCondParams_;
583  SolidEnergyLawParams solidEnergyLawParams_;
584 
585  Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
586  Opm::CompositionalFluidState<Scalar, FluidSystem> outletFluidState_;
587 
588  Scalar temperature_;
589  Scalar eps_;
590 };
591 } // namespace Opm
592 
593 #endif
Problem where liquid water is first stopped by a low-permeability lens and then seeps though it.
Definition: obstacleproblem.hh:170
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:354
void finishInit()
Definition: obstacleproblem.hh:213
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:370
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:461
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:331
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:341
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:419
std::string name() const
Definition: obstacleproblem.hh:317
void endTimeStep()
Definition: obstacleproblem.hh:280
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:445
ObstacleProblem(Simulator &simulator)
Definition: obstacleproblem.hh:206
const ThermalConductionLawParams & thermalConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:398
const SolidEnergyLawParams & solidEnergyLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Return the parameters for the energy storage law of the rock.
Definition: obstacleproblem.hh:388
Definition: obstacleproblem.hh:63