My Project
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 
30 #include <opm/models/pvs/pvsproperties.hh>
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/common/Valgrind.hpp>
38 
39 #include <dune/grid/yaspgrid.hh>
40 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
41 
42 #include <dune/common/version.hh>
43 #include <dune/common/fvector.hh>
44 #include <dune/common/fmatrix.hh>
45 
46 #include <sstream>
47 #include <string>
48 
49 namespace Opm {
50 template <class TypeTag>
51 class InfiltrationProblem;
52 }
53 
54 namespace Opm::Properties {
55 
56 namespace TTag {
58 }
59 
60 // Set the grid type
61 template<class TypeTag>
62 struct Grid<TypeTag, TTag::InfiltrationBaseProblem> { using type = Dune::YaspGrid<2>; };
63 
64 // Set the problem property
65 template<class TypeTag>
66 struct Problem<TypeTag, TTag::InfiltrationBaseProblem> { using type = Opm::InfiltrationProblem<TypeTag>; };
67 
68 // Set the fluid system
69 template<class TypeTag>
70 struct FluidSystem<TypeTag, TTag::InfiltrationBaseProblem>
71 { using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
72 
73 // Enable gravity?
74 template<class TypeTag>
75 struct EnableGravity<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr bool value = true; };
76 
77 // Write newton convergence?
78 template<class TypeTag>
79 struct NewtonWriteConvergence<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr bool value = false; };
80 
81 // -1 backward differences, 0: central differences, +1: forward differences
82 template<class TypeTag>
83 struct NumericDifferenceMethod<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr int value = 1; };
84 
85 // Set the material Law
86 template<class TypeTag>
87 struct MaterialLaw<TypeTag, TTag::InfiltrationBaseProblem>
88 {
89 private:
90  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
92 
93  using Traits= Opm::ThreePhaseMaterialTraits<
94  Scalar,
95  /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
96  /*nonWettingPhaseIdx=*/FluidSystem::naplPhaseIdx,
97  /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
98 
99 public:
100  using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
101 };
102 
103 // The default for the end time of the simulation
104 template<class TypeTag>
105 struct EndTime<TypeTag, TTag::InfiltrationBaseProblem>
106 {
107  using type = GetPropType<TypeTag, Scalar>;
108  static constexpr type value = 6e3;
109 };
110 
111 // The default for the initial time step size of the simulation
112 template<class TypeTag>
113 struct InitialTimeStepSize<TypeTag, TTag::InfiltrationBaseProblem>
114 {
115  using type = GetPropType<TypeTag, Scalar>;
116  static constexpr type value = 60;
117 };
118 
119 // The default DGF file to load
120 template<class TypeTag>
121 struct GridFile<TypeTag, TTag::InfiltrationBaseProblem>
122 { static constexpr auto value = "./data/infiltration_50x3.dgf"; };
123 
124 } // namespace Opm::Properties
125 
126 namespace Opm {
149 template <class TypeTag>
150 class InfiltrationProblem : public GetPropType<TypeTag, Properties::BaseProblem>
151 {
152  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
153 
154  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
155  using GridView = GetPropType<TypeTag, Properties::GridView>;
156  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
157  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
158  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
159  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
160  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
161  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
162  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
163  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
164  using Model = GetPropType<TypeTag, Properties::Model>;
165 
166  // copy some indices for convenience
167  using Indices = GetPropType<TypeTag, Properties::Indices>;
168  enum {
169  // equation indices
170  conti0EqIdx = Indices::conti0EqIdx,
171 
172  // number of phases/components
173  numPhases = FluidSystem::numPhases,
174 
175  // component indices
176  NAPLIdx = FluidSystem::NAPLIdx,
177  H2OIdx = FluidSystem::H2OIdx,
178  airIdx = FluidSystem::airIdx,
179 
180  // phase indices
181  waterPhaseIdx = FluidSystem::waterPhaseIdx,
182  gasPhaseIdx = FluidSystem::gasPhaseIdx,
183  naplPhaseIdx = FluidSystem::naplPhaseIdx,
184 
185  // Grid and world dimension
186  dim = GridView::dimension,
187  dimWorld = GridView::dimensionworld
188  };
189 
190  using CoordScalar = typename GridView::ctype;
191  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
192  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
193 
194 public:
198  InfiltrationProblem(Simulator& simulator)
199  : ParentType(simulator)
200  , eps_(1e-6)
201  { }
202 
206  void finishInit()
207  {
208  ParentType::finishInit();
209 
210  temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
211  FluidSystem::init(/*tempMin=*/temperature_ - 1,
212  /*tempMax=*/temperature_ + 1,
213  /*nTemp=*/3,
214  /*pressMin=*/0.8 * 1e5,
215  /*pressMax=*/3 * 1e5,
216  /*nPress=*/200);
217 
218  // intrinsic permeabilities
219  fineK_ = this->toDimMatrix_(1e-11);
220  coarseK_ = this->toDimMatrix_(1e-11);
221 
222  // porosities
223  porosity_ = 0.40;
224 
225  // residual saturations
226  materialParams_.setSwr(0.12);
227  materialParams_.setSwrx(0.12);
228  materialParams_.setSnr(0.07);
229  materialParams_.setSgr(0.03);
230 
231  // parameters for the three-phase van Genuchten law
232  materialParams_.setVgAlpha(0.0005);
233  materialParams_.setVgN(4.);
234  materialParams_.setkrRegardsSnr(false);
235 
236  materialParams_.finalize();
237  materialParams_.checkDefined();
238  }
239 
244 
251  { return true; }
252 
256  std::string name() const
257  {
258  std::ostringstream oss;
259  oss << "infiltration_" << Model::name();
260  return oss.str();
261  }
262 
266  void endTimeStep()
267  {
268 #ifndef NDEBUG
269  this->model().checkConservativeness();
270 
271  // Calculate storage terms
272  EqVector storage;
273  this->model().globalStorage(storage);
274 
275  // Write mass balance information for rank 0
276  if (this->gridView().comm().rank() == 0) {
277  std::cout << "Storage: " << storage << std::endl << std::flush;
278  }
279 #endif // NDEBUG
280  }
281 
285  template <class Context>
286  Scalar temperature(const Context& /*context*/,
287  unsigned /*spaceIdx*/,
288  unsigned /*timeIdx*/) const
289  { return temperature_; }
290 
294  template <class Context>
295  const DimMatrix&
296  intrinsicPermeability(const Context& context,
297  unsigned spaceIdx,
298  unsigned timeIdx) const
299  {
300  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
301  if (isFineMaterial_(pos))
302  return fineK_;
303  return coarseK_;
304  }
305 
309  template <class Context>
310  Scalar porosity(const Context& /*context*/,
311  unsigned /*spaceIdx*/,
312  unsigned /*timeIdx*/) const
313  { return porosity_; }
314 
318  template <class Context>
319  const MaterialLawParams&
320  materialLawParams(const Context& /*context*/,
321  unsigned /*spaceIdx*/,
322  unsigned /*timeIdx*/) const
323  { return materialParams_; }
324 
326 
331 
335  template <class Context>
336  void boundary(BoundaryRateVector& values,
337  const Context& context,
338  unsigned spaceIdx,
339  unsigned timeIdx) const
340  {
341  const auto& pos = context.pos(spaceIdx, timeIdx);
342 
343  if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
344  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
345 
346  initialFluidState_(fs, context, spaceIdx, timeIdx);
347 
348  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
349  }
350  else if (onInlet_(pos)) {
351  RateVector molarRate(0.0);
352  molarRate[conti0EqIdx + NAPLIdx] = -0.001;
353 
354  values.setMolarRate(molarRate);
355  Opm::Valgrind::CheckDefined(values);
356  }
357  else
358  values.setNoFlow();
359  }
360 
362 
367 
371  template <class Context>
372  void initial(PrimaryVariables& values,
373  const Context& context,
374  unsigned spaceIdx,
375  unsigned timeIdx) const
376  {
377  Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
378 
379  initialFluidState_(fs, context, spaceIdx, timeIdx);
380 
381  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
382  values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true);
383  Opm::Valgrind::CheckDefined(values);
384  }
385 
392  template <class Context>
393  void source(RateVector& rate,
394  const Context& /*context*/,
395  unsigned /*spaceIdx*/,
396  unsigned /*timeIdx*/) const
397  { rate = Scalar(0.0); }
398 
400 
401 private:
402  bool onLeftBoundary_(const GlobalPosition& pos) const
403  { return pos[0] < eps_; }
404 
405  bool onRightBoundary_(const GlobalPosition& pos) const
406  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
407 
408  bool onLowerBoundary_(const GlobalPosition& pos) const
409  { return pos[1] < eps_; }
410 
411  bool onUpperBoundary_(const GlobalPosition& pos) const
412  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
413 
414  bool onInlet_(const GlobalPosition& pos) const
415  { return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
416 
417  template <class FluidState, class Context>
418  void initialFluidState_(FluidState& fs, const Context& context,
419  unsigned spaceIdx, unsigned timeIdx) const
420  {
421  const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
422  Scalar y = pos[1];
423  Scalar x = pos[0];
424 
425  Scalar densityW = 1000.0;
426  Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
427  if (pc < 0.0)
428  pc = 0.0;
429 
430  // set pressures
431  const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
432  Scalar Sw = matParams.Swr();
433  Scalar Swr = matParams.Swr();
434  Scalar Sgr = matParams.Sgr();
435  if (Sw < Swr)
436  Sw = Swr;
437  if (Sw > 1 - Sgr)
438  Sw = 1 - Sgr;
439  Scalar Sg = 1 - Sw;
440 
441  Opm::Valgrind::CheckDefined(Sw);
442  Opm::Valgrind::CheckDefined(Sg);
443 
444  fs.setSaturation(waterPhaseIdx, Sw);
445  fs.setSaturation(gasPhaseIdx, Sg);
446  fs.setSaturation(naplPhaseIdx, 0);
447 
448  // set temperature of all phases
449  fs.setTemperature(temperature_);
450 
451  // compute pressures
452  Scalar pcAll[numPhases];
453  Scalar pg = 1e5;
454  if (onLeftBoundary_(pos))
455  pg += 10e3;
456  MaterialLaw::capillaryPressures(pcAll, matParams, fs);
457  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
458  fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
459 
460  // set composition of gas phase
461  fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
462  fs.setMoleFraction(gasPhaseIdx, airIdx,
463  1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
464  fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
465 
466  using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
467  typename FluidSystem::template ParameterCache<Scalar> paramCache;
468  CFRP::solve(fs, paramCache, gasPhaseIdx,
469  /*setViscosity=*/true,
470  /*setEnthalpy=*/false);
471 
472  fs.setMoleFraction(waterPhaseIdx, H2OIdx,
473  1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
474  }
475 
476  bool isFineMaterial_(const GlobalPosition& pos) const
477  { return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
478 
479  DimMatrix fineK_;
480  DimMatrix coarseK_;
481 
482  Scalar porosity_;
483 
484  MaterialLawParams materialParams_;
485 
486  Scalar temperature_;
487  Scalar eps_;
488 };
489 } // namespace Opm
490 
491 #endif
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:151
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:393
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:296
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:286
std::string name() const
Definition: infiltrationproblem.hh:256
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:310
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:336
void endTimeStep()
Definition: infiltrationproblem.hh:266
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:198
bool shouldWriteRestartFile() const
Definition: infiltrationproblem.hh:250
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:320
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:372
void finishInit()
Definition: infiltrationproblem.hh:206
Definition: infiltrationproblem.hh:57