84 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
85 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
86 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
88 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
89 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
90 FluidSystem::wettingPhaseIdx,
91 FluidSystem::nonWettingPhaseIdx>;
95 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
99 using type = Opm::EffToAbsLaw<EffectiveLaw>;
165 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
167 using GridView = GetPropType<TypeTag, Properties::GridView>;
168 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
169 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
170 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
171 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
172 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
173 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
174 using Model = GetPropType<TypeTag, Properties::Model>;
175 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
176 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
178 using Indices = GetPropType<TypeTag, Properties::Indices>;
181 pressureWIdx = Indices::pressureWIdx,
182 contiEqIdx = Indices::contiEqIdx,
183 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
184 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
185 numPhases = FluidSystem::numPhases,
188 dimWorld = GridView::dimensionworld
192 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
194 using MaterialLawParams =
typename MaterialLaw::Params;
196 using CoordScalar =
typename GridView::ctype;
197 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
198 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
199 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
206 : ParentType(simulator)
209 dofIsInLens_.resize(simulator.model().numGridDof());
217 ParentType::finishInit();
222 lensLowerLeft_[0] = 1.0;
223 lensLowerLeft_[1] = 2.0;
225 lensUpperRight_[0] = 4.0;
226 lensUpperRight_[1] = 3.0;
230 lensMaterialParams_.setVgAlpha(0.00045);
231 lensMaterialParams_.setVgN(7.3);
232 lensMaterialParams_.finalize();
234 outerMaterialParams_.setVgAlpha(0.0037);
235 outerMaterialParams_.setVgN(4.7);
236 outerMaterialParams_.finalize();
245 lensK_ = this->toDimMatrix_(1e-12);
246 outerK_ = this->toDimMatrix_(5e-12);
249 Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
250 for (
const auto& elem : elements(this->gridView())) {
251 stencil.update(elem);
252 for (
unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
253 unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
254 const auto& dofPos = stencil.subControlVolume(dofIdx).center();
255 dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
270 std::ostringstream oss;
271 oss <<
"lens_richards_"
272 << Model::discretizationName();
282 this->model().checkConservativeness();
286 this->model().globalStorage(storage);
289 if (this->gridView().comm().rank() == 0) {
290 std::cout <<
"Storage: " << storage << std::endl << std::flush;
298 template <
class Context>
299 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
300 {
return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
303 {
return 273.15 + 10; }
308 template <
class Context>
311 unsigned timeIdx)
const
313 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
322 template <
class Context>
331 template <
class Context>
334 unsigned timeIdx)
const
336 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
343 if (dofIsInLens_[globalSpaceIdx])
344 return lensMaterialParams_;
345 return outerMaterialParams_;
353 template <
class Context>
356 unsigned timeIdx)
const
357 {
return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
375 template <
class Context>
377 const Context& context,
379 unsigned timeIdx)
const
381 const auto& pos = context.pos(spaceIdx, timeIdx);
383 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
384 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
387 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
388 fs.setSaturation(wettingPhaseIdx, Sw);
389 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
392 MaterialLaw::capillaryPressures(pC, materialParams, fs);
393 fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
394 fs.setPressure(nonWettingPhaseIdx, pnRef_);
396 typename FluidSystem::template ParameterCache<Scalar> paramCache;
397 paramCache.updateAll(fs);
398 fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx));
401 fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx));
404 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
406 else if (onInlet_(pos)) {
407 RateVector massRate(0.0);
410 massRate[contiEqIdx] = -0.04;
412 values.setMassRate(massRate);
428 template <
class Context>
430 const Context& context,
432 unsigned timeIdx)
const
434 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
437 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
438 fs.setSaturation(wettingPhaseIdx, Sw);
439 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
442 MaterialLaw::capillaryPressures(pC, materialParams, fs);
443 values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
452 template <
class Context>
457 { rate = Scalar(0.0); }
462 bool onLeftBoundary_(
const GlobalPosition& pos)
const
463 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
465 bool onRightBoundary_(
const GlobalPosition& pos)
const
466 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
468 bool onLowerBoundary_(
const GlobalPosition& pos)
const
469 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
471 bool onUpperBoundary_(
const GlobalPosition& pos)
const
472 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
474 bool onInlet_(
const GlobalPosition& pos)
const
476 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
477 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
478 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
481 bool isInLens_(
const GlobalPosition& pos)
const
483 for (
unsigned i = 0; i < dimWorld; ++i) {
484 if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
490 GlobalPosition lensLowerLeft_;
491 GlobalPosition lensUpperRight_;
495 MaterialLawParams lensMaterialParams_;
496 MaterialLawParams outerMaterialParams_;
498 std::vector<bool> dofIsInLens_;