28#ifndef EWOMS_RICHARDS_LENS_PROBLEM_HH
29#define EWOMS_RICHARDS_LENS_PROBLEM_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>
48template <
class TypeTag>
49class RichardsLensProblem;
61template<
class TypeTag>
65template<
class TypeTag>
69template<
class TypeTag>
76 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
80template<
class TypeTag>
85 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
86 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
89 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
90 FluidSystem::wettingPhaseIdx,
91 FluidSystem::nonWettingPhaseIdx>;
95 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
99 using type = Opm::EffToAbsLaw<EffectiveLaw>;
122template <
class TypeTag>
141 pressureWIdx = Indices::pressureWIdx,
142 contiEqIdx = Indices::contiEqIdx,
143 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
144 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
145 numPhases = FluidSystem::numPhases,
148 dimWorld = GridView::dimensionworld
154 using MaterialLawParams =
typename MaterialLaw::Params;
156 using CoordScalar =
typename GridView::ctype;
157 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
158 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
159 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
166 : ParentType(simulator)
169 dofIsInLens_.resize(simulator.model().numGridDof());
177 ParentType::finishInit();
182 lensLowerLeft_[0] = 1.0;
183 lensLowerLeft_[1] = 2.0;
185 lensUpperRight_[0] = 4.0;
186 lensUpperRight_[1] = 3.0;
190 lensMaterialParams_.setVgAlpha(0.00045);
191 lensMaterialParams_.setVgN(7.3);
192 lensMaterialParams_.finalize();
194 outerMaterialParams_.setVgAlpha(0.0037);
195 outerMaterialParams_.setVgN(4.7);
196 outerMaterialParams_.finalize();
205 lensK_ = this->toDimMatrix_(1e-12);
206 outerK_ = this->toDimMatrix_(5e-12);
209 Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
210 for (
const auto& elem : elements(this->gridView())) {
211 stencil.update(elem);
212 for (
unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
213 unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
214 const auto& dofPos = stencil.subControlVolume(dofIdx).center();
215 dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
225 ParentType::registerParameters();
227 Parameters::SetDefault<Parameters::GridFile>(
"./data/richardslens_24x16.dgf");
231 constexpr bool useFD = std::is_same_v<LLS, Properties::TTag::FiniteDifferenceLocalLinearizer>;
232 if constexpr (useFD) {
233 Parameters::SetDefault<Parameters::NumericDifferenceMethod>(0);
236 Parameters::SetDefault<Parameters::EndTime<Scalar>>(3000.0);
237 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(100.0);
238 Parameters::SetDefault<Parameters::NewtonMaxIterations>(28);
239 Parameters::SetDefault<Parameters::NewtonTargetIterations>(18);
240 Parameters::SetDefault<Parameters::EnableGravity>(
true);
253 std::ostringstream oss;
254 oss <<
"lens_richards_"
255 << Model::discretizationName();
265 this->model().checkConservativeness();
269 this->model().globalStorage(storage);
272 if (this->gridView().comm().rank() == 0) {
273 std::cout <<
"Storage: " << storage << std::endl << std::flush;
281 template <
class Context>
282 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
283 {
return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
286 {
return 273.15 + 10; }
291 template <
class Context>
294 unsigned timeIdx)
const
296 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
305 template <
class Context>
314 template <
class Context>
317 unsigned timeIdx)
const
319 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
326 if (dofIsInLens_[globalSpaceIdx])
327 return lensMaterialParams_;
328 return outerMaterialParams_;
336 template <
class Context>
339 unsigned timeIdx)
const
340 {
return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
358 template <
class Context>
360 const Context& context,
362 unsigned timeIdx)
const
364 const auto& pos = context.pos(spaceIdx, timeIdx);
366 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
367 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
370 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
371 fs.setSaturation(wettingPhaseIdx, Sw);
372 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
375 MaterialLaw::capillaryPressures(pC, materialParams, fs);
376 fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
377 fs.setPressure(nonWettingPhaseIdx, pnRef_);
379 typename FluidSystem::template ParameterCache<Scalar> paramCache;
380 paramCache.updateAll(fs);
381 fs.setDensity(wettingPhaseIdx, FluidSystem::density(fs, paramCache, wettingPhaseIdx));
384 fs.setViscosity(wettingPhaseIdx, FluidSystem::viscosity(fs, paramCache, wettingPhaseIdx));
387 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
389 else if (onInlet_(pos)) {
390 RateVector massRate(0.0);
393 massRate[contiEqIdx] = -0.04;
395 values.setMassRate(massRate);
411 template <
class Context>
413 const Context& context,
415 unsigned timeIdx)
const
417 const auto& materialParams = this->
materialLawParams(context, spaceIdx, timeIdx);
420 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
421 fs.setSaturation(wettingPhaseIdx, Sw);
422 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
425 MaterialLaw::capillaryPressures(pC, materialParams, fs);
426 values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
435 template <
class Context>
440 { rate = Scalar(0.0); }
445 bool onLeftBoundary_(
const GlobalPosition& pos)
const
446 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
448 bool onRightBoundary_(
const GlobalPosition& pos)
const
449 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
451 bool onLowerBoundary_(
const GlobalPosition& pos)
const
452 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
454 bool onUpperBoundary_(
const GlobalPosition& pos)
const
455 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
457 bool onInlet_(
const GlobalPosition& pos)
const
459 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
460 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
461 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
464 bool isInLens_(
const GlobalPosition& pos)
const
466 for (
unsigned i = 0; i < dimWorld; ++i) {
467 if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
473 GlobalPosition lensLowerLeft_;
474 GlobalPosition lensUpperRight_;
478 MaterialLawParams lensMaterialParams_;
479 MaterialLawParams outerMaterialParams_;
481 std::vector<bool> dofIsInLens_;
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain.
Definition: richardslensproblem.hh:124
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: richardslensproblem.hh:175
static void registerParameters()
Definition: richardslensproblem.hh:223
std::string name() const
The problem name.
Definition: richardslensproblem.hh:251
const MaterialLawParams & materialLawParams(unsigned globalSpaceIdx, unsigned) const
Definition: richardslensproblem.hh:323
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:315
Scalar temperature(unsigned, unsigned) const
Definition: richardslensproblem.hh:285
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:337
RichardsLensProblem(Simulator &simulator)
Definition: richardslensproblem.hh:165
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:292
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: richardslensproblem.hh:436
Scalar referencePressure(unsigned, unsigned) const
Definition: richardslensproblem.hh:344
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: richardslensproblem.hh:359
void endTimeStep()
Called by the simulator after each time integration.
Definition: richardslensproblem.hh:262
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:282
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: richardslensproblem.hh:412
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: richardslensproblem.hh:306
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
Dune::YaspGrid< 2 > type
Definition: richardslensproblem.hh:62
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::EffToAbsLaw< EffectiveLaw > type
Definition: richardslensproblem.hh:99
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Definition: richardslensproblem.hh:57
std::tuple< Richards > InheritsFrom
Definition: richardslensproblem.hh:57
Opm::LiquidPhase< Scalar, Opm::SimpleH2O< Scalar > > type
Definition: richardslensproblem.hh:76
Definition: richardsproperties.hh:40