All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 
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/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 Ewoms {
50 template <class TypeTag>
52 }
53 
54 namespace Ewoms {
55 namespace Properties {
56 NEW_TYPE_TAG(GroundWaterBaseProblem);
57 
58 NEW_PROP_TAG(LensLowerLeftX);
59 NEW_PROP_TAG(LensLowerLeftY);
60 NEW_PROP_TAG(LensLowerLeftZ);
61 NEW_PROP_TAG(LensUpperRightX);
62 NEW_PROP_TAG(LensUpperRightY);
63 NEW_PROP_TAG(LensUpperRightZ);
64 NEW_PROP_TAG(Permeability);
65 NEW_PROP_TAG(PermeabilityLens);
66 
67 SET_PROP(GroundWaterBaseProblem, Fluid)
68 {
69 private:
70  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
71 
72 public:
73  typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
74 };
75 
76 // Set the grid type
77 SET_TYPE_PROP(GroundWaterBaseProblem, Grid, Dune::YaspGrid<2>);
78 // SET_TYPE_PROP(GroundWaterBaseProblem, Grid, Dune::SGrid<2, 2>);
79 
80 SET_TYPE_PROP(GroundWaterBaseProblem, Problem,
81  Ewoms::GroundWaterProblem<TypeTag>);
82 
83 SET_SCALAR_PROP(GroundWaterBaseProblem, LensLowerLeftX, 0.25);
84 SET_SCALAR_PROP(GroundWaterBaseProblem, LensLowerLeftY, 0.25);
85 SET_SCALAR_PROP(GroundWaterBaseProblem, LensLowerLeftZ, 0.25);
86 SET_SCALAR_PROP(GroundWaterBaseProblem, LensUpperRightX, 0.75);
87 SET_SCALAR_PROP(GroundWaterBaseProblem, LensUpperRightY, 0.75);
88 SET_SCALAR_PROP(GroundWaterBaseProblem, LensUpperRightZ, 0.75);
89 SET_SCALAR_PROP(GroundWaterBaseProblem, Permeability, 1e-10);
90 SET_SCALAR_PROP(GroundWaterBaseProblem, PermeabilityLens, 1e-12);
91 
92 // Enable gravity
93 SET_BOOL_PROP(GroundWaterBaseProblem, EnableGravity, true);
94 
95 // The default for the end time of the simulation
96 SET_SCALAR_PROP(GroundWaterBaseProblem, EndTime, 1);
97 
98 // The default for the initial time step size of the simulation
99 SET_SCALAR_PROP(GroundWaterBaseProblem, InitialTimeStepSize, 1);
100 
101 // The default DGF file to load
102 SET_STRING_PROP(GroundWaterBaseProblem, GridFile, "./data/groundwater_2d.dgf");
103 
104 // Use the conjugated gradient linear solver with the default preconditioner (i.e.,
105 // ILU-0) from dune-istl
106 SET_TAG_PROP(GroundWaterBaseProblem, LinearSolverSplice, ParallelIstlLinearSolver);
107 SET_TYPE_PROP(GroundWaterBaseProblem, LinearSolverWrapper,
108  Ewoms::Linear::SolverWrapperConjugatedGradients<TypeTag>);
109 }} // namespace Properties, Ewoms
110 
111 namespace Ewoms {
124 template <class TypeTag>
125 class GroundWaterProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
126 {
127  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
128 
129  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
130  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
131 
132  // copy some indices for convenience
133  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
134  enum {
135  // Grid and world dimension
136  dim = GridView::dimension,
137  dimWorld = GridView::dimensionworld,
138 
139  // indices of the primary variables
140  pressure0Idx = Indices::pressure0Idx
141  };
142 
143  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
144  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
145  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
146  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
147  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
148  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
149  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
150 
151  typedef typename GridView::ctype CoordScalar;
152  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
153 
154  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
155 
156 public:
160  GroundWaterProblem(Simulator& simulator)
161  : ParentType(simulator)
162  { }
163 
167  void finishInit()
168  {
169  ParentType::finishInit();
170 
171  eps_ = 1.0e-3;
172 
173  lensLowerLeft_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftX);
174  if (dim > 1)
175  lensLowerLeft_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
176  if (dim > 2)
177  lensLowerLeft_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
178 
179  lensUpperRight_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX);
180  if (dim > 1)
181  lensUpperRight_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
182  if (dim > 2)
183  lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
184 
185  intrinsicPerm_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, Permeability));
186  intrinsicPermLens_ = this->toDimMatrix_(EWOMS_GET_PARAM(TypeTag, Scalar, PermeabilityLens));
187  }
188 
192  static void registerParameters()
193  {
194  ParentType::registerParameters();
195 
196  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX,
197  "The x-coordinate of the lens' lower-left corner "
198  "[m].");
199  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX,
200  "The x-coordinate of the lens' upper-right corner "
201  "[m].");
202 
203  if (dimWorld > 1) {
204  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY,
205  "The y-coordinate of the lens' lower-left "
206  "corner [m].");
207  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY,
208  "The y-coordinate of the lens' upper-right "
209  "corner [m].");
210  }
211 
212  if (dimWorld > 2) {
213  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ,
214  "The z-coordinate of the lens' lower-left "
215  "corner [m].");
216  EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ,
217  "The z-coordinate of the lens' upper-right "
218  "corner [m].");
219  }
220 
221  EWOMS_REGISTER_PARAM(TypeTag, Scalar, Permeability,
222  "The intrinsic permeability [m^2] of the ambient "
223  "material.");
224  EWOMS_REGISTER_PARAM(TypeTag, Scalar, PermeabilityLens,
225  "The intrinsic permeability [m^2] of the lens.");
226  }
227 
231  // \{
232 
236  std::string name() const
237  {
238  std::ostringstream oss;
239  oss << "groundwater_" << Model::name();
240  return oss.str();
241  }
242 
246  void endTimeStep()
247  {
248 #ifndef NDEBUG
249  this->model().checkConservativeness();
250 
251  // Calculate storage terms
252  EqVector storage;
253  this->model().globalStorage(storage);
254 
255  // Write mass balance information for rank 0
256  if (this->gridView().comm().rank() == 0) {
257  std::cout << "Storage: " << storage << std::endl << std::flush;
258  }
259 #endif // NDEBUG
260  }
261 
265  template <class Context>
266  Scalar temperature(const Context& context OPM_UNUSED,
267  unsigned spaceIdx OPM_UNUSED,
268  unsigned timeIdx OPM_UNUSED) const
269  { return 273.15 + 10; } // 10C
270 
274  template <class Context>
275  Scalar porosity(const Context& context OPM_UNUSED,
276  unsigned spaceIdx OPM_UNUSED,
277  unsigned timeIdx OPM_UNUSED) const
278  { return 0.4; }
279 
283  template <class Context>
284  const DimMatrix& intrinsicPermeability(const Context& context,
285  unsigned spaceIdx,
286  unsigned timeIdx) const
287  {
288  if (isInLens_(context.pos(spaceIdx, timeIdx)))
289  return intrinsicPermLens_;
290  else
291  return intrinsicPerm_;
292  }
293 
295 
298 
303  template <class Context>
304  void boundary(BoundaryRateVector& values, const Context& context,
305  unsigned spaceIdx, unsigned timeIdx) const
306  {
307  const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
308 
309  if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)) {
310  Scalar pressure;
311  Scalar T = temperature(context, spaceIdx, timeIdx);
312  if (onLowerBoundary_(globalPos))
313  pressure = 2e5;
314  else // on upper boundary
315  pressure = 1e5;
316 
317  Opm::ImmiscibleFluidState<Scalar, FluidSystem,
318  /*storeEnthalpy=*/false> fs;
319  fs.setSaturation(/*phaseIdx=*/0, 1.0);
320  fs.setPressure(/*phaseIdx=*/0, pressure);
321  fs.setTemperature(T);
322 
323  // impose an freeflow boundary condition
324  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
325  }
326  else {
327  // no flow boundary
328  values.setNoFlow();
329  }
330  }
331 
333 
337 
342  template <class Context>
343  void initial(PrimaryVariables& values,
344  const Context& context OPM_UNUSED,
345  unsigned spaceIdx OPM_UNUSED,
346  unsigned timeIdx OPM_UNUSED) const
347  {
348  // const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
349  values[pressure0Idx] = 1.0e+5; // + 9.81*1.23*(20-globalPos[dim-1]);
350  }
351 
355  template <class Context>
356  void source(RateVector& rate,
357  const Context& context OPM_UNUSED,
358  unsigned spaceIdx OPM_UNUSED,
359  unsigned timeIdx OPM_UNUSED) const
360  { rate = Scalar(0.0); }
361 
363 
364 private:
365  bool onLowerBoundary_(const GlobalPosition& pos) const
366  { return pos[dim - 1] < eps_; }
367 
368  bool onUpperBoundary_(const GlobalPosition& pos) const
369  { return pos[dim - 1] > this->boundingBoxMax()[dim - 1] - eps_; }
370 
371  bool isInLens_(const GlobalPosition& pos) const
372  {
373  return lensLowerLeft_[0] <= pos[0] && pos[0] <= lensUpperRight_[0]
374  && lensLowerLeft_[1] <= pos[1] && pos[1] <= lensUpperRight_[1];
375  }
376 
377  GlobalPosition lensLowerLeft_;
378  GlobalPosition lensUpperRight_;
379 
380  DimMatrix intrinsicPerm_;
381  DimMatrix intrinsicPermLens_;
382 
383  Scalar eps_;
384 };
385 } // namespace Ewoms
386 
387 #endif
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:275
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: groundwaterproblem.hh:167
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: groundwaterproblem.hh:304
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the initial value for a control volume.
Definition: groundwaterproblem.hh:343
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Test for the immisicible VCVF discretization with only a single phase.
Definition: groundwaterproblem.hh:51
#define SET_TAG_PROP(EffTypeTagName, PropTagName, ValueTypeTagName)
Define a property containing a type tag.
Definition: propertysystem.hh:436
void endTimeStep()
Called by the simulator after each time integration.
Definition: groundwaterproblem.hh:246
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Provides all unmodified linear solvers from dune-istl.
std::string name() const
The problem name.
Definition: groundwaterproblem.hh:236
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: groundwaterproblem.hh:356
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Defines the properties required for the immiscible multi-phase model.
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: groundwaterproblem.hh:284
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#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
static void registerParameters()
Definition: groundwaterproblem.hh:192
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: groundwaterproblem.hh:266
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394