28#ifndef EWOMS_LENS_PROBLEM_HH
29#define EWOMS_LENS_PROBLEM_HH
31#include <dune/common/fmatrix.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/version.hh>
35#include <opm/material/components/Dnapl.hpp>
36#include <opm/material/components/SimpleH2O.hpp>
37#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
38#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
40#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
41#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
42#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
60template <
class TypeTag>
72template<
class TypeTag>
76template<
class TypeTag>
77struct Grid<TypeTag, TTag::LensBaseProblem> {
using type = Dune::YaspGrid<2>; };
80template<
class TypeTag>
87 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
91template<
class TypeTag>
98 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
102template<
class TypeTag>
107 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
108 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
111 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
112 FluidSystem::wettingPhaseIdx,
113 FluidSystem::nonWettingPhaseIdx>;
117 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
121 using type = Opm::EffToAbsLaw<EffectiveLaw>;
129template<
class Scalar>
130struct LensLowerLeftX {
static constexpr Scalar
value = 1.0; };
132template<
class Scalar>
133struct LensLowerLeftY {
static constexpr Scalar
value = 2.0; };
135template<
class Scalar>
136struct LensLowerLeftZ {
static constexpr Scalar
value = 0.0; };
138template<
class Scalar>
139struct LensUpperRightX {
static constexpr Scalar
value = 4.0; };
141template<
class Scalar>
142struct LensUpperRightY {
static constexpr Scalar
value = 3.0; };
144template<
class Scalar>
145struct LensUpperRightZ {
static constexpr Scalar
value = 1.0; };
174template <
class TypeTag>
191 numPhases = FluidSystem::numPhases,
194 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
195 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
198 contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
201 dim = GridView::dimension,
202 dimWorld = GridView::dimensionworld
211 using CoordScalar =
typename GridView::ctype;
212 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
214 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
221 : ParentType(simulator)
229 ParentType::finishInit();
234 temperature_ = 273.15 + 20;
235 lensLowerLeft_[0] = Parameters::Get<Parameters::LensLowerLeftX<Scalar>>();
236 lensLowerLeft_[1] = Parameters::Get<Parameters::LensLowerLeftY<Scalar>>();
237 lensUpperRight_[0] = Parameters::Get<Parameters::LensUpperRightX<Scalar>>();
238 lensUpperRight_[1] = Parameters::Get<Parameters::LensUpperRightY<Scalar>>();
240 if constexpr (dim == 3) {
241 lensLowerLeft_[2] = Parameters::Get<Parameters::LensLowerLeftZ<Scalar>>();
242 lensUpperRight_[2] = Parameters::Get<Parameters::LensUpperRightZ<Scalar>>();
246 lensMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.18);
247 lensMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
248 outerMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.05);
249 outerMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
252 lensMaterialParams_.setVgAlpha(0.00045);
253 lensMaterialParams_.setVgN(7.3);
254 outerMaterialParams_.setVgAlpha(0.0037);
255 outerMaterialParams_.setVgN(4.7);
257 lensMaterialParams_.finalize();
258 outerMaterialParams_.finalize();
260 lensK_ = this->toDimMatrix_(9.05e-12);
261 outerK_ = this->toDimMatrix_(4.6e-10);
265 this->gravity_[1] = -9.81;
274 ParentType::registerParameters();
276 Parameters::Register<Parameters::LensLowerLeftX<Scalar>>
277 (
"The x-coordinate of the lens' lower-left corner [m].");
278 Parameters::Register<Parameters::LensLowerLeftY<Scalar>>
279 (
"The y-coordinate of the lens' lower-left corner [m].");
280 Parameters::Register<Parameters::LensUpperRightX<Scalar>>
281 (
"The x-coordinate of the lens' upper-right corner [m].");
282 Parameters::Register<Parameters::LensUpperRightY<Scalar>>
283 (
"The y-coordinate of the lens' upper-right corner [m].");
285 if constexpr (dim == 3) {
286 Parameters::Register<Parameters::LensLowerLeftZ<Scalar>>
287 (
"The z-coordinate of the lens' lower-left corner [m].");
288 Parameters::Register<Parameters::LensUpperRightZ<Scalar>>
289 (
"The z-coordinate of the lens' upper-right corner [m].");
292 Parameters::SetDefault<Parameters::CellsX>(48);
293 Parameters::SetDefault<Parameters::CellsY>(32);
294 Parameters::SetDefault<Parameters::DomainSizeX<Scalar>>(6.0);
295 Parameters::SetDefault<Parameters::DomainSizeY<Scalar>>(4.0);
297 if constexpr (dim == 3) {
298 Parameters::SetDefault<Parameters::CellsZ>(16);
299 Parameters::SetDefault<Parameters::DomainSizeZ<Scalar>>(1.0);
304 constexpr bool useFD = std::is_same_v<LLS, Properties::TTag::FiniteDifferenceLocalLinearizer>;
305 if constexpr (useFD) {
306 Parameters::SetDefault<Parameters::NumericDifferenceMethod>(+1);
309 Parameters::SetDefault<Parameters::EndTime<Scalar>>(30e3);
310 Parameters::SetDefault<Parameters::EnableIntensiveQuantityCache>(
true);
311 Parameters::SetDefault<Parameters::EnableStorageCache>(
true);
312 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(250.0);
313 Parameters::SetDefault<Parameters::VtkWriteIntrinsicPermeabilities>(
true);
314 Parameters::SetDefault<Parameters::EnableGravity>(
true);
322 std::string thermal =
"isothermal";
323 constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
324 if constexpr (enableEnergy)
325 thermal =
"non-isothermal";
327 std::string deriv =
"finite difference";
329 constexpr bool useAutoDiff = std::is_same_v<LLS, Properties::TTag::AutoDiffLocalLinearizer>;
330 if constexpr (useAutoDiff) {
331 deriv =
"automatic differentiation";
334 std::string disc =
"vertex centered finite volume";
336 constexpr bool useEcfv = std::is_same<D, Opm::EcfvDiscretization<TypeTag>>::value;
337 if constexpr (useEcfv)
338 disc =
"element centered finite volume";
340 return std::string(
"")+
341 "Ground remediation problem where a dense oil infiltrates "+
342 "an aquifer with an embedded low-permability lens. " +
343 "This is the binary for the "+thermal+
" variant using "+deriv+
344 "and the "+disc+
" discretization";
355 template <
class Context>
357 unsigned timeIdx)
const
359 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
361 if (isInLens_(globalPos))
369 template <
class Context>
378 template <
class Context>
380 unsigned spaceIdx,
unsigned timeIdx)
const
382 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
384 if (isInLens_(globalPos))
385 return lensMaterialParams_;
386 return outerMaterialParams_;
392 template <
class Context>
396 {
return temperature_; }
412 constexpr bool useAutoDiff = std::is_same_v<LLS, Properties::TTag::AutoDiffLocalLinearizer>;
415 constexpr bool useTrans = std::is_same_v<FM, Opm::TransFluxModule<TypeTag>>;
417 std::ostringstream oss;
418 oss <<
"lens_" << Model::name()
419 <<
"_" << Model::discretizationName()
420 <<
"_" << (useAutoDiff?
"ad":
"fd");
449 this->model().globalStorage(storage);
452 if (this->gridView().comm().rank() == 0) {
453 std::cout <<
"Storage: " << storage << std::endl << std::flush;
468 template <
class Context>
470 const Context& context,
472 unsigned timeIdx)
const
474 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
476 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
478 Scalar densityW = WettingPhase::density(temperature_, Scalar(1e5));
479 Scalar densityN = NonwettingPhase::density(temperature_, Scalar(1e5));
481 Scalar T =
temperature(context, spaceIdx, timeIdx);
485 if (onLeftBoundary_(pos)) {
486 Scalar height = this->boundingBoxMax()[1] - this->boundingBoxMin()[1];
487 Scalar depth = this->boundingBoxMax()[1] - pos[1];
488 Scalar alpha = (1 + 1.5 / height);
491 pw = 1e5 - alpha * densityW * this->gravity()[1] * depth;
495 Scalar depth = this->boundingBoxMax()[1] - pos[1];
498 pw = 1e5 - densityW * this->gravity()[1] * depth;
503 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
505 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
507 fs.setSaturation(wettingPhaseIdx, Sw);
508 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
509 fs.setTemperature(T);
511 Scalar pC[numPhases];
512 MaterialLaw::capillaryPressures(pC, matParams, fs);
513 fs.setPressure(wettingPhaseIdx, pw);
514 fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]);
516 fs.setDensity(wettingPhaseIdx, densityW);
517 fs.setDensity(nonWettingPhaseIdx, densityN);
519 fs.setViscosity(wettingPhaseIdx, WettingPhase::viscosity(temperature_, fs.pressure(wettingPhaseIdx)));
520 fs.setViscosity(nonWettingPhaseIdx, NonwettingPhase::viscosity(temperature_, fs.pressure(nonWettingPhaseIdx)));
523 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
525 else if (onInlet_(pos)) {
526 RateVector massRate(0.0);
528 massRate[contiNEqIdx] = -0.04;
531 values.setMassRate(massRate);
549 template <
class Context>
550 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
552 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
553 Scalar depth = this->boundingBoxMax()[1] - pos[1];
555 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
556 fs.setPressure(wettingPhaseIdx, 1e5);
559 fs.setSaturation(wettingPhaseIdx, Sw);
560 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
562 fs.setTemperature(temperature_);
564 typename FluidSystem::template ParameterCache<Scalar> paramCache;
565 paramCache.updatePhase(fs, wettingPhaseIdx);
566 Scalar densityW = FluidSystem::density(fs, paramCache, wettingPhaseIdx);
569 Scalar pw = 1e5 - densityW * this->gravity()[1] * depth;
572 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
573 Scalar pC[numPhases];
574 MaterialLaw::capillaryPressures(pC, matParams, fs);
577 fs.setPressure(wettingPhaseIdx, pw);
578 fs.setPressure(nonWettingPhaseIdx, pw + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]));
581 values.assignNaive(fs);
590 template <
class Context>
595 { rate = Scalar(0.0); }
600 bool isInLens_(
const GlobalPosition& pos)
const
602 for (
unsigned i = 0; i < dim; ++i) {
603 if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i]
610 bool onLeftBoundary_(
const GlobalPosition& pos)
const
611 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
613 bool onRightBoundary_(
const GlobalPosition& pos)
const
614 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
616 bool onLowerBoundary_(
const GlobalPosition& pos)
const
617 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
619 bool onUpperBoundary_(
const GlobalPosition& pos)
const
620 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
622 bool onInlet_(
const GlobalPosition& pos)
const
624 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
625 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
626 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
629 GlobalPosition lensLowerLeft_;
630 GlobalPosition lensUpperRight_;
634 MaterialLawParams lensMaterialParams_;
635 MaterialLawParams outerMaterialParams_;
Soil contamination problem where DNAPL infiltrates a fully water saturated medium.
Definition: lensproblem.hh:176
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: lensproblem.hh:393
static void registerParameters()
Definition: lensproblem.hh:272
void beginTimeStep()
Called by the simulator before each time integration.
Definition: lensproblem.hh:430
static std::string briefDescription()
Returns a human readable description of the problem for the help message.
Definition: lensproblem.hh:320
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:379
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: lensproblem.hh:550
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: lensproblem.hh:370
LensProblem(Simulator &simulator)
Definition: lensproblem.hh:220
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: lensproblem.hh:227
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: lensproblem.hh:469
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:356
std::string name() const
The problem name.
Definition: lensproblem.hh:408
void endTimeStep()
Called by the simulator after each time integration.
Definition: lensproblem.hh:442
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition: lensproblem.hh:436
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: lensproblem.hh:591
Defines the properties required for the immiscible multi-phase model.
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilnewtonmethodparams.hpp:31
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
static constexpr Scalar value
Definition: groundwaterproblem.hh:95
static constexpr Scalar value
Definition: groundwaterproblem.hh:98
static constexpr Scalar value
Definition: groundwaterproblem.hh:101
static constexpr Scalar value
Definition: groundwaterproblem.hh:104
static constexpr Scalar value
Definition: groundwaterproblem.hh:107
static constexpr Scalar value
Definition: groundwaterproblem.hh:110
Dune::YaspGrid< 2 > type
Definition: lensproblem.hh:77
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::EffToAbsLaw< EffectiveLaw > type
Definition: lensproblem.hh:121
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
Opm::LiquidPhase< Scalar, Opm::DNAPL< Scalar > > type
Definition: lensproblem.hh:98
The non-wetting phase for two-phase models.
Definition: immiscibleproperties.hh:44
The type of the problem.
Definition: fvbaseproperties.hh:81
Definition: lensproblem.hh:68
std::tuple< StructuredGridVanguard > InheritsFrom
Definition: lensproblem.hh:68
Opm::LiquidPhase< Scalar, Opm::SimpleH2O< Scalar > > type
Definition: lensproblem.hh:87
The wetting phase for two-phase models.
Definition: immiscibleproperties.hh:41
This file contains the flux module that uses transmissibilities.