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 Parameters::registerParam<TypeTag, Properties::InitialWaterSaturation>
309 (
"The initial saturation in the domain [] of the wetting phase");
317 ParentType::finishInit();
321 temperature_ = 273.15 + 20;
327 micParams_.setVgAlpha(0.0037);
328 micParams_.setVgN(4.7);
329 micParams_.finalize();
331 mdcParams_.setVgAlpha(0.0037);
332 mdcParams_.setVgN(4.7);
333 mdcParams_.finalize();
337 materialParams_.resize();
339 for (
auto it = materialParams_.begin(),
340 end = materialParams_.end(); it != end; ++it ) {
341 std::shared_ptr< MaterialLawParams >& materialParams = *it ;
342 if( ! materialParams )
344 materialParams.reset(
new MaterialLawParams() );
345 materialParams->setMicParams(&micParams_);
346 materialParams->setMdcParams(&mdcParams_);
347 materialParams->setSwr(0.0);
348 materialParams->setSnr(0.1);
349 materialParams->finalize();
350 ParkerLenhard::reset(*materialParams);
354 K_ = this->toDimMatrix_(4.6e-10);
356 setupInitialFluidState_();
371 this->model().globalStorage(storage);
374 if (this->gridView().comm().rank() == 0) {
375 std::cout <<
"Storage: " << storage << std::endl << std::flush;
380 ElementContext elemCtx(this->simulator());
382 for (
const auto& elem : elements(this->gridView())) {
383 elemCtx.updateAll(elem);
384 size_t numDofs = elemCtx.numDof(0);
385 for (
unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
388 const auto& fs = elemCtx.intensiveQuantities(scvIdx, 0).fluidState();
389 ParkerLenhard::update(materialParam, fs);
404 template <
class Context>
406 {
return temperature_; }
411 template <
class Context>
418 template <
class Context>
419 Scalar
porosity(
const Context& ,
unsigned ,
unsigned )
const
425 template <
class Context>
427 unsigned spaceIdx,
unsigned timeIdx)
429 const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
430 assert(materialParams_[entity]);
431 return *materialParams_[entity];
437 template <
class Context>
439 unsigned spaceIdx,
unsigned timeIdx)
const
441 const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
442 assert(materialParams_[entity]);
443 return *materialParams_[entity];
456 template <
class Context>
457 void boundary(BoundaryRateVector& values,
const Context& context,
458 unsigned spaceIdx,
unsigned timeIdx)
const
460 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
462 if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
465 assert(onUpperBoundary_(pos));
467 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
473 values[contiWettingEqIdx] = -0.001;
487 template <
class Context>
488 void initial(PrimaryVariables& values,
const Context& ,
unsigned ,
unsigned )
const
491 values.assignNaive(initialFluidState_);
497 template <
class Context>
499 unsigned spaceIdx,
unsigned timeIdx)
const
501 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
503 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
507 else if (onLowerBoundary_(pos)) {
519 template <
class Context>
520 void source(RateVector& rate,
const Context& ,
521 unsigned ,
unsigned )
const
522 { rate = Scalar(0.0); }
526 bool onLeftBoundary_(
const GlobalPosition& pos)
const
527 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
529 bool onRightBoundary_(
const GlobalPosition& pos)
const
530 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
532 bool onLowerBoundary_(
const GlobalPosition& pos)
const
533 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
535 bool onUpperBoundary_(
const GlobalPosition& pos)
const
536 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
538 bool onInlet_(
const GlobalPosition& pos)
const
540 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
541 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
543 if (!onUpperBoundary_(pos))
546 Scalar xInject[] = { 0.25, 0.75 };
547 Scalar injectLen[] = { 0.1, 0.1 };
548 for (
unsigned i = 0; i <
sizeof(xInject) /
sizeof(Scalar); ++i) {
549 if (xInject[i] - injectLen[i] / 2 < lambda
550 && lambda < xInject[i] + injectLen[i] / 2)
556 void setupInitialFluidState_()
558 auto& fs = initialFluidState_;
559 fs.setPressure(wettingPhaseIdx, 1e5);
561 Scalar Sw = Parameters::get<TypeTag, Properties::InitialWaterSaturation>();
562 fs.setSaturation(wettingPhaseIdx, Sw);
563 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
565 fs.setTemperature(temperature_);
569 fs.setPressure(nonWettingPhaseIdx, pn);
570 fs.setPressure(wettingPhaseIdx, pn);
572 typename FluidSystem::template ParameterCache<Scalar> paramCache;
573 paramCache.updateAll(fs);
574 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
575 fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
576 fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
583 typename MaterialLawParams::VanGenuchtenParams micParams_;
584 typename MaterialLawParams::VanGenuchtenParams mdcParams_;
586 MaterialLawParamsContainer materialParams_;
588 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;