28 #ifndef EWOMS_RICHARDS_LENS_PROBLEM_HH
29 #define EWOMS_RICHARDS_LENS_PROBLEM_HH
31 #include <opm/models/richards/richardsmodel.hh>
33 #include <opm/material/components/SimpleH2O.hpp>
34 #include <opm/material/fluidsystems/LiquidPhase.hpp>
35 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
36 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
38 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <dune/grid/yaspgrid.hh>
41 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
43 #include <dune/common/version.hh>
44 #include <dune/common/fvector.hh>
45 #include <dune/common/fmatrix.hh>
48 template <
class TypeTag>
49 class RichardsLensProblem;
53 namespace Opm::Properties {
61 template<
class TypeTag>
65 template<
class TypeTag>
69 template<
class TypeTag>
73 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
76 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
80 template<
class TypeTag>
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>;
103 template<
class TypeTag>
107 template<
class TypeTag>
111 template<
class TypeTag>
115 template<
class TypeTag>
119 template<
class TypeTag>
120 struct NewtonWriteConvergence<TypeTag, TTag::
RichardsLensProblem> {
static constexpr
bool value =
false; };
123 template<
class TypeTag>
126 using type = GetPropType<TypeTag, Scalar>;
127 static constexpr type value = 3000;
131 template<
class TypeTag>
134 using type = GetPropType<TypeTag, Scalar>;
135 static constexpr type value = 100;
139 template<
class TypeTag>
140 struct GridFile<TypeTag, TTag::
RichardsLensProblem> {
static constexpr
auto value =
"./data/richardslens_24x16.dgf"; };
162 template <
class TypeTag>
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_;
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain.
Definition: richardslensproblem.hh:164
void finishInit()
Definition: richardslensproblem.hh:215
std::string name() const
Definition: richardslensproblem.hh:268
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:354
RichardsLensProblem(Simulator &simulator)
Definition: richardslensproblem.hh:205
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: richardslensproblem.hh:453
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:376
void endTimeStep()
Definition: richardslensproblem.hh:279
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:309
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:299
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:429
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: richardslensproblem.hh:323
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:332
Definition: richardslensproblem.hh:57