All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 
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/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 Ewoms {
49 template <class TypeTag>
51 
52 namespace Properties {
54 
55 // Use 2d YaspGrid
56 SET_TYPE_PROP(RichardsLensProblem, Grid, Dune::YaspGrid<2>);
57 
58 // Set the physical problem to be solved
60 
61 // Set the wetting phase
62 SET_PROP(RichardsLensProblem, WettingFluid)
63 {
64 private:
65  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
66 
67 public:
68  typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
69 };
70 
71 // Set the material Law
72 SET_PROP(RichardsLensProblem, MaterialLaw)
73 {
74 private:
75  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76  enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
77  enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
78 
79  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
80  typedef Opm::TwoPhaseMaterialTraits<Scalar,
81  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
82  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>
83  Traits;
84 
85  // define the material law which is parameterized by effective
86  // saturations
87  typedef Opm::RegularizedVanGenuchten<Traits> EffectiveLaw;
88 
89 public:
90  // define the material law parameterized by absolute saturations
91  typedef Opm::EffToAbsLaw<EffectiveLaw> type;
92 };
93 
94 // Enable gravitational acceleration
95 SET_BOOL_PROP(RichardsLensProblem, EnableGravity, true);
96 
97 // Use central differences to approximate the Jacobian matrix
98 SET_INT_PROP(RichardsLensProblem, NumericDifferenceMethod, 0);
99 
100 // Set the maximum number of newton iterations of a time step
101 SET_INT_PROP(RichardsLensProblem, NewtonMaxIterations, 28);
102 
103 // Set the "desireable" number of newton iterations of a time step
104 SET_INT_PROP(RichardsLensProblem, NewtonTargetIterations, 18);
105 
106 // Do not write the intermediate results of the newton method
107 SET_BOOL_PROP(RichardsLensProblem, NewtonWriteConvergence, false);
108 
109 // The default for the end time of the simulation
110 SET_SCALAR_PROP(RichardsLensProblem, EndTime, 3000);
111 
112 // The default for the initial time step size of the simulation
113 SET_SCALAR_PROP(RichardsLensProblem, InitialTimeStepSize, 100);
114 
115 // The default DGF file to load
116 SET_STRING_PROP(RichardsLensProblem, GridFile, "./data/richardslens_24x16.dgf");
117 } // namespace Properties
118 
135 template <class TypeTag>
136 class RichardsLensProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
137 {
138  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
139 
140  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
141  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
142  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
143  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
144  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
145  typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
146  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
147  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
148  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
149  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
150 
151  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
152  enum {
153  // copy some indices for convenience
154  pressureWIdx = Indices::pressureWIdx,
155  contiEqIdx = Indices::contiEqIdx,
156  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
157  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
158  numPhases = FluidSystem::numPhases,
159 
160  // Grid and world dimension
161  dimWorld = GridView::dimensionworld
162  };
163 
164  // get the material law from the property system
165  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
167  typedef typename MaterialLaw::Params MaterialLawParams;
168 
169  typedef typename GridView::ctype CoordScalar;
170  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
171  typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
172  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
173 
174 public:
178  RichardsLensProblem(Simulator& simulator)
179  : ParentType(simulator)
180  , pnRef_(1e5)
181  {
182  dofIsInLens_.resize(simulator.model().numGridDof());
183  }
184 
188  void finishInit()
189  {
190  ParentType::finishInit();
191 
192  eps_ = 3e-6;
193  pnRef_ = 1e5;
194 
195  lensLowerLeft_[0] = 1.0;
196  lensLowerLeft_[1] = 2.0;
197 
198  lensUpperRight_[0] = 4.0;
199  lensUpperRight_[1] = 3.0;
200 
201  // parameters for the Van Genuchten law
202  // alpha and n
203  lensMaterialParams_.setVgAlpha(0.00045);
204  lensMaterialParams_.setVgN(7.3);
205  lensMaterialParams_.finalize();
206 
207  outerMaterialParams_.setVgAlpha(0.0037);
208  outerMaterialParams_.setVgN(4.7);
209  outerMaterialParams_.finalize();
210 
211  // parameters for the linear law
212  // minimum and maximum pressures
213  // lensMaterialParams_.setEntryPC(0);
214  // outerMaterialParams_.setEntryPC(0);
215  // lensMaterialParams_.setMaxPC(0);
216  // outerMaterialParams_.setMaxPC(0);
217 
218  lensK_ = this->toDimMatrix_(1e-12);
219  outerK_ = this->toDimMatrix_(5e-12);
220 
221  // determine which degrees of freedom are in the lens
222  Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
223  auto elemIt = this->gridView().template begin</*codim=*/0>();
224  auto elemEndIt = this->gridView().template end</*codim=*/0>();
225  for (; elemIt != elemEndIt; ++elemIt) {
226  stencil.update(*elemIt);
227  for (unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
228  unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
229  const auto& dofPos = stencil.subControlVolume(dofIdx).center();
230  dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
231  }
232  }
233  }
234 
238 
243  std::string name() const
244  {
245  std::ostringstream oss;
246  oss << "lens_richards_"
247  << Model::discretizationName();
248  return oss.str();
249  }
250 
254  void endTimeStep()
255  {
256 #ifndef NDEBUG
257  this->model().checkConservativeness();
258 
259  // Calculate storage terms
260  EqVector storage;
261  this->model().globalStorage(storage);
262 
263  // Write mass balance information for rank 0
264  if (this->gridView().comm().rank() == 0) {
265  std::cout << "Storage: " << storage << std::endl << std::flush;
266  }
267 #endif // NDEBUG
268  }
269 
273  template <class Context>
274  Scalar temperature(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
275  { return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
276 
277  Scalar temperature(unsigned globalSpaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
278  { return 273.15 + 10; } // -> 10°C
279 
283  template <class Context>
284  const DimMatrix& intrinsicPermeability(const Context& context,
285  unsigned spaceIdx,
286  unsigned timeIdx) const
287  {
288  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
289  if (isInLens_(pos))
290  return lensK_;
291  return outerK_;
292  }
293 
297  template <class Context>
298  Scalar porosity(const Context& context OPM_UNUSED,
299  unsigned spaceIdx OPM_UNUSED,
300  unsigned timeIdx OPM_UNUSED) const
301  { return 0.4; }
302 
306  template <class Context>
307  const MaterialLawParams& materialLawParams(const Context& context,
308  unsigned spaceIdx,
309  unsigned timeIdx) const
310  {
311  unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
312  return materialLawParams(globalSpaceIdx, timeIdx);
313  }
314 
315  const MaterialLawParams& materialLawParams(unsigned globalSpaceIdx,
316  unsigned timeIdx OPM_UNUSED) const
317  {
318  if (dofIsInLens_[globalSpaceIdx])
319  return lensMaterialParams_;
320  return outerMaterialParams_;
321  }
322 
328  template <class Context>
329  Scalar referencePressure(const Context& context,
330  unsigned spaceIdx,
331  unsigned timeIdx) const
332  { return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
333 
334  // the Richards model does not have an element context available at all places
335  // where the reference pressure is required...
336  Scalar referencePressure(unsigned globalSpaceIdx OPM_UNUSED,
337  unsigned timeIdx OPM_UNUSED) const
338  { return pnRef_; }
339 
341 
345 
350  template <class Context>
351  void boundary(BoundaryRateVector& values,
352  const Context& context,
353  unsigned spaceIdx,
354  unsigned timeIdx) const
355  {
356  const auto& pos = context.pos(spaceIdx, timeIdx);
357 
358  if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
359  const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
360 
361  Scalar Sw = 0.0;
362  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
363  fs.setSaturation(wettingPhaseIdx, Sw);
364  fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
365 
366  PhaseVector pC;
367  MaterialLaw::capillaryPressures(pC, materialParams, fs);
368  fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
369  fs.setPressure(nonWettingPhaseIdx, pnRef_);
370 
371  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
372  }
373  else if (onInlet_(pos)) {
374  RateVector massRate(0.0);
375 
376  // inflow of water
377  massRate[contiEqIdx] = -0.04; // kg / (m * s)
378 
379  values.setMassRate(massRate);
380  }
381  else
382  values.setNoFlow();
383  }
384 
386 
390 
395  template <class Context>
396  void initial(PrimaryVariables& values,
397  const Context& context,
398  unsigned spaceIdx,
399  unsigned timeIdx) const
400  {
401  const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
402 
403  Scalar Sw = 0.0;
404  Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
405  fs.setSaturation(wettingPhaseIdx, Sw);
406  fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
407 
408  PhaseVector pC;
409  MaterialLaw::capillaryPressures(pC, materialParams, fs);
410  values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
411  }
412 
419  template <class Context>
420  void source(RateVector& rate,
421  const Context& context OPM_UNUSED,
422  unsigned spaceIdx OPM_UNUSED,
423  unsigned timeIdx OPM_UNUSED) const
424  { rate = Scalar(0.0); }
425 
427 
428 private:
429  bool onLeftBoundary_(const GlobalPosition& pos) const
430  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
431 
432  bool onRightBoundary_(const GlobalPosition& pos) const
433  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
434 
435  bool onLowerBoundary_(const GlobalPosition& pos) const
436  { return pos[1] < this->boundingBoxMin()[1] + eps_; }
437 
438  bool onUpperBoundary_(const GlobalPosition& pos) const
439  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
440 
441  bool onInlet_(const GlobalPosition& pos) const
442  {
443  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
444  Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
445  return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
446  }
447 
448  bool isInLens_(const GlobalPosition& pos) const
449  {
450  for (unsigned i = 0; i < dimWorld; ++i) {
451  if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
452  return false;
453  }
454  return true;
455  }
456 
457  GlobalPosition lensLowerLeft_;
458  GlobalPosition lensUpperRight_;
459 
460  DimMatrix lensK_;
461  DimMatrix outerK_;
462  MaterialLawParams lensMaterialParams_;
463  MaterialLawParams outerMaterialParams_;
464 
465  std::vector<bool> dofIsInLens_;
466 
467  Scalar eps_;
468  Scalar pnRef_;
469 };
470 } // namespace Ewoms
471 
472 #endif
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
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: richardslensproblem.hh:420
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:298
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: richardslensproblem.hh:188
#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
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:274
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
void endTimeStep()
Called by the simulator after each time integration.
Definition: richardslensproblem.hh:254
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:329
std::string name() const
The problem name.
Definition: richardslensproblem.hh:243
This model implements a variant of the Richards equation for quasi-twophase flow. ...
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:284
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:307
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain...
Definition: richardslensproblem.hh:50
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: richardslensproblem.hh:396
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#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: richardslensproblem.hh:351
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394