My Project
Loading...
Searching...
No Matches
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
48namespace Opm {
49template <class TypeTag>
50class GroundWaterProblem;
51}
52
53namespace Opm::Properties {
54
55namespace TTag {
57}
58
59template<class TypeTag, class MyTypeTag>
60struct LensLowerLeftX { using type = UndefinedProperty; };
61template<class TypeTag, class MyTypeTag>
62struct LensLowerLeftY { using type = UndefinedProperty; };
63template<class TypeTag, class MyTypeTag>
64struct LensLowerLeftZ { using type = UndefinedProperty; };
65template<class TypeTag, class MyTypeTag>
66struct LensUpperRightX { using type = UndefinedProperty; };
67template<class TypeTag, class MyTypeTag>
68struct LensUpperRightY { using type = UndefinedProperty; };
69template<class TypeTag, class MyTypeTag>
70struct LensUpperRightZ { using type = UndefinedProperty; };
71template<class TypeTag, class MyTypeTag>
72struct Permeability { using type = UndefinedProperty; };
73template<class TypeTag, class MyTypeTag>
74struct PermeabilityLens { using type = UndefinedProperty; };
75
76template<class TypeTag>
77struct Fluid<TypeTag, TTag::GroundWaterBaseProblem>
78{
79private:
80 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
81
82public:
83 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
84};
85
86// Set the grid type
87template<class TypeTag>
88struct Grid<TypeTag, TTag::GroundWaterBaseProblem> { using type = Dune::YaspGrid<2>; };
89// struct Grid<TypeTag, TTag::GroundWaterBaseProblem> { using type = Dune::SGrid<2, 2>; };
90
91template<class TypeTag>
92struct Problem<TypeTag, TTag::GroundWaterBaseProblem>
93{ using type = Opm::GroundWaterProblem<TypeTag>; };
94
95template<class TypeTag>
96struct LensLowerLeftX<TypeTag, TTag::GroundWaterBaseProblem>
97{
98 using type = GetPropType<TypeTag, Scalar>;
99 static constexpr type value = 0.25;
100};
101template<class TypeTag>
102struct LensLowerLeftY<TypeTag, TTag::GroundWaterBaseProblem>
103{
104 using type = GetPropType<TypeTag, Scalar>;
105 static constexpr type value = 0.25;
106};
107template<class TypeTag>
108struct LensLowerLeftZ<TypeTag, TTag::GroundWaterBaseProblem>
109{
110 using type = GetPropType<TypeTag, Scalar>;
111 static constexpr type value = 0.25;
112};
113template<class TypeTag>
114struct LensUpperRightX<TypeTag, TTag::GroundWaterBaseProblem>
115{
116 using type = GetPropType<TypeTag, Scalar>;
117 static constexpr type value = 0.75;
118};
119template<class TypeTag>
120struct LensUpperRightY<TypeTag, TTag::GroundWaterBaseProblem>
121{
122 using type = GetPropType<TypeTag, Scalar>;
123 static constexpr type value = 0.75;
124};
125template<class TypeTag>
126struct LensUpperRightZ<TypeTag, TTag::GroundWaterBaseProblem>
127{
128 using type = GetPropType<TypeTag, Scalar>;
129 static constexpr type value = 0.75;
130};
131template<class TypeTag>
132struct Permeability<TypeTag, TTag::GroundWaterBaseProblem>
133{
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 1e-10;
136};
137template<class TypeTag>
138struct PermeabilityLens<TypeTag, TTag::GroundWaterBaseProblem>
139{
140 using type = GetPropType<TypeTag, Scalar>;
141 static constexpr type value = 1e-12;
142};
143
144// Enable gravity
145template<class TypeTag>
146struct EnableGravity<TypeTag, TTag::GroundWaterBaseProblem> { static constexpr bool value = true; };
147
148// The default for the end time of the simulation
149template<class TypeTag>
150struct 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
157template<class TypeTag>
158struct 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
165template<class TypeTag>
166struct 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
170template<class TypeTag>
171struct LinearSolverSplice<TypeTag, TTag::GroundWaterBaseProblem> { using type = TTag::ParallelIstlLinearSolver; };
172
173template<class TypeTag>
174struct LinearSolverWrapper<TypeTag, TTag::GroundWaterBaseProblem>
175{ using type = Opm::Linear::SolverWrapperConjugatedGradients<TypeTag>; };
176
177} // namespace Opm::Properties
178
179namespace Opm {
192template <class TypeTag>
193class 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
226public:
230 GroundWaterProblem(Simulator& simulator)
231 : ParentType(simulator)
232 { }
233
238 {
239 ParentType::finishInit();
240
241 eps_ = 1.0e-3;
242
243 lensLowerLeft_[0] = Parameters::get<TypeTag, Properties::LensLowerLeftX>();
244 if (dim > 1)
245 lensLowerLeft_[1] = Parameters::get<TypeTag, Properties::LensLowerLeftY>();
246 if (dim > 2)
247 lensLowerLeft_[2] = Parameters::get<TypeTag, Properties::LensLowerLeftY>();
248
249 lensUpperRight_[0] = Parameters::get<TypeTag, Properties::LensUpperRightX>();
250 if (dim > 1)
251 lensUpperRight_[1] = Parameters::get<TypeTag, Properties::LensUpperRightY>();
252 if (dim > 2)
253 lensUpperRight_[2] = Parameters::get<TypeTag, Properties::LensUpperRightY>();
254
255 intrinsicPerm_ = this->toDimMatrix_(Parameters::get<TypeTag, Properties::Permeability>());
256 intrinsicPermLens_ = this->toDimMatrix_(Parameters::get<TypeTag, Properties::PermeabilityLens>());
257 }
258
262 static void registerParameters()
263 {
264 ParentType::registerParameters();
265
266 Parameters::registerParam<TypeTag, Properties::LensLowerLeftX>
267 ("The x-coordinate of the lens' lower-left corner [m].");
268 Parameters::registerParam<TypeTag, Properties::LensUpperRightX>
269 ("The x-coordinate of the lens' upper-right corner [m].");
270
271 if (dimWorld > 1) {
272 Parameters::registerParam<TypeTag, Properties::LensLowerLeftY>
273 ("The y-coordinate of the lens' lower-left corner [m].");
274 Parameters::registerParam<TypeTag, Properties::LensUpperRightY>
275 ("The y-coordinate of the lens' upper-right corner [m].");
276 }
277
278 if (dimWorld > 2) {
279 Parameters::registerParam<TypeTag, Properties::LensLowerLeftZ>
280 ("The z-coordinate of the lens' lower-left corner [m].");
281 Parameters::registerParam<TypeTag, Properties::LensUpperRightZ>
282 ("The z-coordinate of the lens' upper-right corner [m].");
283 }
284
285 Parameters::registerParam<TypeTag, Properties::Permeability>
286 ("The intrinsic permeability [m^2] of the ambient material.");
287 Parameters::registerParam<TypeTag, Properties::PermeabilityLens>
288 ("The intrinsic permeability [m^2] of the lens.");
289 }
290
294 // \{
295
299 std::string name() const
300 {
301 std::ostringstream oss;
302 oss << "groundwater_" << Model::name();
303 return oss.str();
304 }
305
310 {
311#ifndef NDEBUG
312 this->model().checkConservativeness();
313
314 // Calculate storage terms
315 EqVector storage;
316 this->model().globalStorage(storage);
317
318 // Write mass balance information for rank 0
319 if (this->gridView().comm().rank() == 0) {
320 std::cout << "Storage: " << storage << std::endl << std::flush;
321 }
322#endif // NDEBUG
323 }
324
328 template <class Context>
329 Scalar temperature(const Context& /*context*/,
330 unsigned /*spaceIdx*/,
331 unsigned /*timeIdx*/) const
332 { return 273.15 + 10; } // 10C
333
337 template <class Context>
338 Scalar porosity(const Context& /*context*/,
339 unsigned /*spaceIdx*/,
340 unsigned /*timeIdx*/) const
341 { return 0.4; }
342
346 template <class Context>
347 const DimMatrix& intrinsicPermeability(const Context& context,
348 unsigned spaceIdx,
349 unsigned timeIdx) const
350 {
351 if (isInLens_(context.pos(spaceIdx, timeIdx)))
352 return intrinsicPermLens_;
353 else
354 return intrinsicPerm_;
355 }
356
358
362
366 template <class Context>
367 void boundary(BoundaryRateVector& values, const Context& context,
368 unsigned spaceIdx, unsigned timeIdx) const
369 {
370 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
371
372 if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos)) {
373 Scalar pressure;
374 Scalar T = temperature(context, spaceIdx, timeIdx);
375 if (onLowerBoundary_(globalPos))
376 pressure = 2e5;
377 else // on upper boundary
378 pressure = 1e5;
379
380 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
381 /*storeEnthalpy=*/false> fs;
382 fs.setSaturation(/*phaseIdx=*/0, 1.0);
383 fs.setPressure(/*phaseIdx=*/0, pressure);
384 fs.setTemperature(T);
385
386 typename FluidSystem::template ParameterCache<Scalar> paramCache;
387 paramCache.updateAll(fs);
388 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
389 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
390 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
391 }
392
393 // impose an freeflow boundary condition
394 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
395 }
396 else {
397 // no flow boundary
398 values.setNoFlow();
399 }
400 }
401
403
408
412 template <class Context>
413 void initial(PrimaryVariables& values,
414 const Context& /*context*/,
415 unsigned /*spaceIdx*/,
416 unsigned /*timeIdx*/) const
417 {
418 // const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
419 values[pressure0Idx] = 1.0e+5; // + 9.81*1.23*(20-globalPos[dim-1]);
420 }
421
425 template <class Context>
426 void source(RateVector& rate,
427 const Context& /*context*/,
428 unsigned /*spaceIdx*/,
429 unsigned /*timeIdx*/) const
430 { rate = Scalar(0.0); }
431
433
434private:
435 bool onLowerBoundary_(const GlobalPosition& pos) const
436 { return pos[dim - 1] < eps_; }
437
438 bool onUpperBoundary_(const GlobalPosition& pos) const
439 { return pos[dim - 1] > this->boundingBoxMax()[dim - 1] - eps_; }
440
441 bool isInLens_(const GlobalPosition& pos) const
442 {
443 return lensLowerLeft_[0] <= pos[0] && pos[0] <= lensUpperRight_[0]
444 && lensLowerLeft_[1] <= pos[1] && pos[1] <= lensUpperRight_[1];
445 }
446
447 GlobalPosition lensLowerLeft_;
448 GlobalPosition lensUpperRight_;
449
450 DimMatrix intrinsicPerm_;
451 DimMatrix intrinsicPermLens_;
452
453 Scalar eps_;
454};
455} // namespace Opm
456
457#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:338
void endTimeStep()
Definition groundwaterproblem.hh:309
GroundWaterProblem(Simulator &simulator)
Definition groundwaterproblem.hh:230
static void registerParameters()
Definition groundwaterproblem.hh:262
std::string name() const
Definition groundwaterproblem.hh:299
Scalar temperature(const Context &, unsigned, unsigned) const
Definition groundwaterproblem.hh:329
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition groundwaterproblem.hh:367
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition groundwaterproblem.hh:426
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition groundwaterproblem.hh:347
void finishInit()
Definition groundwaterproblem.hh:237
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition groundwaterproblem.hh:413
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