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