infiltrationproblem.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 */
27 #ifndef EWOMS_INFILTRATION_PROBLEM_HH
28 #define EWOMS_INFILTRATION_PROBLEM_HH
29 
31 
32 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
33 #include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>
34 #include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
35 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
36 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
37 #include <opm/material/heatconduction/Somerton.hpp>
38 #include <opm/common/Valgrind.hpp>
39 #include <opm/common/Unused.hpp>
40 
41 #include <dune/grid/yaspgrid.hh>
42 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
43 
44 #include <dune/common/version.hh>
45 #include <dune/common/fvector.hh>
46 #include <dune/common/fmatrix.hh>
47 
48 #include <sstream>
49 #include <string>
50 
51 namespace Ewoms {
52 template <class TypeTag>
54 }
55 
56 namespace Ewoms {
57 namespace Properties {
58 NEW_TYPE_TAG(InfiltrationBaseProblem);
59 
60 // Set the grid type
61 SET_TYPE_PROP(InfiltrationBaseProblem, Grid, Dune::YaspGrid<2>);
62 
63 // Set the problem property
64 SET_TYPE_PROP(InfiltrationBaseProblem, Problem,
66 
67 // Set the fluid system
69  InfiltrationBaseProblem, FluidSystem,
70  Opm::FluidSystems::H2OAirMesitylene<typename GET_PROP_TYPE(TypeTag, Scalar)>);
71 
72 // Enable gravity?
73 SET_BOOL_PROP(InfiltrationBaseProblem, EnableGravity, true);
74 
75 // Write newton convergence?
76 SET_BOOL_PROP(InfiltrationBaseProblem, NewtonWriteConvergence, false);
77 
78 // -1 backward differences, 0: central differences, +1: forward differences
79 SET_INT_PROP(InfiltrationBaseProblem, NumericDifferenceMethod, 1);
80 
81 // Set the material Law
82 SET_PROP(InfiltrationBaseProblem, MaterialLaw)
83 {
84 private:
85  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
86  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
87 
88  typedef Opm::ThreePhaseMaterialTraits<
89  Scalar,
90  /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
91  /*nonWettingPhaseIdx=*/FluidSystem::naplPhaseIdx,
92  /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> Traits;
93 
94 public:
95  typedef Opm::ThreePhaseParkerVanGenuchten<Traits> type;
96 };
97 
98 // Set the heat conduction law
99 SET_PROP(InfiltrationBaseProblem, HeatConductionLaw)
100 {
101 private:
102  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
103  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
104 
105 public:
106  // define the material law parameterized by absolute saturations
107  typedef Opm::Somerton<FluidSystem, Scalar> type;
108 };
109 
110 // The default for the end time of the simulation
111 SET_SCALAR_PROP(InfiltrationBaseProblem, EndTime, 6e3);
112 
113 // The default for the initial time step size of the simulation
114 SET_SCALAR_PROP(InfiltrationBaseProblem, InitialTimeStepSize, 60);
115 
116 // The default DGF file to load
117 SET_STRING_PROP(InfiltrationBaseProblem, GridFile,
118  "./data/infiltration_50x3.dgf");
119 } // namespace Properties
120 } // namespace Ewoms
121 
122 namespace Ewoms {
145 template <class TypeTag>
146 class InfiltrationProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
147 {
148  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
149 
150  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
151  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
152  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
153  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
154  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
155  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
156  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
157  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
158  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
159  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
160  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
161 
162  // copy some indices for convenience
163  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
164  enum {
165  // equation indices
166  conti0EqIdx = Indices::conti0EqIdx,
167 
168  // number of phases/components
169  numPhases = FluidSystem::numPhases,
170 
171  // component indices
172  NAPLIdx = FluidSystem::NAPLIdx,
173  H2OIdx = FluidSystem::H2OIdx,
174  airIdx = FluidSystem::airIdx,
175 
176  // phase indices
177  waterPhaseIdx = FluidSystem::waterPhaseIdx,
178  gasPhaseIdx = FluidSystem::gasPhaseIdx,
179  naplPhaseIdx = FluidSystem::naplPhaseIdx,
180 
181  // Grid and world dimension
182  dim = GridView::dimension,
183  dimWorld = GridView::dimensionworld
184  };
185 
186  typedef typename GridView::ctype CoordScalar;
187  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
188  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
189 
190 public:
195  : ParentType(simulator)
196  , eps_(1e-6)
197  { }
198 
202  void finishInit()
203  {
204  ParentType::finishInit();
205 
206  temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
207  FluidSystem::init(/*tempMin=*/temperature_ - 1,
208  /*tempMax=*/temperature_ + 1,
209  /*nTemp=*/3,
210  /*pressMin=*/0.8 * 1e5,
211  /*pressMax=*/3 * 1e5,
212  /*nPress=*/200);
213 
214  // intrinsic permeabilities
215  fineK_ = this->toDimMatrix_(1e-11);
216  coarseK_ = this->toDimMatrix_(1e-11);
217 
218  // porosities
219  porosity_ = 0.40;
220 
221  // residual saturations
222  materialParams_.setSwr(0.12);
223  materialParams_.setSwrx(0.12);
224  materialParams_.setSnr(0.07);
225  materialParams_.setSgr(0.03);
226 
227  // parameters for the three-phase van Genuchten law
228  materialParams_.setVgAlpha(0.0005);
229  materialParams_.setVgN(4.);
230  materialParams_.setkrRegardsSnr(false);
231 
232  materialParams_.finalize();
233  materialParams_.checkDefined();
234  }
235 
239 
247  { return true; }
248 
252  std::string name() const
253  {
254  std::ostringstream oss;
255  oss << "infiltration_" << Model::name();
256  return oss.str();
257  }
258 
262  void endTimeStep()
263  {
264 #ifndef NDEBUG
265  this->model().checkConservativeness();
266 
267  // Calculate storage terms
268  EqVector storage;
269  this->model().globalStorage(storage);
270 
271  // Write mass balance information for rank 0
272  if (this->gridView().comm().rank() == 0) {
273  std::cout << "Storage: " << storage << std::endl << std::flush;
274  }
275 #endif // NDEBUG
276  }
277 
281  template <class Context>
282  Scalar temperature(const Context& context OPM_UNUSED,
283  unsigned spaceIdx OPM_UNUSED,
284  unsigned timeIdx OPM_UNUSED) const
285  { return temperature_; }
286 
290  template <class Context>
291  const DimMatrix&
292  intrinsicPermeability(const Context& context,
293  unsigned spaceIdx,
294  unsigned timeIdx) const
295  {
296  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
297  if (isFineMaterial_(pos))
298  return fineK_;
299  return coarseK_;
300  }
301 
305  template <class Context>
306  Scalar porosity(const Context& context OPM_UNUSED,
307  unsigned spaceIdx OPM_UNUSED,
308  unsigned timeIdx OPM_UNUSED) const
309  { return porosity_; }
310 
314  template <class Context>
315  const MaterialLawParams&
316  materialLawParams(const Context& context OPM_UNUSED,
317  unsigned spaceIdx OPM_UNUSED,
318  unsigned timeIdx OPM_UNUSED) const
319  { return materialParams_; }
320 
326  template <class Context>
327  Scalar heatCapacitySolid(const Context& context OPM_UNUSED,
328  unsigned spaceIdx OPM_UNUSED,
329  unsigned timeIdx OPM_UNUSED) const
330  {
331  return 850. // specific heat capacity [J / (kg K)]
332  * 2650.; // density of sand [kg/m^3]
333  }
334 
336 
340 
345  template <class Context>
346  void boundary(BoundaryRateVector& values,
347  const Context& context,
348  unsigned spaceIdx,
349  unsigned timeIdx) const
350  {
351  const auto& pos = context.pos(spaceIdx, timeIdx);
352 
353  if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
354  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
355 
356  initialFluidState_(fs, context, spaceIdx, timeIdx);
357 
358  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
359  }
360  else if (onInlet_(pos)) {
361  RateVector molarRate(0.0);
362  molarRate[conti0EqIdx + NAPLIdx] = -0.001;
363 
364  values.setMolarRate(molarRate);
365  Opm::Valgrind::CheckDefined(values);
366  }
367  else
368  values.setNoFlow();
369  }
370 
372 
376 
381  template <class Context>
382  void initial(PrimaryVariables& values,
383  const Context& context,
384  unsigned spaceIdx,
385  unsigned timeIdx) const
386  {
387  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
388 
389  initialFluidState_(fs, context, spaceIdx, timeIdx);
390 
391  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
392  values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true);
393  Opm::Valgrind::CheckDefined(values);
394  }
395 
402  template <class Context>
403  void source(RateVector& rate,
404  const Context& context OPM_UNUSED,
405  unsigned spaceIdx OPM_UNUSED,
406  unsigned timeIdx OPM_UNUSED) const
407  { rate = Scalar(0.0); }
408 
410 
411 private:
412  bool onLeftBoundary_(const GlobalPosition& pos) const
413  { return pos[0] < eps_; }
414 
415  bool onRightBoundary_(const GlobalPosition& pos) const
416  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
417 
418  bool onLowerBoundary_(const GlobalPosition& pos) const
419  { return pos[1] < eps_; }
420 
421  bool onUpperBoundary_(const GlobalPosition& pos) const
422  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
423 
424  bool onInlet_(const GlobalPosition& pos) const
425  { return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
426 
427  template <class FluidState, class Context>
428  void initialFluidState_(FluidState& fs, const Context& context,
429  unsigned spaceIdx, unsigned timeIdx) const
430  {
431  const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
432  Scalar y = pos[1];
433  Scalar x = pos[0];
434 
435  Scalar densityW = 1000.0;
436  Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
437  if (pc < 0.0)
438  pc = 0.0;
439 
440  // set pressures
441  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
442  Scalar Sw = matParams.Swr();
443  Scalar Swr = matParams.Swr();
444  Scalar Sgr = matParams.Sgr();
445  if (Sw < Swr)
446  Sw = Swr;
447  if (Sw > 1 - Sgr)
448  Sw = 1 - Sgr;
449  Scalar Sg = 1 - Sw;
450 
451  Opm::Valgrind::CheckDefined(Sw);
452  Opm::Valgrind::CheckDefined(Sg);
453 
454  fs.setSaturation(waterPhaseIdx, Sw);
455  fs.setSaturation(gasPhaseIdx, Sg);
456  fs.setSaturation(naplPhaseIdx, 0);
457 
458  // set temperature of all phases
459  fs.setTemperature(temperature_);
460 
461  // compute pressures
462  Scalar pcAll[numPhases];
463  Scalar pg = 1e5;
464  if (onLeftBoundary_(pos))
465  pg += 10e3;
466  MaterialLaw::capillaryPressures(pcAll, matParams, fs);
467  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
468  fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
469 
470  // set composition of gas phase
471  fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
472  fs.setMoleFraction(gasPhaseIdx, airIdx,
473  1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
474  fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
475 
476  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
477  typename FluidSystem::template ParameterCache<Scalar> paramCache;
478  CFRP::solve(fs, paramCache, gasPhaseIdx,
479  /*setViscosity=*/false,
480  /*setEnthalpy=*/false);
481 
482  fs.setMoleFraction(waterPhaseIdx, H2OIdx,
483  1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
484  }
485 
486  bool isFineMaterial_(const GlobalPosition& pos) const
487  { return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
488 
489  DimMatrix fineK_;
490  DimMatrix coarseK_;
491 
492  Scalar porosity_;
493 
494  MaterialLawParams materialParams_;
495 
496  Scalar temperature_;
497  Scalar eps_;
498 };
499 } // namespace Ewoms
500 
501 #endif
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: infiltrationproblem.hh:403
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void endTimeStep()
Called by the simulator after each time integration.
Definition: infiltrationproblem.hh:262
std::string name() const
The problem name.
Definition: infiltrationproblem.hh:252
Definition: baseauxiliarymodule.hh:37
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:194
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:327
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: infiltrationproblem.hh:246
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:53
#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 finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: infiltrationproblem.hh:202
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: infiltrationproblem.hh:382
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:306
#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: infiltrationproblem.hh:282
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:292
Declares the properties required for the compositional multi-phase primary variable switching model...
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:316
#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
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: infiltrationproblem.hh:346
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394