My Project
richardslensproblem.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_RICHARDS_LENS_PROBLEM_HH
29 #define EWOMS_RICHARDS_LENS_PROBLEM_HH
30 
31 #include <opm/models/richards/richardsmodel.hh>
32 
33 #include <opm/material/components/SimpleH2O.hpp>
34 #include <opm/material/fluidsystems/LiquidPhase.hpp>
35 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
36 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
39 #include <opm/material/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 namespace Opm {
49 template <class TypeTag>
50 class RichardsLensProblem;
51 
52 } // namespace Opm
53 
54 namespace Opm::Properties {
55 
56 // Create new type tags
57 namespace TTag {
58 struct RichardsLensProblem { using InheritsFrom = std::tuple<Richards>; };
59 } // end namespace TTag
60 
61 // Use 2d YaspGrid
62 template<class TypeTag>
63 struct Grid<TypeTag, TTag::RichardsLensProblem> { using type = Dune::YaspGrid<2>; };
64 
65 // Set the physical problem to be solved
66 template<class TypeTag>
67 struct Problem<TypeTag, TTag::RichardsLensProblem> { using type = Opm::RichardsLensProblem<TypeTag>; };
68 
69 // Set the wetting phase
70 template<class TypeTag>
71 struct WettingFluid<TypeTag, TTag::RichardsLensProblem>
72 {
73 private:
74  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
75 
76 public:
77  using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
78 };
79 
80 // Set the material Law
81 template<class TypeTag>
82 struct MaterialLaw<TypeTag, TTag::RichardsLensProblem>
83 {
84 private:
85  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
86  enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
87  enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
88 
89  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
90  using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
91  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
92  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
93 
94  // define the material law which is parameterized by effective
95  // saturations
96  using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
97 
98 public:
99  // define the material law parameterized by absolute saturations
100  using type = Opm::EffToAbsLaw<EffectiveLaw>;
101 };
102 
103 // Enable gravitational acceleration
104 template<class TypeTag>
105 struct EnableGravity<TypeTag, TTag::RichardsLensProblem> { static constexpr bool value = true; };
106 
107 // Use central differences to approximate the Jacobian matrix
108 template<class TypeTag>
109 struct NumericDifferenceMethod<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 0; };
110 
111 // Set the maximum number of newton iterations of a time step
112 template<class TypeTag>
113 struct NewtonMaxIterations<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 28; };
114 
115 // Set the "desireable" number of newton iterations of a time step
116 template<class TypeTag>
117 struct NewtonTargetIterations<TypeTag, TTag::RichardsLensProblem> { static constexpr int value = 18; };
118 
119 // Do not write the intermediate results of the newton method
120 template<class TypeTag>
121 struct NewtonWriteConvergence<TypeTag, TTag::RichardsLensProblem> { static constexpr bool value = false; };
122 
123 // The default for the end time of the simulation
124 template<class TypeTag>
125 struct EndTime<TypeTag, TTag::RichardsLensProblem>
126 {
127  using type = GetPropType<TypeTag, Scalar>;
128  static constexpr type value = 3000;
129 };
130 
131 // The default for the initial time step size of the simulation
132 template<class TypeTag>
133 struct InitialTimeStepSize<TypeTag, TTag::RichardsLensProblem>
134 {
135  using type = GetPropType<TypeTag, Scalar>;
136  static constexpr type value = 100;
137 };
138 
139 // The default DGF file to load
140 template<class TypeTag>
141 struct GridFile<TypeTag, TTag::RichardsLensProblem> { static constexpr auto value = "./data/richardslens_24x16.dgf"; };
142 
143 } // namespace Opm::Properties
144 
145 namespace Opm {
146 
163 template <class TypeTag>
164 class RichardsLensProblem : public GetPropType<TypeTag, Properties::BaseProblem>
165 {
166  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
167 
168  using GridView = GetPropType<TypeTag, Properties::GridView>;
169  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
170  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
171  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
172  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
173  using Stencil = GetPropType<TypeTag, Properties::Stencil>;
174  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
175  using Model = GetPropType<TypeTag, Properties::Model>;
176  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
177  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
178 
179  using Indices = GetPropType<TypeTag, Properties::Indices>;
180  enum {
181  // copy some indices for convenience
182  pressureWIdx = Indices::pressureWIdx,
183  contiEqIdx = Indices::contiEqIdx,
184  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
185  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
186  numPhases = FluidSystem::numPhases,
187 
188  // Grid and world dimension
189  dimWorld = GridView::dimensionworld
190  };
191 
192  // get the material law from the property system
193  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
195  using MaterialLawParams = typename MaterialLaw::Params;
196 
197  using CoordScalar = typename GridView::ctype;
198  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
199  using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
200  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
201 
202 public:
206  RichardsLensProblem(Simulator& simulator)
207  : ParentType(simulator)
208  , pnRef_(1e5)
209  {
210  dofIsInLens_.resize(simulator.model().numGridDof());
211  }
212 
216  void finishInit()
217  {
218  ParentType::finishInit();
219 
220  eps_ = 3e-6;
221  pnRef_ = 1e5;
222 
223  lensLowerLeft_[0] = 1.0;
224  lensLowerLeft_[1] = 2.0;
225 
226  lensUpperRight_[0] = 4.0;
227  lensUpperRight_[1] = 3.0;
228 
229  // parameters for the Van Genuchten law
230  // alpha and n
231  lensMaterialParams_.setVgAlpha(0.00045);
232  lensMaterialParams_.setVgN(7.3);
233  lensMaterialParams_.finalize();
234 
235  outerMaterialParams_.setVgAlpha(0.0037);
236  outerMaterialParams_.setVgN(4.7);
237  outerMaterialParams_.finalize();
238 
239  // parameters for the linear law
240  // minimum and maximum pressures
241  // lensMaterialParams_.setEntryPC(0);
242  // outerMaterialParams_.setEntryPC(0);
243  // lensMaterialParams_.setMaxPC(0);
244  // outerMaterialParams_.setMaxPC(0);
245 
246  lensK_ = this->toDimMatrix_(1e-12);
247  outerK_ = this->toDimMatrix_(5e-12);
248 
249  // determine which degrees of freedom are in the lens
250  Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
251  auto elemIt = this->gridView().template begin</*codim=*/0>();
252  auto elemEndIt = this->gridView().template end</*codim=*/0>();
253  for (; elemIt != elemEndIt; ++elemIt) {
254  stencil.update(*elemIt);
255  for (unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
256  unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
257  const auto& dofPos = stencil.subControlVolume(dofIdx).center();
258  dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
259  }
260  }
261  }
262 
267 
271  std::string name() const
272  {
273  std::ostringstream oss;
274  oss << "lens_richards_"
275  << Model::discretizationName();
276  return oss.str();
277  }
278 
282  void endTimeStep()
283  {
284 #ifndef NDEBUG
285  this->model().checkConservativeness();
286 
287  // Calculate storage terms
288  EqVector storage;
289  this->model().globalStorage(storage);
290 
291  // Write mass balance information for rank 0
292  if (this->gridView().comm().rank() == 0) {
293  std::cout << "Storage: " << storage << std::endl << std::flush;
294  }
295 #endif // NDEBUG
296  }
297 
301  template <class Context>
302  Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
303  { return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
304 
305  Scalar temperature(unsigned globalSpaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
306  { return 273.15 + 10; } // -> 10°C
307 
311  template <class Context>
312  const DimMatrix& intrinsicPermeability(const Context& context,
313  unsigned spaceIdx,
314  unsigned timeIdx) const
315  {
316  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
317  if (isInLens_(pos))
318  return lensK_;
319  return outerK_;
320  }
321 
325  template <class Context>
326  Scalar porosity(const Context& context OPM_UNUSED,
327  unsigned spaceIdx OPM_UNUSED,
328  unsigned timeIdx OPM_UNUSED) const
329  { return 0.4; }
330 
334  template <class Context>
335  const MaterialLawParams& materialLawParams(const Context& context,
336  unsigned spaceIdx,
337  unsigned timeIdx) const
338  {
339  unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
340  return materialLawParams(globalSpaceIdx, timeIdx);
341  }
342 
343  const MaterialLawParams& materialLawParams(unsigned globalSpaceIdx,
344  unsigned timeIdx OPM_UNUSED) const
345  {
346  if (dofIsInLens_[globalSpaceIdx])
347  return lensMaterialParams_;
348  return outerMaterialParams_;
349  }
350 
356  template <class Context>
357  Scalar referencePressure(const Context& context,
358  unsigned spaceIdx,
359  unsigned timeIdx) const
360  { return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
361 
362  // the Richards model does not have an element context available at all places
363  // where the reference pressure is required...
364  Scalar referencePressure(unsigned globalSpaceIdx OPM_UNUSED,
365  unsigned timeIdx OPM_UNUSED) const
366  { return pnRef_; }
367 
369 
374 
378  template <class Context>
379  void boundary(BoundaryRateVector& values,
380  const Context& context,
381  unsigned spaceIdx,
382  unsigned timeIdx) const
383  {
384  const auto& pos = context.pos(spaceIdx, timeIdx);
385 
386  if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
387  const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
388 
389  Scalar Sw = 0.0;
390  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
391  fs.setSaturation(wettingPhaseIdx, Sw);
392  fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
393 
394  PhaseVector pC;
395  MaterialLaw::capillaryPressures(pC, materialParams, fs);
396  fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
397  fs.setPressure(nonWettingPhaseIdx, pnRef_);
398 
399  typename FluidSystem::template ParameterCache<Scalar> paramCache;
400  paramCache.updateAll(fs);
401  fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx));
402  //fs.setDensity(nonWettingPhaseIdx, FluidSystem::density(fs, paramCache, nonWettingPhaseIdx));
403 
404  fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx));
405  //fs.setViscosity(nonWettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, nonWettingPhaseIdx));
406 
407  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
408  }
409  else if (onInlet_(pos)) {
410  RateVector massRate(0.0);
411 
412  // inflow of water
413  massRate[contiEqIdx] = -0.04; // kg / (m * s)
414 
415  values.setMassRate(massRate);
416  }
417  else
418  values.setNoFlow();
419  }
420 
422 
427 
431  template <class Context>
432  void initial(PrimaryVariables& values,
433  const Context& context,
434  unsigned spaceIdx,
435  unsigned timeIdx) const
436  {
437  const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
438 
439  Scalar Sw = 0.0;
440  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
441  fs.setSaturation(wettingPhaseIdx, Sw);
442  fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
443 
444  PhaseVector pC;
445  MaterialLaw::capillaryPressures(pC, materialParams, fs);
446  values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
447  }
448 
455  template <class Context>
456  void source(RateVector& rate,
457  const Context& context OPM_UNUSED,
458  unsigned spaceIdx OPM_UNUSED,
459  unsigned timeIdx OPM_UNUSED) const
460  { rate = Scalar(0.0); }
461 
463 
464 private:
465  bool onLeftBoundary_(const GlobalPosition& pos) const
466  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
467 
468  bool onRightBoundary_(const GlobalPosition& pos) const
469  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
470 
471  bool onLowerBoundary_(const GlobalPosition& pos) const
472  { return pos[1] < this->boundingBoxMin()[1] + eps_; }
473 
474  bool onUpperBoundary_(const GlobalPosition& pos) const
475  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
476 
477  bool onInlet_(const GlobalPosition& pos) const
478  {
479  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
480  Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
481  return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
482  }
483 
484  bool isInLens_(const GlobalPosition& pos) const
485  {
486  for (unsigned i = 0; i < dimWorld; ++i) {
487  if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
488  return false;
489  }
490  return true;
491  }
492 
493  GlobalPosition lensLowerLeft_;
494  GlobalPosition lensUpperRight_;
495 
496  DimMatrix lensK_;
497  DimMatrix outerK_;
498  MaterialLawParams lensMaterialParams_;
499  MaterialLawParams outerMaterialParams_;
500 
501  std::vector<bool> dofIsInLens_;
502 
503  Scalar eps_;
504  Scalar pnRef_;
505 };
506 } // namespace Opm
507 
508 #endif
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain.
Definition: richardslensproblem.hh:165
void finishInit()
Definition: richardslensproblem.hh:216
std::string name() const
Definition: richardslensproblem.hh:271
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:326
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:357
RichardsLensProblem(Simulator &simulator)
Definition: richardslensproblem.hh:206
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:379
void endTimeStep()
Definition: richardslensproblem.hh:282
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:456
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:312
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:302
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:432
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:335
Definition: richardslensproblem.hh:58