28 #ifndef EWOMS_LENS_PROBLEM_HH
29 #define EWOMS_LENS_PROBLEM_HH
31 #include <opm/models/io/structuredgridvanguard.hh>
32 #include <opm/models/immiscible/immiscibleproperties.hh>
33 #include <opm/models/discretization/common/fvbaseadlocallinearizer.hh>
34 #include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
35 #include <opm/models/common/transfluxmodule.hh>
36 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
37 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
38 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
39 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40 #include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
41 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
42 #include <opm/material/components/SimpleH2O.hpp>
43 #include <opm/material/components/Dnapl.hpp>
45 #include <dune/common/version.hh>
46 #include <dune/common/fvector.hh>
47 #include <dune/common/fmatrix.hh>
54 template <
class TypeTag>
58 namespace Opm::Properties {
62 struct LensBaseProblem {
using InheritsFrom = std::tuple<StructuredGridVanguard>; };
66 template<
class TypeTag,
class MyTypeTag>
68 template<
class TypeTag,
class MyTypeTag>
69 struct LensLowerLeftY {
using type = UndefinedProperty; };
70 template<
class TypeTag,
class MyTypeTag>
71 struct LensLowerLeftZ {
using type = UndefinedProperty; };
72 template<
class TypeTag,
class MyTypeTag>
73 struct LensUpperRightX {
using type = UndefinedProperty; };
74 template<
class TypeTag,
class MyTypeTag>
75 struct LensUpperRightY {
using type = UndefinedProperty; };
76 template<
class TypeTag,
class MyTypeTag>
77 struct LensUpperRightZ {
using type = UndefinedProperty; };
80 template<
class TypeTag>
84 template<
class TypeTag>
85 struct Grid<TypeTag, TTag::LensBaseProblem> {
using type = Dune::YaspGrid<2>; };
88 template<
class TypeTag>
89 struct WettingPhase<TypeTag, TTag::LensBaseProblem>
92 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
95 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
99 template<
class TypeTag>
100 struct NonwettingPhase<TypeTag, TTag::LensBaseProblem>
103 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
106 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
110 template<
class TypeTag>
111 struct MaterialLaw<TypeTag, TTag::LensBaseProblem>
114 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
115 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
116 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
118 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
119 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
120 FluidSystem::wettingPhaseIdx,
121 FluidSystem::nonWettingPhaseIdx>;
125 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
129 using type = Opm::EffToAbsLaw<EffectiveLaw>;
133 template<
class TypeTag>
134 struct NewtonWriteConvergence<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
false; };
137 template<
class TypeTag>
138 struct NumericDifferenceMethod<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = +1; };
141 template<
class TypeTag>
142 struct EnableGravity<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
145 template<
class TypeTag>
148 using type = GetPropType<TypeTag, Scalar>;
149 static constexpr type value = 1.0;
151 template<
class TypeTag>
154 using type = GetPropType<TypeTag, Scalar>;
155 static constexpr type value = 2.0;
157 template<
class TypeTag>
160 using type = GetPropType<TypeTag, Scalar>;
161 static constexpr type value = 0.0;
163 template<
class TypeTag>
166 using type = GetPropType<TypeTag, Scalar>;
167 static constexpr type value = 4.0;
169 template<
class TypeTag>
172 using type = GetPropType<TypeTag, Scalar>;
173 static constexpr type value = 3.0;
175 template<
class TypeTag>
178 using type = GetPropType<TypeTag, Scalar>;
179 static constexpr type value = 1.0;
182 template<
class TypeTag>
183 struct DomainSizeX<TypeTag, TTag::LensBaseProblem>
185 using type = GetPropType<TypeTag, Scalar>;
186 static constexpr type value = 6.0;
188 template<
class TypeTag>
189 struct DomainSizeY<TypeTag, TTag::LensBaseProblem>
191 using type = GetPropType<TypeTag, Scalar>;
192 static constexpr type value = 4.0;
194 template<
class TypeTag>
195 struct DomainSizeZ<TypeTag, TTag::LensBaseProblem>
197 using type = GetPropType<TypeTag, Scalar>;
198 static constexpr type value = 1.0;
201 template<
class TypeTag>
202 struct CellsX<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 48; };
203 template<
class TypeTag>
204 struct CellsY<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 32; };
205 template<
class TypeTag>
206 struct CellsZ<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 16; };
209 template<
class TypeTag>
210 struct EndTime<TypeTag, TTag::LensBaseProblem>
212 using type = GetPropType<TypeTag, Scalar>;
213 static constexpr type value = 30e3;
217 template<
class TypeTag>
218 struct InitialTimeStepSize<TypeTag, TTag::LensBaseProblem>
220 using type = GetPropType<TypeTag, Scalar>;
221 static constexpr type value = 250;
225 template<
class TypeTag>
226 struct VtkWriteIntrinsicPermeabilities<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
229 template<
class TypeTag>
230 struct EnableStorageCache<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
233 template<
class TypeTag>
234 struct EnableIntensiveQuantityCache<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
263 template <
class TypeTag>
264 class LensProblem :
public GetPropType<TypeTag, Properties::BaseProblem>
266 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
268 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
269 using GridView = GetPropType<TypeTag, Properties::GridView>;
270 using Indices = GetPropType<TypeTag, Properties::Indices>;
271 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
272 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
273 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
274 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
275 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
276 using Model = GetPropType<TypeTag, Properties::Model>;
280 numPhases = FluidSystem::numPhases,
283 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
284 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
287 contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
290 dim = GridView::dimension,
291 dimWorld = GridView::dimensionworld
294 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
295 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
296 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
297 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
298 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
300 using CoordScalar =
typename GridView::ctype;
301 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
303 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
310 : ParentType(simulator)
318 ParentType::finishInit();
323 temperature_ = 273.15 + 20;
324 lensLowerLeft_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftX);
325 lensLowerLeft_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
326 lensUpperRight_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX);
327 lensUpperRight_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
330 lensLowerLeft_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftZ);
331 lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightZ);
335 lensMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.18);
336 lensMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
337 outerMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.05);
338 outerMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
341 lensMaterialParams_.setVgAlpha(0.00045);
342 lensMaterialParams_.setVgN(7.3);
343 outerMaterialParams_.setVgAlpha(0.0037);
344 outerMaterialParams_.setVgN(4.7);
346 lensMaterialParams_.finalize();
347 outerMaterialParams_.finalize();
349 lensK_ = this->toDimMatrix_(9.05e-12);
350 outerK_ = this->toDimMatrix_(4.6e-10);
354 this->gravity_[1] = -9.81;
363 ParentType::registerParameters();
365 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX,
366 "The x-coordinate of the lens' lower-left corner "
368 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY,
369 "The y-coordinate of the lens' lower-left corner "
371 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX,
372 "The x-coordinate of the lens' upper-right corner "
374 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY,
375 "The y-coordinate of the lens' upper-right corner "
379 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ,
380 "The z-coordinate of the lens' lower-left "
382 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ,
383 "The z-coordinate of the lens' upper-right "
393 std::string thermal =
"isothermal";
394 bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
396 thermal =
"non-isothermal";
398 std::string deriv =
"finite difference";
399 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
400 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
402 deriv =
"automatic differentiation";
404 std::string disc =
"vertex centered finite volume";
405 using D = GetPropType<TypeTag, Properties::Discretization>;
406 bool useEcfv = std::is_same<D, Opm::EcfvDiscretization<TypeTag>>::value;
408 disc =
"element centered finite volume";
410 return std::string(
"")+
411 "Ground remediation problem where a dense oil infiltrates "+
412 "an aquifer with an embedded low-permability lens. " +
413 "This is the binary for the "+thermal+
" variant using "+deriv+
414 "and the "+disc+
" discretization";
425 template <
class Context>
427 unsigned timeIdx)
const
429 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
431 if (isInLens_(globalPos))
439 template <
class Context>
448 template <
class Context>
450 unsigned spaceIdx,
unsigned timeIdx)
const
452 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
454 if (isInLens_(globalPos))
455 return lensMaterialParams_;
456 return outerMaterialParams_;
462 template <
class Context>
466 {
return temperature_; }
480 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
482 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
484 using FM = GetPropType<TypeTag, Properties::FluxModule>;
485 bool useTrans = std::is_same<FM, Opm::TransFluxModule<TypeTag>>::value;
487 std::ostringstream oss;
488 oss <<
"lens_" << Model::name()
489 <<
"_" << Model::discretizationName()
490 <<
"_" << (useAutoDiff?
"ad":
"fd");
519 this->model().globalStorage(storage);
522 if (this->gridView().comm().rank() == 0) {
523 std::cout <<
"Storage: " << storage << std::endl << std::flush;
538 template <
class Context>
540 const Context& context,
542 unsigned timeIdx)
const
544 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
546 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
548 Scalar densityW = WettingPhase::density(temperature_, Scalar(1e5));
549 Scalar densityN = NonwettingPhase::density(temperature_, Scalar(1e5));
551 Scalar T =
temperature(context, spaceIdx, timeIdx);
555 if (onLeftBoundary_(pos)) {
556 Scalar height = this->boundingBoxMax()[1] - this->boundingBoxMin()[1];
557 Scalar depth = this->boundingBoxMax()[1] - pos[1];
558 Scalar alpha = (1 + 1.5 / height);
561 pw = 1e5 - alpha * densityW * this->gravity()[1] * depth;
565 Scalar depth = this->boundingBoxMax()[1] - pos[1];
568 pw = 1e5 - densityW * this->gravity()[1] * depth;
573 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
575 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
577 fs.setSaturation(wettingPhaseIdx, Sw);
578 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
579 fs.setTemperature(T);
581 Scalar pC[numPhases];
582 MaterialLaw::capillaryPressures(pC, matParams, fs);
583 fs.setPressure(wettingPhaseIdx, pw);
584 fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]);
586 fs.setDensity(wettingPhaseIdx, densityW);
587 fs.setDensity(nonWettingPhaseIdx, densityN);
589 fs.setViscosity(wettingPhaseIdx, WettingPhase::viscosity(temperature_, fs.pressure(wettingPhaseIdx)));
590 fs.setViscosity(nonWettingPhaseIdx, NonwettingPhase::viscosity(temperature_, fs.pressure(nonWettingPhaseIdx)));
593 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
595 else if (onInlet_(pos)) {
596 RateVector massRate(0.0);
598 massRate[contiNEqIdx] = -0.04;
601 values.setMassRate(massRate);
619 template <
class Context>
620 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
622 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
623 Scalar depth = this->boundingBoxMax()[1] - pos[1];
625 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
626 fs.setPressure(wettingPhaseIdx, 1e5);
629 fs.setSaturation(wettingPhaseIdx, Sw);
630 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
632 fs.setTemperature(temperature_);
634 typename FluidSystem::template ParameterCache<Scalar> paramCache;
635 paramCache.updatePhase(fs, wettingPhaseIdx);
636 Scalar densityW = FluidSystem::density(fs, paramCache, wettingPhaseIdx);
639 Scalar pw = 1e5 - densityW * this->gravity()[1] * depth;
642 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
643 Scalar pC[numPhases];
644 MaterialLaw::capillaryPressures(pC, matParams, fs);
647 fs.setPressure(wettingPhaseIdx, pw);
648 fs.setPressure(nonWettingPhaseIdx, pw + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]));
651 values.assignNaive(fs);
660 template <
class Context>
665 { rate = Scalar(0.0); }
670 bool isInLens_(
const GlobalPosition& pos)
const
672 for (
unsigned i = 0; i < dim; ++i) {
673 if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i]
680 bool onLeftBoundary_(
const GlobalPosition& pos)
const
681 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
683 bool onRightBoundary_(
const GlobalPosition& pos)
const
684 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
686 bool onLowerBoundary_(
const GlobalPosition& pos)
const
687 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
689 bool onUpperBoundary_(
const GlobalPosition& pos)
const
690 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
692 bool onInlet_(
const GlobalPosition& pos)
const
694 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
695 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
696 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
699 GlobalPosition lensLowerLeft_;
700 GlobalPosition lensUpperRight_;
704 MaterialLawParams lensMaterialParams_;
705 MaterialLawParams outerMaterialParams_;
Soil contamination problem where DNAPL infiltrates a fully water saturated medium.
Definition: lensproblem.hh:265
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: lensproblem.hh:463
static void registerParameters()
Definition: lensproblem.hh:361
void beginTimeStep()
Definition: lensproblem.hh:500
static std::string briefDescription()
Definition: lensproblem.hh:391
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:620
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:449
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: lensproblem.hh:440
LensProblem(Simulator &simulator)
Definition: lensproblem.hh:309
void finishInit()
Definition: lensproblem.hh:316
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:539
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:426
std::string name() const
Definition: lensproblem.hh:478
void endTimeStep()
Definition: lensproblem.hh:512
void beginIteration()
Definition: lensproblem.hh:506
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: lensproblem.hh:661
Definition: groundwaterproblem.hh:60
Definition: groundwaterproblem.hh:62
Definition: groundwaterproblem.hh:64
Definition: groundwaterproblem.hh:66
Definition: groundwaterproblem.hh:68
Definition: groundwaterproblem.hh:70
Definition: lensproblem.hh:62