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