My Project
Loading...
Searching...
No Matches
fingerproblem.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_FINGER_PROBLEM_HH
29#define EWOMS_FINGER_PROBLEM_HH
30
31#include <opm/models/io/structuredgridvanguard.hh>
32
33#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
34#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
35#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
36#include <opm/material/fluidmatrixinteractions/ParkerLenhard.hpp>
37#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
38
39#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
40#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
41#include <opm/material/components/SimpleH2O.hpp>
42#include <opm/material/components/Air.hpp>
43
44#include <opm/models/immiscible/immiscibleproperties.hh>
45#include <opm/models/discretization/common/restrictprolong.hh>
46
47#if HAVE_DUNE_ALUGRID
48#include <dune/alugrid/grid.hh>
49#endif
50
51#include <dune/common/version.hh>
52#include <dune/common/fvector.hh>
53#include <dune/common/fmatrix.hh>
54#include <dune/grid/utility/persistentcontainer.hh>
55
56#include <vector>
57#include <string>
58
59namespace Opm {
60template <class TypeTag>
61class FingerProblem;
62
63} // namespace Opm
64
65namespace Opm::Properties {
66
67// Create new type tags
68namespace TTag {
69struct FingerBaseProblem { using InheritsFrom = std::tuple<StructuredGridVanguard>; };
70} // end namespace TTag
71
72#if HAVE_DUNE_ALUGRID
73// use dune-alugrid if available
74template<class TypeTag>
75struct Grid<TypeTag, TTag::FingerBaseProblem>
76{ using type = Dune::ALUGrid</*dim=*/2,
77 /*dimWorld=*/2,
78 Dune::cube,
79 Dune::nonconforming>; };
80#endif
81
82// declare the properties used by the finger problem
83template<class TypeTag, class MyTypeTag>
84struct InitialWaterSaturation { using type = UndefinedProperty; };
85
86// Set the problem property
87template<class TypeTag>
88struct Problem<TypeTag, TTag::FingerBaseProblem> { using type = Opm::FingerProblem<TypeTag>; };
89
90// Set the wetting phase
91template<class TypeTag>
92struct WettingPhase<TypeTag, TTag::FingerBaseProblem>
93{
94private:
95 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
96
97public:
98 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
99};
100
101// Set the non-wetting phase
102template<class TypeTag>
103struct NonwettingPhase<TypeTag, TTag::FingerBaseProblem>
104{
105private:
106 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
107
108public:
109 using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
110};
111
112// Set the material Law
113template<class TypeTag>
114struct MaterialLaw<TypeTag, TTag::FingerBaseProblem>
115{
116 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
117 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
118 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
119 /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
120 /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
121
122 // use the parker-lenhard hysteresis law
123 using ParkerLenhard = Opm::ParkerLenhard<Traits>;
124 using type = ParkerLenhard;
125};
126
127// Write the solutions of individual newton iterations?
128template<class TypeTag>
129struct NewtonWriteConvergence<TypeTag, TTag::FingerBaseProblem> { static constexpr bool value = false; };
130
131// Use forward differences instead of central differences
132template<class TypeTag>
133struct NumericDifferenceMethod<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = +1; };
134
135// Enable constraints
136template<class TypeTag>
137struct EnableConstraints<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = true; };
138
139// Enable gravity
140template<class TypeTag>
141struct EnableGravity<TypeTag, TTag::FingerBaseProblem> { static constexpr bool value = true; };
142
143// define the properties specific for the finger problem
144template<class TypeTag>
145struct DomainSizeX<TypeTag, TTag::FingerBaseProblem>
146{
147 using type = GetPropType<TypeTag, Scalar>;
148 static constexpr type value = 0.1;
149};
150template<class TypeTag>
151struct DomainSizeY<TypeTag, TTag::FingerBaseProblem>
152{
153 using type = GetPropType<TypeTag, Scalar>;
154 static constexpr type value = 0.3;
155};
156template<class TypeTag>
157struct DomainSizeZ<TypeTag, TTag::FingerBaseProblem>
158{
159 using type = GetPropType<TypeTag, Scalar>;
160 static constexpr type value = 0.1;
161};
162
163template<class TypeTag>
164struct InitialWaterSaturation<TypeTag, TTag::FingerBaseProblem>
165{
166 using type = GetPropType<TypeTag, Scalar>;
167 static constexpr type value = 0.01;
168};
169
170template<class TypeTag>
171struct CellsX<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 20; };
172template<class TypeTag>
173struct CellsY<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 70; };
174template<class TypeTag>
175struct CellsZ<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 1; };
176
177// The default for the end time of the simulation
178template<class TypeTag>
179struct EndTime<TypeTag, TTag::FingerBaseProblem>
180{
181 using type = GetPropType<TypeTag, Scalar>;
182 static constexpr type value = 215;
183};
184
185// The default for the initial time step size of the simulation
186template<class TypeTag>
187struct InitialTimeStepSize<TypeTag, TTag::FingerBaseProblem>
188{
189 using type = GetPropType<TypeTag, Scalar>;
190 static constexpr type value = 10;
191};
192
193} // namespace Opm::Properties
194
195namespace Opm {
196
212template <class TypeTag>
213class FingerProblem : public GetPropType<TypeTag, Properties::BaseProblem>
214{
216 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
217
218 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
219 using GridView = GetPropType<TypeTag, Properties::GridView>;
220 using Indices = GetPropType<TypeTag, Properties::Indices>;
221 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
222 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
223 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
224 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
225 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
226 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
227 using Model = GetPropType<TypeTag, Properties::Model>;
228
229 enum {
230 // number of phases
231 numPhases = FluidSystem::numPhases,
232
233 // phase indices
234 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
235 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
236
237 // equation indices
238 contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
239
240 // Grid and world dimension
241 dim = GridView::dimension,
242 dimWorld = GridView::dimensionworld
243 };
244
245 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
246 using Stencil = GetPropType<TypeTag, Properties::Stencil> ;
247 enum { codim = Stencil::Entity::codimension };
248 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
249 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
250 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
251
252 using ParkerLenhard = typename GetProp<TypeTag, Properties::MaterialLaw>::ParkerLenhard;
253 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
254 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
255
256 using CoordScalar = typename GridView::ctype;
257 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
258 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
259
260 using Grid = typename GridView :: Grid;
261
262 using MaterialLawParamsContainer = Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > > ;
264
265public:
266 using RestrictProlongOperator = CopyRestrictProlong< Grid, MaterialLawParamsContainer >;
267
271 FingerProblem(Simulator& simulator)
272 : ParentType(simulator),
273 materialParams_( simulator.vanguard().grid(), codim )
274 {
275 }
276
281
285 RestrictProlongOperator restrictProlongOperator()
286 {
287 return RestrictProlongOperator( materialParams_ );
288 }
289
293 std::string name() const
294 { return
295 std::string("finger") +
296 "_" + Model::name() +
297 "_" + Model::discretizationName() +
298 (this->model().enableGridAdaptation()?"_adaptive":"");
299 }
300
304 static void registerParameters()
305 {
306 ParentType::registerParameters();
307
308 EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation,
309 "The initial saturation in the domain [] of the "
310 "wetting phase");
311 }
312
317 {
318 ParentType::finishInit();
319
320 eps_ = 3e-6;
321
322 temperature_ = 273.15 + 20; // -> 20°C
323
324 FluidSystem::init();
325
326 // parameters for the Van Genuchten law of the main imbibition
327 // and the main drainage curves.
328 micParams_.setVgAlpha(0.0037);
329 micParams_.setVgN(4.7);
330 micParams_.finalize();
331
332 mdcParams_.setVgAlpha(0.0037);
333 mdcParams_.setVgN(4.7);
334 mdcParams_.finalize();
335
336 // initialize the material parameter objects of the individual
337 // finite volumes, resize will resize the container to the number of elements
338 materialParams_.resize();
339
340 for (auto it = materialParams_.begin(),
341 end = materialParams_.end(); it != end; ++it ) {
342 std::shared_ptr< MaterialLawParams >& materialParams = *it ;
343 if( ! materialParams )
344 {
345 materialParams.reset( new MaterialLawParams() );
346 materialParams->setMicParams(&micParams_);
347 materialParams->setMdcParams(&mdcParams_);
348 materialParams->setSwr(0.0);
349 materialParams->setSnr(0.1);
350 materialParams->finalize();
351 ParkerLenhard::reset(*materialParams);
352 }
353 }
354
355 K_ = this->toDimMatrix_(4.6e-10);
356
357 setupInitialFluidState_();
358 }
359
364 {
365#ifndef NDEBUG
366 // checkConservativeness() does not include the effect of constraints, so we
367 // disable it for this problem...
368 //this->model().checkConservativeness();
369
370 // Calculate storage terms
371 EqVector storage;
372 this->model().globalStorage(storage);
373
374 // Write mass balance information for rank 0
375 if (this->gridView().comm().rank() == 0) {
376 std::cout << "Storage: " << storage << std::endl << std::flush;
377 }
378#endif // NDEBUG
379
380 // update the history of the hysteresis law
381 ElementContext elemCtx(this->simulator());
382
383 for (const auto& elem : elements(this->gridView())) {
384 elemCtx.updateAll(elem);
385 size_t numDofs = elemCtx.numDof(/*timeIdx=*/0);
386 for (unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
387 {
388 MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, /*timeIdx=*/0 );
389 const auto& fs = elemCtx.intensiveQuantities(scvIdx, /*timeIdx=*/0).fluidState();
390 ParkerLenhard::update(materialParam, fs);
391 }
392 }
393 }
394
396
401
405 template <class Context>
406 Scalar temperature(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
407 { return temperature_; }
408
412 template <class Context>
413 const DimMatrix& intrinsicPermeability(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
414 { return K_; }
415
419 template <class Context>
420 Scalar porosity(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
421 { return 0.4; }
422
426 template <class Context>
427 MaterialLawParams& materialLawParams(const Context& context,
428 unsigned spaceIdx, unsigned timeIdx)
429 {
430 const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
431 assert(materialParams_[entity]);
432 return *materialParams_[entity];
433 }
434
438 template <class Context>
439 const MaterialLawParams& materialLawParams(const Context& context,
440 unsigned spaceIdx, unsigned timeIdx) const
441 {
442 const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
443 assert(materialParams_[entity]);
444 return *materialParams_[entity];
445 }
446
448
453
457 template <class Context>
458 void boundary(BoundaryRateVector& values, const Context& context,
459 unsigned spaceIdx, unsigned timeIdx) const
460 {
461 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
462
463 if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
464 values.setNoFlow();
465 else {
466 assert(onUpperBoundary_(pos));
467
468 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
469 }
470
471 // override the value for the liquid phase by forced
472 // imbibition of water on inlet boundary segments
473 if (onInlet_(pos)) {
474 values[contiWettingEqIdx] = -0.001; // [kg/(m^2 s)]
475 }
476 }
477
479
484
488 template <class Context>
489 void initial(PrimaryVariables& values, const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
490 {
491 // assign the primary variables
492 values.assignNaive(initialFluidState_);
493 }
494
498 template <class Context>
499 void constraints(Constraints& constraints, const Context& context,
500 unsigned spaceIdx, unsigned timeIdx) const
501 {
502 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
503
504 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
505 constraints.setActive(true);
506 constraints.assignNaive(initialFluidState_);
507 }
508 else if (onLowerBoundary_(pos)) {
509 constraints.setActive(true);
510 constraints.assignNaive(initialFluidState_);
511 }
512 }
513
520 template <class Context>
521 void source(RateVector& rate, const Context& /*context*/,
522 unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
523 { rate = Scalar(0.0); }
525
526private:
527 bool onLeftBoundary_(const GlobalPosition& pos) const
528 { return pos[0] < this->boundingBoxMin()[0] + eps_; }
529
530 bool onRightBoundary_(const GlobalPosition& pos) const
531 { return pos[0] > this->boundingBoxMax()[0] - eps_; }
532
533 bool onLowerBoundary_(const GlobalPosition& pos) const
534 { return pos[1] < this->boundingBoxMin()[1] + eps_; }
535
536 bool onUpperBoundary_(const GlobalPosition& pos) const
537 { return pos[1] > this->boundingBoxMax()[1] - eps_; }
538
539 bool onInlet_(const GlobalPosition& pos) const
540 {
541 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
542 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
543
544 if (!onUpperBoundary_(pos))
545 return false;
546
547 Scalar xInject[] = { 0.25, 0.75 };
548 Scalar injectLen[] = { 0.1, 0.1 };
549 for (unsigned i = 0; i < sizeof(xInject) / sizeof(Scalar); ++i) {
550 if (xInject[i] - injectLen[i] / 2 < lambda
551 && lambda < xInject[i] + injectLen[i] / 2)
552 return true;
553 }
554 return false;
555 }
556
557 void setupInitialFluidState_()
558 {
559 auto& fs = initialFluidState_;
560 fs.setPressure(wettingPhaseIdx, /*pressure=*/1e5);
561
562 Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
563 fs.setSaturation(wettingPhaseIdx, Sw);
564 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
565
566 fs.setTemperature(temperature_);
567
568 // set the absolute pressures
569 Scalar pn = 1e5;
570 fs.setPressure(nonWettingPhaseIdx, pn);
571 fs.setPressure(wettingPhaseIdx, pn);
572
573 typename FluidSystem::template ParameterCache<Scalar> paramCache;
574 paramCache.updateAll(fs);
575 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
576 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
577 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
578 }
579
580 }
581
582 DimMatrix K_;
583
584 typename MaterialLawParams::VanGenuchtenParams micParams_;
585 typename MaterialLawParams::VanGenuchtenParams mdcParams_;
586
587 MaterialLawParamsContainer materialParams_;
588
589 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
590
591 Scalar temperature_;
592 Scalar eps_;
593};
594
595} // namespace Opm
596
597#endif
Two-phase problem featuring some gravity-driven saturation fingers.
Definition fingerproblem.hh:214
MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx)
Definition fingerproblem.hh:427
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition fingerproblem.hh:458
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition fingerproblem.hh:439
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition fingerproblem.hh:413
Scalar porosity(const Context &, unsigned, unsigned) const
Definition fingerproblem.hh:420
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition fingerproblem.hh:521
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition fingerproblem.hh:499
void finishInit()
Definition fingerproblem.hh:316
Scalar temperature(const Context &, unsigned, unsigned) const
Definition fingerproblem.hh:406
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition fingerproblem.hh:489
void endTimeStep()
Definition fingerproblem.hh:363
RestrictProlongOperator restrictProlongOperator()
Definition fingerproblem.hh:285
static void registerParameters()
Definition fingerproblem.hh:304
std::string name() const
Definition fingerproblem.hh:293
FingerProblem(Simulator &simulator)
Definition fingerproblem.hh:271
Definition fingerproblem.hh:84
Definition fingerproblem.hh:69