All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 
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/heatconduction/Somerton.hpp>
41 #include <opm/common/Unused.hpp>
42 
43 #include <dune/grid/yaspgrid.hh>
44 #include <dune/grid/io/file/dgfparser/dgfyasp.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 <iostream>
53 
54 namespace Ewoms {
55 template <class TypeTag>
57 }
58 
59 namespace Ewoms {
60 namespace Properties {
61 NEW_TYPE_TAG(ObstacleBaseProblem);
62 
63 // Set the grid type
64 SET_TYPE_PROP(ObstacleBaseProblem, Grid, Dune::YaspGrid<2>);
65 
66 // Set the problem property
67 SET_TYPE_PROP(ObstacleBaseProblem, Problem, Ewoms::ObstacleProblem<TypeTag>);
68 
69 // Set fluid configuration
70 SET_TYPE_PROP(ObstacleBaseProblem, FluidSystem,
71  Opm::FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar)>);
72 
73 // Set the material Law
74 SET_PROP(ObstacleBaseProblem, MaterialLaw)
75 {
76 private:
77  // define the material law
78  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
79  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
80  typedef Opm::TwoPhaseMaterialTraits<Scalar,
81  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
82  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>
83  MaterialTraits;
84 
85  typedef Opm::LinearMaterial<MaterialTraits> EffMaterialLaw;
86 
87 public:
88  typedef Opm::EffToAbsLaw<EffMaterialLaw> type;
89 };
90 
91 // Set the heat conduction law
92 SET_PROP(ObstacleBaseProblem, HeatConductionLaw)
93 {
94 private:
95  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
96  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
97 
98 public:
99  // define the material law parameterized by absolute saturations
100  typedef Opm::Somerton<FluidSystem, Scalar> type;
101 };
102 
103 // Enable gravity
104 SET_BOOL_PROP(ObstacleBaseProblem, EnableGravity, true);
105 
106 // The default for the end time of the simulation
107 SET_SCALAR_PROP(ObstacleBaseProblem, EndTime, 1e4);
108 
109 // The default for the initial time step size of the simulation
110 SET_SCALAR_PROP(ObstacleBaseProblem, InitialTimeStepSize, 250);
111 
112 // The default DGF file to load
113 SET_STRING_PROP(ObstacleBaseProblem, GridFile, "./data/obstacle_24x16.dgf");
114 } // namespace Properties
115 } // namespace Ewoms
116 
117 namespace Ewoms {
144 template <class TypeTag>
145 class ObstacleProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
146 {
147  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
148 
149  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
150  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
151  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
152  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
153  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
154  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
155  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
156  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
157  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
158  typedef typename GET_PROP_TYPE(TypeTag, HeatConductionLaw) HeatConductionLaw;
159  typedef typename HeatConductionLaw::Params HeatConductionLawParams;
160 
161  enum {
162  // Grid and world dimension
163  dim = GridView::dimension,
164  dimWorld = GridView::dimensionworld,
165  numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
166  gasPhaseIdx = FluidSystem::gasPhaseIdx,
167  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
168  H2OIdx = FluidSystem::H2OIdx,
169  N2Idx = FluidSystem::N2Idx
170  };
171 
172  typedef Dune::FieldVector<typename GridView::ctype, dimWorld> GlobalPosition;
173  typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
174  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
175  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
176  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
177 
178 public:
182  ObstacleProblem(Simulator& simulator)
183  : ParentType(simulator)
184  { }
185 
189  void finishInit()
190  {
191  ParentType::finishInit();
192 
193  eps_ = 1e-6;
194  temperature_ = 273.15 + 25; // -> 25°C
195 
196  // initialize the tables of the fluid system
197  Scalar Tmin = temperature_ - 1.0;
198  Scalar Tmax = temperature_ + 1.0;
199  unsigned nT = 3;
200 
201  Scalar pmin = 1.0e5 * 0.75;
202  Scalar pmax = 2.0e5 * 1.25;
203  unsigned np = 1000;
204 
205  FluidSystem::init(Tmin, Tmax, nT, pmin, pmax, np);
206 
207  // intrinsic permeabilities
208  coarseK_ = this->toDimMatrix_(1e-12);
209  fineK_ = this->toDimMatrix_(1e-15);
210 
211  // the porosity
212  finePorosity_ = 0.3;
213  coarsePorosity_ = 0.3;
214 
215  // residual saturations
216  fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
217  fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
218  coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.0);
219  coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
220 
221  // parameters for the linear law, i.e. minimum and maximum
222  // pressures
223  fineMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
224  fineMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
225  coarseMaterialParams_.setPcMinSat(liquidPhaseIdx, 0.0);
226  coarseMaterialParams_.setPcMaxSat(liquidPhaseIdx, 0.0);
227 
228  /*
229  // entry pressures for Brooks-Corey
230  fineMaterialParams_.setEntryPressure(5e3);
231  coarseMaterialParams_.setEntryPressure(1e3);
232 
233  // Brooks-Corey shape parameters
234  fineMaterialParams_.setLambda(2);
235  coarseMaterialParams_.setLambda(2);
236  */
237 
238  fineMaterialParams_.finalize();
239  coarseMaterialParams_.finalize();
240 
241  // parameters for the somerton law of heat conduction
242  computeHeatCondParams_(fineHeatCondParams_, finePorosity_);
243  computeHeatCondParams_(coarseHeatCondParams_, coarsePorosity_);
244 
245  initFluidStates_();
246  }
247 
251  void endTimeStep()
252  {
253 #ifndef NDEBUG
254  this->model().checkConservativeness();
255 
256  // Calculate storage terms of the individual phases
257  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
258  PrimaryVariables phaseStorage;
259  this->model().globalPhaseStorage(phaseStorage, phaseIdx);
260 
261  if (this->gridView().comm().rank() == 0) {
262  std::cout << "Storage in " << FluidSystem::phaseName(phaseIdx)
263  << "Phase: [" << phaseStorage << "]"
264  << "\n" << std::flush;
265  }
266  }
267 
268  // Calculate total storage terms
269  EqVector storage;
270  this->model().globalStorage(storage);
271 
272  // Write mass balance information for rank 0
273  if (this->gridView().comm().rank() == 0) {
274  std::cout << "Storage total: [" << storage << "]"
275  << "\n" << std::flush;
276  }
277 #endif // NDEBUG
278  }
279 
283 
288  std::string name() const
289  {
290  std::ostringstream oss;
291  oss << "obstacle"
292  << "_" << Model::name();
293  return oss.str();
294  }
295 
301  template <class Context>
302  Scalar temperature(const Context& context OPM_UNUSED,
303  unsigned spaceIdx OPM_UNUSED,
304  unsigned timeIdx OPM_UNUSED) const
305  { return temperature_; }
306 
310  template <class Context>
311  const DimMatrix&
312  intrinsicPermeability(const Context& context,
313  unsigned spaceIdx,
314  unsigned timeIdx) const
315  {
316  if (isFineMaterial_(context.pos(spaceIdx, timeIdx)))
317  return fineK_;
318  return coarseK_;
319  }
320 
324  template <class Context>
325  Scalar porosity(const Context& context,
326  unsigned spaceIdx,
327  unsigned timeIdx) const
328  {
329  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
330  if (isFineMaterial_(pos))
331  return finePorosity_;
332  else
333  return coarsePorosity_;
334  }
335 
339  template <class Context>
340  const MaterialLawParams&
341  materialLawParams(const Context& context,
342  unsigned spaceIdx,
343  unsigned timeIdx) const
344  {
345  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
346  if (isFineMaterial_(pos))
347  return fineMaterialParams_;
348  else
349  return coarseMaterialParams_;
350  }
351 
358  template <class Context>
359  Scalar heatCapacitySolid(const Context& context OPM_UNUSED,
360  unsigned spaceIdx OPM_UNUSED,
361  unsigned timeIdx OPM_UNUSED) const
362  {
363  return 790 // specific heat capacity of granite [J / (kg K)]
364  * 2700; // density of granite [kg/m^3]
365  }
366 
370  template <class Context>
371  const HeatConductionLawParams &
372  heatConductionParams(const Context& context,
373  unsigned spaceIdx,
374  unsigned timeIdx) const
375  {
376  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
377  if (isFineMaterial_(pos))
378  return fineHeatCondParams_;
379  return coarseHeatCondParams_;
380  }
381 
383 
387 
392  template <class Context>
393  void boundary(BoundaryRateVector& values,
394  const Context& context,
395  unsigned spaceIdx,
396  unsigned timeIdx) const
397  {
398  const auto& pos = context.pos(spaceIdx, timeIdx);
399 
400  if (onInlet_(pos))
401  values.setFreeFlow(context, spaceIdx, timeIdx, inletFluidState_);
402  else if (onOutlet_(pos))
403  values.setFreeFlow(context, spaceIdx, timeIdx, outletFluidState_);
404  else
405  values.setNoFlow();
406  }
407 
409 
413 
418  template <class Context>
419  void initial(PrimaryVariables& values,
420  const Context& context,
421  unsigned spaceIdx,
422  unsigned timeIdx) const
423  {
424  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
425  values.assignMassConservative(outletFluidState_, matParams);
426  }
427 
434  template <class Context>
435  void source(RateVector& rate,
436  const Context& context OPM_UNUSED,
437  unsigned spaceIdx OPM_UNUSED,
438  unsigned timeIdx OPM_UNUSED) const
439  { rate = 0.0; }
440 
442 
443 private:
448  bool isFineMaterial_(const GlobalPosition& pos) const
449  { return 10 <= pos[0] && pos[0] <= 20 && 0 <= pos[1] && pos[1] <= 35; }
450 
451  bool onInlet_(const GlobalPosition& globalPos) const
452  {
453  Scalar x = globalPos[0];
454  Scalar y = globalPos[1];
455  return x >= 60 - eps_ && y <= 10;
456  }
457 
458  bool onOutlet_(const GlobalPosition& globalPos) const
459  {
460  Scalar x = globalPos[0];
461  Scalar y = globalPos[1];
462  return x < eps_ && y <= 10;
463  }
464 
465  void initFluidStates_()
466  {
467  initFluidState_(inletFluidState_, coarseMaterialParams_,
468  /*isInlet=*/true);
469  initFluidState_(outletFluidState_, coarseMaterialParams_,
470  /*isInlet=*/false);
471  }
472 
473  template <class FluidState>
474  void initFluidState_(FluidState& fs, const MaterialLawParams& matParams, bool isInlet)
475  {
476  unsigned refPhaseIdx;
477  unsigned otherPhaseIdx;
478 
479  // set the fluid temperatures
480  fs.setTemperature(temperature_);
481 
482  if (isInlet) {
483  // only liquid on inlet
484  refPhaseIdx = liquidPhaseIdx;
485  otherPhaseIdx = gasPhaseIdx;
486 
487  // set liquid saturation
488  fs.setSaturation(liquidPhaseIdx, 1.0);
489 
490  // set pressure of the liquid phase
491  fs.setPressure(liquidPhaseIdx, 2e5);
492 
493  // set the liquid composition to pure water
494  fs.setMoleFraction(liquidPhaseIdx, N2Idx, 0.0);
495  fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
496  }
497  else {
498  // elsewhere, only gas
499  refPhaseIdx = gasPhaseIdx;
500  otherPhaseIdx = liquidPhaseIdx;
501 
502  // set gas saturation
503  fs.setSaturation(gasPhaseIdx, 1.0);
504 
505  // set pressure of the gas phase
506  fs.setPressure(gasPhaseIdx, 1e5);
507 
508  // set the gas composition to 99% nitrogen and 1% steam
509  fs.setMoleFraction(gasPhaseIdx, N2Idx, 0.99);
510  fs.setMoleFraction(gasPhaseIdx, H2OIdx, 0.01);
511  }
512 
513  // set the other saturation
514  fs.setSaturation(otherPhaseIdx, 1.0 - fs.saturation(refPhaseIdx));
515 
516  // calulate the capillary pressure
517  PhaseVector pC;
518  MaterialLaw::capillaryPressures(pC, matParams, fs);
519  fs.setPressure(otherPhaseIdx, fs.pressure(refPhaseIdx)
520  + (pC[otherPhaseIdx] - pC[refPhaseIdx]));
521 
522  // make the fluid state consistent with local thermodynamic
523  // equilibrium
524  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem>
525  ComputeFromReferencePhase;
526 
527  typename FluidSystem::template ParameterCache<Scalar> paramCache;
528  ComputeFromReferencePhase::solve(fs, paramCache, refPhaseIdx,
529  /*setViscosity=*/false,
530  /*setEnthalpy=*/false);
531  }
532 
533  void computeHeatCondParams_(HeatConductionLawParams& params, Scalar poro)
534  {
535  Scalar lambdaWater = 0.6;
536  Scalar lambdaGranite = 2.8;
537 
538  Scalar lambdaWet = std::pow(lambdaGranite, (1 - poro))
539  * std::pow(lambdaWater, poro);
540  Scalar lambdaDry = std::pow(lambdaGranite, (1 - poro));
541 
542  params.setFullySaturatedLambda(gasPhaseIdx, lambdaDry);
543  params.setFullySaturatedLambda(liquidPhaseIdx, lambdaWet);
544  params.setVacuumLambda(lambdaDry);
545  }
546 
547  DimMatrix coarseK_;
548  DimMatrix fineK_;
549 
550  Scalar coarsePorosity_;
551  Scalar finePorosity_;
552 
553  MaterialLawParams fineMaterialParams_;
554  MaterialLawParams coarseMaterialParams_;
555 
556  HeatConductionLawParams fineHeatCondParams_;
557  HeatConductionLawParams coarseHeatCondParams_;
558 
559  Opm::CompositionalFluidState<Scalar, FluidSystem> inletFluidState_;
560  Opm::CompositionalFluidState<Scalar, FluidSystem> outletFluidState_;
561 
562  Scalar temperature_;
563  Scalar eps_;
564 };
565 } // namespace Ewoms
566 
567 #endif
std::string name() const
The problem name.
Definition: obstacleproblem.hh:288
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:359
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:312
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: obstacleproblem.hh:393
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: obstacleproblem.hh:189
#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.
void endTimeStep()
Called by the simulator after each time integration.
Definition: obstacleproblem.hh:251
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:325
const HeatConductionLawParams & heatConductionParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:372
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: obstacleproblem.hh:419
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: obstacleproblem.hh:302
Problem where liquid water is first stopped by a low-permeability lens and then seeps though it...
Definition: obstacleproblem.hh:56
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: obstacleproblem.hh:341
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: obstacleproblem.hh:435
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394