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>
39 #include <opm/material/common/Unused.hpp>
41 #include <dune/grid/yaspgrid.hh>
42 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
44 #include <dune/common/version.hh>
45 #include <dune/common/fvector.hh>
46 #include <dune/common/fmatrix.hh>
49 template <
class TypeTag>
50 class RichardsLensProblem;
54 namespace Opm::Properties {
62 template<
class TypeTag>
66 template<
class TypeTag>
70 template<
class TypeTag>
74 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
77 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
81 template<
class TypeTag>
85 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
86 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
87 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
89 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
90 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
91 FluidSystem::wettingPhaseIdx,
92 FluidSystem::nonWettingPhaseIdx>;
96 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
100 using type = Opm::EffToAbsLaw<EffectiveLaw>;
104 template<
class TypeTag>
108 template<
class TypeTag>
112 template<
class TypeTag>
116 template<
class TypeTag>
120 template<
class TypeTag>
121 struct NewtonWriteConvergence<TypeTag, TTag::
RichardsLensProblem> {
static constexpr
bool value =
false; };
124 template<
class TypeTag>
127 using type = GetPropType<TypeTag, Scalar>;
128 static constexpr type value = 3000;
132 template<
class TypeTag>
135 using type = GetPropType<TypeTag, Scalar>;
136 static constexpr type value = 100;
140 template<
class TypeTag>
141 struct GridFile<TypeTag, TTag::
RichardsLensProblem> {
static constexpr
auto value =
"./data/richardslens_24x16.dgf"; };
163 template <
class TypeTag>
166 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
168 using GridView = GetPropType<TypeTag, Properties::GridView>;
169 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
170 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
171 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
172 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
173 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
174 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
175 using Model = GetPropType<TypeTag, Properties::Model>;
176 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
177 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
179 using Indices = GetPropType<TypeTag, Properties::Indices>;
182 pressureWIdx = Indices::pressureWIdx,
183 contiEqIdx = Indices::contiEqIdx,
184 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
185 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
186 numPhases = FluidSystem::numPhases,
189 dimWorld = GridView::dimensionworld
193 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
195 using MaterialLawParams =
typename MaterialLaw::Params;
197 using CoordScalar =
typename GridView::ctype;
198 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
199 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
200 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
207 : ParentType(simulator)
210 dofIsInLens_.resize(simulator.model().numGridDof());
218 ParentType::finishInit();
223 lensLowerLeft_[0] = 1.0;
224 lensLowerLeft_[1] = 2.0;
226 lensUpperRight_[0] = 4.0;
227 lensUpperRight_[1] = 3.0;
231 lensMaterialParams_.setVgAlpha(0.00045);
232 lensMaterialParams_.setVgN(7.3);
233 lensMaterialParams_.finalize();
235 outerMaterialParams_.setVgAlpha(0.0037);
236 outerMaterialParams_.setVgN(4.7);
237 outerMaterialParams_.finalize();
246 lensK_ = this->toDimMatrix_(1e-12);
247 outerK_ = this->toDimMatrix_(5e-12);
250 Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
251 auto elemIt = this->gridView().template begin<0>();
252 auto elemEndIt = this->gridView().template end<0>();
253 for (; elemIt != elemEndIt; ++elemIt) {
254 stencil.update(*elemIt);
255 for (
unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
256 unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
257 const auto& dofPos = stencil.subControlVolume(dofIdx).center();
258 dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
273 std::ostringstream oss;
274 oss <<
"lens_richards_"
275 << Model::discretizationName();
285 this->model().checkConservativeness();
289 this->model().globalStorage(storage);
292 if (this->gridView().comm().rank() == 0) {
293 std::cout <<
"Storage: " << storage << std::endl << std::flush;
301 template <
class Context>
302 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
303 {
return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
305 Scalar
temperature(
unsigned globalSpaceIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED)
const
306 {
return 273.15 + 10; }
311 template <
class Context>
314 unsigned timeIdx)
const
316 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
325 template <
class Context>
327 unsigned spaceIdx OPM_UNUSED,
328 unsigned timeIdx OPM_UNUSED)
const
334 template <
class Context>
337 unsigned timeIdx)
const
339 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
344 unsigned timeIdx OPM_UNUSED)
const
346 if (dofIsInLens_[globalSpaceIdx])
347 return lensMaterialParams_;
348 return outerMaterialParams_;
356 template <
class Context>
359 unsigned timeIdx)
const
360 {
return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
365 unsigned timeIdx OPM_UNUSED)
const
378 template <
class Context>
380 const Context& context,
382 unsigned timeIdx)
const
384 const auto& pos = context.pos(spaceIdx, timeIdx);
386 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
387 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
390 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
391 fs.setSaturation(wettingPhaseIdx, Sw);
392 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
395 MaterialLaw::capillaryPressures(pC, materialParams, fs);
396 fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
397 fs.setPressure(nonWettingPhaseIdx, pnRef_);
399 typename FluidSystem::template ParameterCache<Scalar> paramCache;
400 paramCache.updateAll(fs);
401 fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx));
404 fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx));
407 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
409 else if (onInlet_(pos)) {
410 RateVector massRate(0.0);
413 massRate[contiEqIdx] = -0.04;
415 values.setMassRate(massRate);
431 template <
class Context>
433 const Context& context,
435 unsigned timeIdx)
const
437 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
440 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
441 fs.setSaturation(wettingPhaseIdx, Sw);
442 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
445 MaterialLaw::capillaryPressures(pC, materialParams, fs);
446 values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
455 template <
class Context>
457 const Context& context OPM_UNUSED,
458 unsigned spaceIdx OPM_UNUSED,
459 unsigned timeIdx OPM_UNUSED)
const
460 { rate = Scalar(0.0); }
465 bool onLeftBoundary_(
const GlobalPosition& pos)
const
466 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
468 bool onRightBoundary_(
const GlobalPosition& pos)
const
469 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
471 bool onLowerBoundary_(
const GlobalPosition& pos)
const
472 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
474 bool onUpperBoundary_(
const GlobalPosition& pos)
const
475 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
477 bool onInlet_(
const GlobalPosition& pos)
const
479 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
480 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
481 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
484 bool isInLens_(
const GlobalPosition& pos)
const
486 for (
unsigned i = 0; i < dimWorld; ++i) {
487 if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
493 GlobalPosition lensLowerLeft_;
494 GlobalPosition lensUpperRight_;
498 MaterialLawParams lensMaterialParams_;
499 MaterialLawParams outerMaterialParams_;
501 std::vector<bool> dofIsInLens_;
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain.
Definition: richardslensproblem.hh:165
void finishInit()
Definition: richardslensproblem.hh:216
std::string name() const
Definition: richardslensproblem.hh:271
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:326
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:357
RichardsLensProblem(Simulator &simulator)
Definition: richardslensproblem.hh:206
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:379
void endTimeStep()
Definition: richardslensproblem.hh:282
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:456
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:312
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:302
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:432
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:335
Definition: richardslensproblem.hh:58