My Project
groundwaterproblem.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_GROUND_WATER_PROBLEM_HH
29 #define EWOMS_GROUND_WATER_PROBLEM_HH
30 
31 #include <opm/models/immiscible/immiscibleproperties.hh>
32 #include <opm/simulators/linalg/parallelistlbackend.hh>
33 
34 #include <opm/material/components/SimpleH2O.hpp>
35 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
36 #include <opm/material/fluidsystems/LiquidPhase.hpp>
37 #include <opm/material/common/Unused.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/fmatrix.hh>
44 #include <dune/common/fvector.hh>
45 
46 #include <sstream>
47 #include <string>
48 
49 namespace Opm {
50 template <class TypeTag>
51 class GroundWaterProblem;
52 }
53 
54 namespace Opm::Properties {
55 
56 namespace TTag {
58 }
59 
60 template<class TypeTag, class MyTypeTag>
61 struct LensLowerLeftX { using type = UndefinedProperty; };
62 template<class TypeTag, class MyTypeTag>
63 struct LensLowerLeftY { using type = UndefinedProperty; };
64 template<class TypeTag, class MyTypeTag>
65 struct LensLowerLeftZ { using type = UndefinedProperty; };
66 template<class TypeTag, class MyTypeTag>
67 struct LensUpperRightX { using type = UndefinedProperty; };
68 template<class TypeTag, class MyTypeTag>
69 struct LensUpperRightY { using type = UndefinedProperty; };
70 template<class TypeTag, class MyTypeTag>
71 struct LensUpperRightZ { using type = UndefinedProperty; };
72 template<class TypeTag, class MyTypeTag>
73 struct Permeability { using type = UndefinedProperty; };
74 template<class TypeTag, class MyTypeTag>
75 struct PermeabilityLens { using type = UndefinedProperty; };
76 
77 template<class TypeTag>
78 struct Fluid<TypeTag, TTag::GroundWaterBaseProblem>
79 {
80 private:
81  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
82 
83 public:
84  using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
85 };
86 
87 // Set the grid type
88 template<class TypeTag>
89 struct Grid<TypeTag, TTag::GroundWaterBaseProblem> { using type = Dune::YaspGrid<2>; };
90 // struct Grid<TypeTag, TTag::GroundWaterBaseProblem> { using type = Dune::SGrid<2, 2>; };
91 
92 template<class TypeTag>
93 struct Problem<TypeTag, TTag::GroundWaterBaseProblem>
95 
96 template<class TypeTag>
97 struct LensLowerLeftX<TypeTag, TTag::GroundWaterBaseProblem>
98 {
99  using type = GetPropType<TypeTag, Scalar>;
100  static constexpr type value = 0.25;
101 };
102 template<class TypeTag>
103 struct LensLowerLeftY<TypeTag, TTag::GroundWaterBaseProblem>
104 {
105  using type = GetPropType<TypeTag, Scalar>;
106  static constexpr type value = 0.25;
107 };
108 template<class TypeTag>
109 struct LensLowerLeftZ<TypeTag, TTag::GroundWaterBaseProblem>
110 {
111  using type = GetPropType<TypeTag, Scalar>;
112  static constexpr type value = 0.25;
113 };
114 template<class TypeTag>
115 struct LensUpperRightX<TypeTag, TTag::GroundWaterBaseProblem>
116 {
117  using type = GetPropType<TypeTag, Scalar>;
118  static constexpr type value = 0.75;
119 };
120 template<class TypeTag>
121 struct LensUpperRightY<TypeTag, TTag::GroundWaterBaseProblem>
122 {
123  using type = GetPropType<TypeTag, Scalar>;
124  static constexpr type value = 0.75;
125 };
126 template<class TypeTag>
127 struct LensUpperRightZ<TypeTag, TTag::GroundWaterBaseProblem>
128 {
129  using type = GetPropType<TypeTag, Scalar>;
130  static constexpr type value = 0.75;
131 };
132 template<class TypeTag>
133 struct Permeability<TypeTag, TTag::GroundWaterBaseProblem>
134 {
135  using type = GetPropType<TypeTag, Scalar>;
136  static constexpr type value = 1e-10;
137 };
138 template<class TypeTag>
139 struct PermeabilityLens<TypeTag, TTag::GroundWaterBaseProblem>
140 {
141  using type = GetPropType<TypeTag, Scalar>;
142  static constexpr type value = 1e-12;
143 };
144 
145 // Enable gravity
146 template<class TypeTag>
147 struct EnableGravity<TypeTag, TTag::GroundWaterBaseProblem> { static constexpr bool value = true; };
148 
149 // The default for the end time of the simulation
150 template<class TypeTag>
151 struct EndTime<TypeTag, TTag::GroundWaterBaseProblem>
152 {
153  using type = GetPropType<TypeTag, Scalar>;
154  static constexpr type value = 1;
155 };
156 
157 // The default for the initial time step size of the simulation
158 template<class TypeTag>
159 struct InitialTimeStepSize<TypeTag, TTag::GroundWaterBaseProblem>
160 {
161  using type = GetPropType<TypeTag, Scalar>;
162  static constexpr type value = 1;
163 };
164 
165 // The default DGF file to load
166 template<class TypeTag>
167 struct GridFile<TypeTag, TTag::GroundWaterBaseProblem> { static constexpr auto value = "./data/groundwater_2d.dgf"; };
168 
169 // Use the conjugated gradient linear solver with the default preconditioner (i.e.,
170 // ILU-0) from dune-istl
171 template<class TypeTag>
172 struct LinearSolverSplice<TypeTag, TTag::GroundWaterBaseProblem> { using type = TTag::ParallelIstlLinearSolver; };
173 
174 template<class TypeTag>
175 struct LinearSolverWrapper<TypeTag, TTag::GroundWaterBaseProblem>
176 { using type = Opm::Linear::SolverWrapperConjugatedGradients<TypeTag>; };
177 
178 } // namespace Opm::Properties
179 
180 namespace Opm {
193 template <class TypeTag>
194 class GroundWaterProblem : public GetPropType<TypeTag, Properties::BaseProblem>
195 {
196  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
197 
198  using GridView = GetPropType<TypeTag, Properties::GridView>;
199  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
200  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
201 
202  // copy some indices for convenience
203  using Indices = GetPropType<TypeTag, Properties::Indices>;
204  enum {
205  numPhases = FluidSystem::numPhases,
206 
207  // Grid and world dimension
208  dim = GridView::dimension,
209  dimWorld = GridView::dimensionworld,
210 
211  // indices of the primary variables
212  pressure0Idx = Indices::pressure0Idx
213  };
214 
215  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
216  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
217  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
218  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
219  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
220  using Model = GetPropType<TypeTag, Properties::Model>;
221 
222  using CoordScalar = typename GridView::ctype;
223  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
224 
225  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
226 
227 public:
231  GroundWaterProblem(Simulator& simulator)
232  : ParentType(simulator)
233  { }
234 
238  void finishInit()
239  {
240  ParentType::finishInit();
241 
242  eps_ = 1.0e-3;
243 
244  lensLowerLeft_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftX);
245  if (dim > 1)
246  lensLowerLeft_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
247  if (dim > 2)
248  lensLowerLeft_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
249 
250  lensUpperRight_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX);
251  if (dim > 1)
252  lensUpperRight_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
253  if (dim > 2)
254  lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
255 
256  intrinsicPerm_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, Permeability));
257  intrinsicPermLens_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, PermeabilityLens));
258  }
259 
263  static void registerParameters()
264  {
265  ParentType::registerParameters();
266 
267  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX,
268  "The x-coordinate of the lens' lower-left corner "
269  "[m].");
270  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX,
271  "The x-coordinate of the lens' upper-right corner "
272  "[m].");
273 
274  if (dimWorld > 1) {
275  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY,
276  "The y-coordinate of the lens' lower-left "
277  "corner [m].");
278  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY,
279  "The y-coordinate of the lens' upper-right "
280  "corner [m].");
281  }
282 
283  if (dimWorld > 2) {
284  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ,
285  "The z-coordinate of the lens' lower-left "
286  "corner [m].");
287  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ,
288  "The z-coordinate of the lens' upper-right "
289  "corner [m].");
290  }
291 
292  EWOMS_REGISTER_PARAM(TypeTag, Scalar, Permeability,
293  "The intrinsic permeability [m^2] of the ambient "
294  "material.");
295  EWOMS_REGISTER_PARAM(TypeTag, Scalar, PermeabilityLens,
296  "The intrinsic permeability [m^2] of the lens.");
297  }
298 
302  // \{
303 
307  std::string name() const
308  {
309  std::ostringstream oss;
310  oss << "groundwater_" << Model::name();
311  return oss.str();
312  }
313 
317  void endTimeStep()
318  {
319 #ifndef NDEBUG
320  this->model().checkConservativeness();
321 
322  // Calculate storage terms
323  EqVector storage;
324  this->model().globalStorage(storage);
325 
326  // Write mass balance information for rank 0
327  if (this->gridView().comm().rank() == 0) {
328  std::cout << "Storage: " << storage << std::endl << std::flush;
329  }
330 #endif // NDEBUG
331  }
332 
336  template <class Context>
337  Scalar temperature(const Context& context OPM_UNUSED,
338  unsigned spaceIdx OPM_UNUSED,
339  unsigned timeIdx OPM_UNUSED) const
340  { return 273.15 + 10; } // 10C
341 
345  template <class Context>
346  Scalar porosity(const Context& context OPM_UNUSED,
347  unsigned spaceIdx OPM_UNUSED,
348  unsigned timeIdx OPM_UNUSED) const
349  { return 0.4; }
350 
354  template <class Context>
355  const DimMatrix& intrinsicPermeability(const Context& context,
356  unsigned spaceIdx,
357  unsigned timeIdx) const
358  {
359  if (isInLens_(context.pos(spaceIdx, timeIdx)))
360  return intrinsicPermLens_;
361  else
362  return intrinsicPerm_;
363  }
364 
366 
370 
374  template <class Context>
375  void boundary(BoundaryRateVector& values, const Context& context,
376  unsigned spaceIdx, unsigned timeIdx) const
377  {
378  const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
379 
380  if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)) {
381  Scalar pressure;
382  Scalar T = temperature(context, spaceIdx, timeIdx);
383  if (onLowerBoundary_(globalPos))
384  pressure = 2e5;
385  else // on upper boundary
386  pressure = 1e5;
387 
388  Opm::ImmiscibleFluidState<Scalar, FluidSystem,
389  /*storeEnthalpy=*/false> fs;
390  fs.setSaturation(/*phaseIdx=*/0, 1.0);
391  fs.setPressure(/*phaseIdx=*/0, pressure);
392  fs.setTemperature(T);
393 
394  typename FluidSystem::template ParameterCache<Scalar> paramCache;
395  paramCache.updateAll(fs);
396  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
397  fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
398  fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
399  }
400 
401  // impose an freeflow boundary condition
402  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
403  }
404  else {
405  // no flow boundary
406  values.setNoFlow();
407  }
408  }
409 
411 
416 
420  template <class Context>
421  void initial(PrimaryVariables& values,
422  const Context& context OPM_UNUSED,
423  unsigned spaceIdx OPM_UNUSED,
424  unsigned timeIdx OPM_UNUSED) const
425  {
426  // const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
427  values[pressure0Idx] = 1.0e+5; // + 9.81*1.23*(20-globalPos[dim-1]);
428  }
429 
433  template <class Context>
434  void source(RateVector& rate,
435  const Context& context OPM_UNUSED,
436  unsigned spaceIdx OPM_UNUSED,
437  unsigned timeIdx OPM_UNUSED) const
438  { rate = Scalar(0.0); }
439 
441 
442 private:
443  bool onLowerBoundary_(const GlobalPosition& pos) const
444  { return pos[dim - 1] < eps_; }
445 
446  bool onUpperBoundary_(const GlobalPosition& pos) const
447  { return pos[dim - 1] > this->boundingBoxMax()[dim - 1] - eps_; }
448 
449  bool isInLens_(const GlobalPosition& pos) const
450  {
451  return lensLowerLeft_[0] <= pos[0] && pos[0] <= lensUpperRight_[0]
452  && lensLowerLeft_[1] <= pos[1] && pos[1] <= lensUpperRight_[1];
453  }
454 
455  GlobalPosition lensLowerLeft_;
456  GlobalPosition lensUpperRight_;
457 
458  DimMatrix intrinsicPerm_;
459  DimMatrix intrinsicPermLens_;
460 
461  Scalar eps_;
462 };
463 } // namespace Opm
464 
465 #endif
Test for the immisicible VCVF discretization with only a single phase.
Definition: groundwaterproblem.hh:195
void endTimeStep()
Definition: groundwaterproblem.hh:317
GroundWaterProblem(Simulator &simulator)
Definition: groundwaterproblem.hh:231
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:434
static void registerParameters()
Definition: groundwaterproblem.hh:263
std::string name() const
Definition: groundwaterproblem.hh:307
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:421
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:346
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: groundwaterproblem.hh:375
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:337
void finishInit()
Definition: groundwaterproblem.hh:238
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: groundwaterproblem.hh:355
Definition: groundwaterproblem.hh:61
Definition: groundwaterproblem.hh:63
Definition: groundwaterproblem.hh:65
Definition: groundwaterproblem.hh:67
Definition: groundwaterproblem.hh:69
Definition: groundwaterproblem.hh:71
Definition: groundwaterproblem.hh:75
Definition: groundwaterproblem.hh:73
Definition: groundwaterproblem.hh:57