My Project
Loading...
Searching...
No Matches
infiltrationproblem.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*/
27#ifndef EWOMS_INFILTRATION_PROBLEM_HH
28#define EWOMS_INFILTRATION_PROBLEM_HH
29
30#include <opm/models/pvs/pvsproperties.hh>
31
32#include <opm/material/fluidstates/CompositionalFluidState.hpp>
33#include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>
34#include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
35#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
36#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
37#include <opm/material/common/Valgrind.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/fvector.hh>
44#include <dune/common/fmatrix.hh>
45
46#include <sstream>
47#include <string>
48
49namespace Opm {
50template <class TypeTag>
51class InfiltrationProblem;
52}
53
54namespace Opm::Properties {
55
56namespace TTag {
58}
59
60// Set the grid type
61template<class TypeTag>
62struct Grid<TypeTag, TTag::InfiltrationBaseProblem> { using type = Dune::YaspGrid<2>; };
63
64// Set the problem property
65template<class TypeTag>
66struct Problem<TypeTag, TTag::InfiltrationBaseProblem> { using type = Opm::InfiltrationProblem<TypeTag>; };
67
68// Set the fluid system
69template<class TypeTag>
70struct FluidSystem<TypeTag, TTag::InfiltrationBaseProblem>
71{ using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
72
73// Enable gravity?
74template<class TypeTag>
75struct EnableGravity<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr bool value = true; };
76
77// Write newton convergence?
78template<class TypeTag>
79struct NewtonWriteConvergence<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr bool value = false; };
80
81// -1 backward differences, 0: central differences, +1: forward differences
82template<class TypeTag>
83struct NumericDifferenceMethod<TypeTag, TTag::InfiltrationBaseProblem> { static constexpr int value = 1; };
84
85// Set the material Law
86template<class TypeTag>
87struct MaterialLaw<TypeTag, TTag::InfiltrationBaseProblem>
88{
89private:
90 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
92
93 using Traits= Opm::ThreePhaseMaterialTraits<
94 Scalar,
95 /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
96 /*nonWettingPhaseIdx=*/FluidSystem::naplPhaseIdx,
97 /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
98
99public:
100 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
101};
102
103// The default for the end time of the simulation
104template<class TypeTag>
105struct EndTime<TypeTag, TTag::InfiltrationBaseProblem>
106{
107 using type = GetPropType<TypeTag, Scalar>;
108 static constexpr type value = 6e3;
109};
110
111// The default for the initial time step size of the simulation
112template<class TypeTag>
113struct InitialTimeStepSize<TypeTag, TTag::InfiltrationBaseProblem>
114{
115 using type = GetPropType<TypeTag, Scalar>;
116 static constexpr type value = 60;
117};
118
119// The default DGF file to load
120template<class TypeTag>
121struct GridFile<TypeTag, TTag::InfiltrationBaseProblem>
122{ static constexpr auto value = "./data/infiltration_50x3.dgf"; };
123
124} // namespace Opm::Properties
125
126namespace Opm {
149template <class TypeTag>
150class InfiltrationProblem : public GetPropType<TypeTag, Properties::BaseProblem>
151{
152 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
153
154 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
155 using GridView = GetPropType<TypeTag, Properties::GridView>;
156 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
157 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
158 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
159 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
160 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
161 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
162 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
163 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
164 using Model = GetPropType<TypeTag, Properties::Model>;
165
166 // copy some indices for convenience
167 using Indices = GetPropType<TypeTag, Properties::Indices>;
168 enum {
169 // equation indices
170 conti0EqIdx = Indices::conti0EqIdx,
171
172 // number of phases/components
173 numPhases = FluidSystem::numPhases,
174
175 // component indices
176 NAPLIdx = FluidSystem::NAPLIdx,
177 H2OIdx = FluidSystem::H2OIdx,
178 airIdx = FluidSystem::airIdx,
179
180 // phase indices
181 waterPhaseIdx = FluidSystem::waterPhaseIdx,
182 gasPhaseIdx = FluidSystem::gasPhaseIdx,
183 naplPhaseIdx = FluidSystem::naplPhaseIdx,
184
185 // Grid and world dimension
186 dim = GridView::dimension,
187 dimWorld = GridView::dimensionworld
188 };
189
190 using CoordScalar = typename GridView::ctype;
191 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
192 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
193
194public:
198 InfiltrationProblem(Simulator& simulator)
199 : ParentType(simulator)
200 , eps_(1e-6)
201 { }
202
207 {
208 ParentType::finishInit();
209
210 temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius
211 FluidSystem::init(/*tempMin=*/temperature_ - 1,
212 /*tempMax=*/temperature_ + 1,
213 /*nTemp=*/3,
214 /*pressMin=*/0.8 * 1e5,
215 /*pressMax=*/3 * 1e5,
216 /*nPress=*/200);
217
218 // intrinsic permeabilities
219 fineK_ = this->toDimMatrix_(1e-11);
220 coarseK_ = this->toDimMatrix_(1e-11);
221
222 // porosities
223 porosity_ = 0.40;
224
225 // residual saturations
226 materialParams_.setSwr(0.12);
227 materialParams_.setSwrx(0.12);
228 materialParams_.setSnr(0.07);
229 materialParams_.setSgr(0.03);
230
231 // parameters for the three-phase van Genuchten law
232 materialParams_.setVgAlpha(0.0005);
233 materialParams_.setVgN(4.);
234 materialParams_.setkrRegardsSnr(false);
235
236 materialParams_.finalize();
237 materialParams_.checkDefined();
238 }
239
244
251 { return true; }
252
256 std::string name() const
257 {
258 std::ostringstream oss;
259 oss << "infiltration_" << Model::name();
260 return oss.str();
261 }
262
267 {
268#ifndef NDEBUG
269 this->model().checkConservativeness();
270
271 // Calculate storage terms
272 EqVector storage;
273 this->model().globalStorage(storage);
274
275 // Write mass balance information for rank 0
276 if (this->gridView().comm().rank() == 0) {
277 std::cout << "Storage: " << storage << std::endl << std::flush;
278 }
279#endif // NDEBUG
280 }
281
285 template <class Context>
286 Scalar temperature(const Context& /*context*/,
287 unsigned /*spaceIdx*/,
288 unsigned /*timeIdx*/) const
289 { return temperature_; }
290
294 template <class Context>
295 const DimMatrix&
296 intrinsicPermeability(const Context& context,
297 unsigned spaceIdx,
298 unsigned timeIdx) const
299 {
300 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
301 if (isFineMaterial_(pos))
302 return fineK_;
303 return coarseK_;
304 }
305
309 template <class Context>
310 Scalar porosity(const Context& /*context*/,
311 unsigned /*spaceIdx*/,
312 unsigned /*timeIdx*/) const
313 { return porosity_; }
314
318 template <class Context>
319 const MaterialLawParams&
320 materialLawParams(const Context& /*context*/,
321 unsigned /*spaceIdx*/,
322 unsigned /*timeIdx*/) const
323 { return materialParams_; }
324
326
331
335 template <class Context>
336 void boundary(BoundaryRateVector& values,
337 const Context& context,
338 unsigned spaceIdx,
339 unsigned timeIdx) const
340 {
341 const auto& pos = context.pos(spaceIdx, timeIdx);
342
343 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
344 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
345
346 initialFluidState_(fs, context, spaceIdx, timeIdx);
347
348 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
349 }
350 else if (onInlet_(pos)) {
351 RateVector molarRate(0.0);
352 molarRate[conti0EqIdx + NAPLIdx] = -0.001;
353
354 values.setMolarRate(molarRate);
355 Opm::Valgrind::CheckDefined(values);
356 }
357 else
358 values.setNoFlow();
359 }
360
362
367
371 template <class Context>
372 void initial(PrimaryVariables& values,
373 const Context& context,
374 unsigned spaceIdx,
375 unsigned timeIdx) const
376 {
377 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
378
379 initialFluidState_(fs, context, spaceIdx, timeIdx);
380
381 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
382 values.assignMassConservative(fs, matParams, /*inEquilibrium=*/true);
383 Opm::Valgrind::CheckDefined(values);
384 }
385
392 template <class Context>
393 void source(RateVector& rate,
394 const Context& /*context*/,
395 unsigned /*spaceIdx*/,
396 unsigned /*timeIdx*/) const
397 { rate = Scalar(0.0); }
398
400
401private:
402 bool onLeftBoundary_(const GlobalPosition& pos) const
403 { return pos[0] < eps_; }
404
405 bool onRightBoundary_(const GlobalPosition& pos) const
406 { return pos[0] > this->boundingBoxMax()[0] - eps_; }
407
408 bool onLowerBoundary_(const GlobalPosition& pos) const
409 { return pos[1] < eps_; }
410
411 bool onUpperBoundary_(const GlobalPosition& pos) const
412 { return pos[1] > this->boundingBoxMax()[1] - eps_; }
413
414 bool onInlet_(const GlobalPosition& pos) const
415 { return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
416
417 template <class FluidState, class Context>
418 void initialFluidState_(FluidState& fs, const Context& context,
419 unsigned spaceIdx, unsigned timeIdx) const
420 {
421 const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
422 Scalar y = pos[1];
423 Scalar x = pos[0];
424
425 Scalar densityW = 1000.0;
426 Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
427 if (pc < 0.0)
428 pc = 0.0;
429
430 // set pressures
431 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
432 Scalar Sw = matParams.Swr();
433 Scalar Swr = matParams.Swr();
434 Scalar Sgr = matParams.Sgr();
435 if (Sw < Swr)
436 Sw = Swr;
437 if (Sw > 1 - Sgr)
438 Sw = 1 - Sgr;
439 Scalar Sg = 1 - Sw;
440
441 Opm::Valgrind::CheckDefined(Sw);
442 Opm::Valgrind::CheckDefined(Sg);
443
444 fs.setSaturation(waterPhaseIdx, Sw);
445 fs.setSaturation(gasPhaseIdx, Sg);
446 fs.setSaturation(naplPhaseIdx, 0);
447
448 // set temperature of all phases
449 fs.setTemperature(temperature_);
450
451 // compute pressures
452 Scalar pcAll[numPhases];
453 Scalar pg = 1e5;
454 if (onLeftBoundary_(pos))
455 pg += 10e3;
456 MaterialLaw::capillaryPressures(pcAll, matParams, fs);
457 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
458 fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
459
460 // set composition of gas phase
461 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
462 fs.setMoleFraction(gasPhaseIdx, airIdx,
463 1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
464 fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
465
466 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
467 typename FluidSystem::template ParameterCache<Scalar> paramCache;
468 CFRP::solve(fs, paramCache, gasPhaseIdx,
469 /*setViscosity=*/true,
470 /*setEnthalpy=*/false);
471
472 fs.setMoleFraction(waterPhaseIdx, H2OIdx,
473 1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
474 }
475
476 bool isFineMaterial_(const GlobalPosition& pos) const
477 { return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
478
479 DimMatrix fineK_;
480 DimMatrix coarseK_;
481
482 Scalar porosity_;
483
484 MaterialLawParams materialParams_;
485
486 Scalar temperature_;
487 Scalar eps_;
488};
489} // namespace Opm
490
491#endif
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition infiltrationproblem.hh:151
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition infiltrationproblem.hh:393
Scalar temperature(const Context &, unsigned, unsigned) const
Definition infiltrationproblem.hh:286
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition infiltrationproblem.hh:296
std::string name() const
Definition infiltrationproblem.hh:256
Scalar porosity(const Context &, unsigned, unsigned) const
Definition infiltrationproblem.hh:310
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition infiltrationproblem.hh:336
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition infiltrationproblem.hh:320
void endTimeStep()
Definition infiltrationproblem.hh:266
InfiltrationProblem(Simulator &simulator)
Definition infiltrationproblem.hh:198
bool shouldWriteRestartFile() const
Definition infiltrationproblem.hh:250
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition infiltrationproblem.hh:372
void finishInit()
Definition infiltrationproblem.hh:206
Definition infiltrationproblem.hh:57