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 for (
const auto& elem : elements(this->gridView())) {
384 elemCtx.updateAll(elem);
385 size_t numDofs = elemCtx.numDof(0);
386 for (
unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
389 const auto& fs = elemCtx.intensiveQuantities(scvIdx, 0).fluidState();
390 ParkerLenhard::update(materialParam, fs);
405 template <
class Context>
407 {
return temperature_; }
412 template <
class Context>
419 template <
class Context>
420 Scalar
porosity(
const Context& ,
unsigned ,
unsigned )
const
426 template <
class Context>
428 unsigned spaceIdx,
unsigned timeIdx)
430 const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
431 assert(materialParams_[entity]);
432 return *materialParams_[entity];
438 template <
class Context>
440 unsigned spaceIdx,
unsigned timeIdx)
const
442 const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
443 assert(materialParams_[entity]);
444 return *materialParams_[entity];
457 template <
class Context>
458 void boundary(BoundaryRateVector& values,
const Context& context,
459 unsigned spaceIdx,
unsigned timeIdx)
const
461 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
463 if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
466 assert(onUpperBoundary_(pos));
468 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
474 values[contiWettingEqIdx] = -0.001;
488 template <
class Context>
489 void initial(PrimaryVariables& values,
const Context& ,
unsigned ,
unsigned )
const
492 values.assignNaive(initialFluidState_);
498 template <
class Context>
500 unsigned spaceIdx,
unsigned timeIdx)
const
502 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
504 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
508 else if (onLowerBoundary_(pos)) {
520 template <
class Context>
521 void source(RateVector& rate,
const Context& ,
522 unsigned ,
unsigned )
const
523 { rate = Scalar(0.0); }
527 bool onLeftBoundary_(
const GlobalPosition& pos)
const
528 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
530 bool onRightBoundary_(
const GlobalPosition& pos)
const
531 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
533 bool onLowerBoundary_(
const GlobalPosition& pos)
const
534 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
536 bool onUpperBoundary_(
const GlobalPosition& pos)
const
537 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
539 bool onInlet_(
const GlobalPosition& pos)
const
541 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
542 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
544 if (!onUpperBoundary_(pos))
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)
557 void setupInitialFluidState_()
559 auto& fs = initialFluidState_;
560 fs.setPressure(wettingPhaseIdx, 1e5);
562 Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
563 fs.setSaturation(wettingPhaseIdx, Sw);
564 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
566 fs.setTemperature(temperature_);
570 fs.setPressure(nonWettingPhaseIdx, pn);
571 fs.setPressure(wettingPhaseIdx, pn);
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));
584 typename MaterialLawParams::VanGenuchtenParams micParams_;
585 typename MaterialLawParams::VanGenuchtenParams mdcParams_;
587 MaterialLawParamsContainer materialParams_;
589 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;