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