28 #ifndef EWOMS_FINGER_PROBLEM_HH
29 #define EWOMS_FINGER_PROBLEM_HH
31 #include <opm/models/io/structuredgridvanguard.hh>
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>
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>
44 #include <opm/models/immiscible/immiscibleproperties.hh>
45 #include <opm/models/discretization/common/restrictprolong.hh>
48 #include <dune/alugrid/grid.hh>
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>
60 template <
class TypeTag>
65 namespace Opm::Properties {
74 template<
class TypeTag>
75 struct Grid<TypeTag, TTag::FingerBaseProblem>
76 {
using type = Dune::ALUGrid<2,
79 Dune::nonconforming>; };
83 template<
class TypeTag,
class MyTypeTag>
87 template<
class TypeTag>
91 template<
class TypeTag>
92 struct WettingPhase<TypeTag, TTag::FingerBaseProblem>
95 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
98 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
102 template<
class TypeTag>
103 struct NonwettingPhase<TypeTag, TTag::FingerBaseProblem>
106 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
109 using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
113 template<
class TypeTag>
114 struct MaterialLaw<TypeTag, TTag::FingerBaseProblem>
116 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
117 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
118 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
119 FluidSystem::wettingPhaseIdx,
120 FluidSystem::nonWettingPhaseIdx>;
123 using ParkerLenhard = Opm::ParkerLenhard<Traits>;
124 using type = ParkerLenhard;
128 template<
class TypeTag>
129 struct NewtonWriteConvergence<TypeTag, TTag::FingerBaseProblem> {
static constexpr
bool value =
false; };
132 template<
class TypeTag>
133 struct NumericDifferenceMethod<TypeTag, TTag::FingerBaseProblem> {
static constexpr
int value = +1; };
136 template<
class TypeTag>
137 struct EnableConstraints<TypeTag, TTag::FingerBaseProblem> {
static constexpr
int value =
true; };
140 template<
class TypeTag>
141 struct EnableGravity<TypeTag, TTag::FingerBaseProblem> {
static constexpr
bool value =
true; };
144 template<
class TypeTag>
145 struct DomainSizeX<TypeTag, TTag::FingerBaseProblem>
147 using type = GetPropType<TypeTag, Scalar>;
148 static constexpr type value = 0.1;
150 template<
class TypeTag>
151 struct DomainSizeY<TypeTag, TTag::FingerBaseProblem>
153 using type = GetPropType<TypeTag, Scalar>;
154 static constexpr type value = 0.3;
156 template<
class TypeTag>
157 struct DomainSizeZ<TypeTag, TTag::FingerBaseProblem>
159 using type = GetPropType<TypeTag, Scalar>;
160 static constexpr type value = 0.1;
163 template<
class TypeTag>
166 using type = GetPropType<TypeTag, Scalar>;
167 static constexpr type value = 0.01;
170 template<
class TypeTag>
171 struct CellsX<TypeTag, TTag::FingerBaseProblem> {
static constexpr
int value = 20; };
172 template<
class TypeTag>
173 struct CellsY<TypeTag, TTag::FingerBaseProblem> {
static constexpr
int value = 70; };
174 template<
class TypeTag>
175 struct CellsZ<TypeTag, TTag::FingerBaseProblem> {
static constexpr
int value = 1; };
178 template<
class TypeTag>
179 struct EndTime<TypeTag, TTag::FingerBaseProblem>
181 using type = GetPropType<TypeTag, Scalar>;
182 static constexpr type value = 215;
186 template<
class TypeTag>
187 struct InitialTimeStepSize<TypeTag, TTag::FingerBaseProblem>
189 using type = GetPropType<TypeTag, Scalar>;
190 static constexpr type value = 10;
212 template <
class TypeTag>
216 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
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>;
231 numPhases = FluidSystem::numPhases,
234 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
235 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
238 contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
241 dim = GridView::dimension,
242 dimWorld = GridView::dimensionworld
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>;
252 using ParkerLenhard =
typename GetProp<TypeTag, Properties::MaterialLaw>::ParkerLenhard;
253 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
254 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
256 using CoordScalar =
typename GridView::ctype;
257 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
258 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
260 using Grid =
typename GridView :: Grid;
262 using MaterialLawParamsContainer = Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > > ;
266 using RestrictProlongOperator = CopyRestrictProlong< Grid, MaterialLawParamsContainer >;
272 : ParentType(simulator),
273 materialParams_( simulator.vanguard().grid(), codim )
287 return RestrictProlongOperator( materialParams_ );
295 std::string(
"finger") +
296 "_" + Model::name() +
297 "_" + Model::discretizationName() +
298 (this->model().enableGridAdaptation()?
"_adaptive":
"");
306 ParentType::registerParameters();
308 EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation,
309 "The initial saturation in the domain [] of the "
318 ParentType::finishInit();
322 temperature_ = 273.15 + 20;
328 micParams_.setVgAlpha(0.0037);
329 micParams_.setVgN(4.7);
330 micParams_.finalize();
332 mdcParams_.setVgAlpha(0.0037);
333 mdcParams_.setVgN(4.7);
334 mdcParams_.finalize();
338 materialParams_.resize();
340 for (
auto it = materialParams_.begin(),
341 end = materialParams_.end(); it != end; ++it ) {
342 std::shared_ptr< MaterialLawParams >& materialParams = *it ;
343 if( ! materialParams )
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);
355 K_ = this->toDimMatrix_(4.6e-10);
357 setupInitialFluidState_();
372 this->model().globalStorage(storage);
375 if (this->gridView().comm().rank() == 0) {
376 std::cout <<
"Storage: " << storage << std::endl << std::flush;
381 ElementContext elemCtx(this->simulator());
383 auto elemIt = this->gridView().template begin<0>();
384 const auto& elemEndIt = this->gridView().template end<0>();
385 for (; elemIt != elemEndIt; ++elemIt) {
386 const auto& elem = *elemIt;
387 elemCtx.updateAll( elem );
388 size_t numDofs = elemCtx.numDof(0);
389 for (
unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
392 const auto& fs = elemCtx.intensiveQuantities(scvIdx, 0).fluidState();
393 ParkerLenhard::update(materialParam, fs);
408 template <
class Context>
410 {
return temperature_; }
415 template <
class Context>
422 template <
class Context>
423 Scalar
porosity(
const Context& ,
unsigned ,
unsigned )
const
429 template <
class Context>
431 unsigned spaceIdx,
unsigned timeIdx)
433 const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
434 assert(materialParams_[entity]);
435 return *materialParams_[entity];
441 template <
class Context>
443 unsigned spaceIdx,
unsigned timeIdx)
const
445 const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
446 assert(materialParams_[entity]);
447 return *materialParams_[entity];
460 template <
class Context>
461 void boundary(BoundaryRateVector& values,
const Context& context,
462 unsigned spaceIdx,
unsigned timeIdx)
const
464 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
466 if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
469 assert(onUpperBoundary_(pos));
471 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
477 values[contiWettingEqIdx] = -0.001;
491 template <
class Context>
492 void initial(PrimaryVariables& values,
const Context& ,
unsigned ,
unsigned )
const
495 values.assignNaive(initialFluidState_);
501 template <
class Context>
503 unsigned spaceIdx,
unsigned timeIdx)
const
505 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
507 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
511 else if (onLowerBoundary_(pos)) {
523 template <
class Context>
524 void source(RateVector& rate,
const Context& ,
525 unsigned ,
unsigned )
const
526 { rate = Scalar(0.0); }
530 bool onLeftBoundary_(
const GlobalPosition& pos)
const
531 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
533 bool onRightBoundary_(
const GlobalPosition& pos)
const
534 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
536 bool onLowerBoundary_(
const GlobalPosition& pos)
const
537 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
539 bool onUpperBoundary_(
const GlobalPosition& pos)
const
540 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
542 bool onInlet_(
const GlobalPosition& pos)
const
544 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
545 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
547 if (!onUpperBoundary_(pos))
550 Scalar xInject[] = { 0.25, 0.75 };
551 Scalar injectLen[] = { 0.1, 0.1 };
552 for (
unsigned i = 0; i <
sizeof(xInject) /
sizeof(Scalar); ++i) {
553 if (xInject[i] - injectLen[i] / 2 < lambda
554 && lambda < xInject[i] + injectLen[i] / 2)
560 void setupInitialFluidState_()
562 auto& fs = initialFluidState_;
563 fs.setPressure(wettingPhaseIdx, 1e5);
565 Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
566 fs.setSaturation(wettingPhaseIdx, Sw);
567 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
569 fs.setTemperature(temperature_);
573 fs.setPressure(nonWettingPhaseIdx, pn);
574 fs.setPressure(wettingPhaseIdx, pn);
576 typename FluidSystem::template ParameterCache<Scalar> paramCache;
577 paramCache.updateAll(fs);
578 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
579 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
580 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
587 typename MaterialLawParams::VanGenuchtenParams micParams_;
588 typename MaterialLawParams::VanGenuchtenParams mdcParams_;
590 MaterialLawParamsContainer materialParams_;
592 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
Two-phase problem featuring some gravity-driven saturation fingers.
Definition: fingerproblem.hh:214
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:416
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:461
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:423
MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx)
Definition: fingerproblem.hh:430
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:524
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:502
void finishInit()
Definition: fingerproblem.hh:316
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:409
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:492
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
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:442
Definition: fingerproblem.hh:84
Definition: fingerproblem.hh:69