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>
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>
44 #include <opm/material/common/Unused.hpp>
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
55 template <
class TypeTag>
59 namespace Opm::Properties {
63 struct LensBaseProblem {
using InheritsFrom = std::tuple<StructuredGridVanguard>; };
67 template<
class TypeTag,
class MyTypeTag>
69 template<
class TypeTag,
class MyTypeTag>
70 struct LensLowerLeftY {
using type = UndefinedProperty; };
71 template<
class TypeTag,
class MyTypeTag>
72 struct LensLowerLeftZ {
using type = UndefinedProperty; };
73 template<
class TypeTag,
class MyTypeTag>
74 struct LensUpperRightX {
using type = UndefinedProperty; };
75 template<
class TypeTag,
class MyTypeTag>
76 struct LensUpperRightY {
using type = UndefinedProperty; };
77 template<
class TypeTag,
class MyTypeTag>
78 struct LensUpperRightZ {
using type = UndefinedProperty; };
81 template<
class TypeTag>
85 template<
class TypeTag>
86 struct Grid<TypeTag, TTag::LensBaseProblem> {
using type = Dune::YaspGrid<2>; };
89 template<
class TypeTag>
90 struct WettingPhase<TypeTag, TTag::LensBaseProblem>
93 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
96 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
100 template<
class TypeTag>
101 struct NonwettingPhase<TypeTag, TTag::LensBaseProblem>
104 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
107 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
111 template<
class TypeTag>
112 struct MaterialLaw<TypeTag, TTag::LensBaseProblem>
115 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
116 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
117 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
119 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
120 using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
121 FluidSystem::wettingPhaseIdx,
122 FluidSystem::nonWettingPhaseIdx>;
126 using EffectiveLaw = Opm::RegularizedVanGenuchten<Traits>;
130 using type = Opm::EffToAbsLaw<EffectiveLaw>;
134 template<
class TypeTag>
135 struct NewtonWriteConvergence<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
false; };
138 template<
class TypeTag>
139 struct NumericDifferenceMethod<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = +1; };
142 template<
class TypeTag>
143 struct EnableGravity<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
146 template<
class TypeTag>
149 using type = GetPropType<TypeTag, Scalar>;
150 static constexpr type value = 1.0;
152 template<
class TypeTag>
155 using type = GetPropType<TypeTag, Scalar>;
156 static constexpr type value = 2.0;
158 template<
class TypeTag>
161 using type = GetPropType<TypeTag, Scalar>;
162 static constexpr type value = 0.0;
164 template<
class TypeTag>
167 using type = GetPropType<TypeTag, Scalar>;
168 static constexpr type value = 4.0;
170 template<
class TypeTag>
173 using type = GetPropType<TypeTag, Scalar>;
174 static constexpr type value = 3.0;
176 template<
class TypeTag>
179 using type = GetPropType<TypeTag, Scalar>;
180 static constexpr type value = 1.0;
183 template<
class TypeTag>
184 struct DomainSizeX<TypeTag, TTag::LensBaseProblem>
186 using type = GetPropType<TypeTag, Scalar>;
187 static constexpr type value = 6.0;
189 template<
class TypeTag>
190 struct DomainSizeY<TypeTag, TTag::LensBaseProblem>
192 using type = GetPropType<TypeTag, Scalar>;
193 static constexpr type value = 4.0;
195 template<
class TypeTag>
196 struct DomainSizeZ<TypeTag, TTag::LensBaseProblem>
198 using type = GetPropType<TypeTag, Scalar>;
199 static constexpr type value = 1.0;
202 template<
class TypeTag>
203 struct CellsX<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 48; };
204 template<
class TypeTag>
205 struct CellsY<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 32; };
206 template<
class TypeTag>
207 struct CellsZ<TypeTag, TTag::LensBaseProblem> {
static constexpr
int value = 16; };
210 template<
class TypeTag>
211 struct EndTime<TypeTag, TTag::LensBaseProblem>
213 using type = GetPropType<TypeTag, Scalar>;
214 static constexpr type value = 30e3;
218 template<
class TypeTag>
219 struct InitialTimeStepSize<TypeTag, TTag::LensBaseProblem>
221 using type = GetPropType<TypeTag, Scalar>;
222 static constexpr type value = 250;
226 template<
class TypeTag>
227 struct VtkWriteIntrinsicPermeabilities<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
230 template<
class TypeTag>
231 struct EnableStorageCache<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
234 template<
class TypeTag>
235 struct EnableIntensiveQuantityCache<TypeTag, TTag::LensBaseProblem> {
static constexpr
bool value =
true; };
264 template <
class TypeTag>
265 class LensProblem :
public GetPropType<TypeTag, Properties::BaseProblem>
267 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
269 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
270 using GridView = GetPropType<TypeTag, Properties::GridView>;
271 using Indices = GetPropType<TypeTag, Properties::Indices>;
272 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
273 using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
274 using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
275 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
276 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
277 using Model = GetPropType<TypeTag, Properties::Model>;
281 numPhases = FluidSystem::numPhases,
284 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
285 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
288 contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
291 dim = GridView::dimension,
292 dimWorld = GridView::dimensionworld
295 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
296 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
297 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
298 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
299 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
301 using CoordScalar =
typename GridView::ctype;
302 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
304 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
311 : ParentType(simulator)
319 ParentType::finishInit();
324 temperature_ = 273.15 + 20;
325 lensLowerLeft_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftX);
326 lensLowerLeft_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftY);
327 lensUpperRight_[0] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX);
328 lensUpperRight_[1] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
331 lensLowerLeft_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensLowerLeftZ);
332 lensUpperRight_[2] = EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightZ);
336 lensMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.18);
337 lensMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
338 outerMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.05);
339 outerMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
342 lensMaterialParams_.setVgAlpha(0.00045);
343 lensMaterialParams_.setVgN(7.3);
344 outerMaterialParams_.setVgAlpha(0.0037);
345 outerMaterialParams_.setVgN(4.7);
347 lensMaterialParams_.finalize();
348 outerMaterialParams_.finalize();
350 lensK_ = this->toDimMatrix_(9.05e-12);
351 outerK_ = this->toDimMatrix_(4.6e-10);
355 this->gravity_[1] = -9.81;
364 ParentType::registerParameters();
366 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftX,
367 "The x-coordinate of the lens' lower-left corner "
369 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftY,
370 "The y-coordinate of the lens' lower-left corner "
372 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightX,
373 "The x-coordinate of the lens' upper-right corner "
375 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightY,
376 "The y-coordinate of the lens' upper-right corner "
380 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensLowerLeftZ,
381 "The z-coordinate of the lens' lower-left "
383 EWOMS_REGISTER_PARAM(TypeTag, Scalar, LensUpperRightZ,
384 "The z-coordinate of the lens' upper-right "
394 std::string thermal =
"isothermal";
395 bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
397 thermal =
"non-isothermal";
399 std::string deriv =
"finite difference";
400 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
401 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
403 deriv =
"automatic differentiation";
405 std::string disc =
"vertex centered finite volume";
406 using D = GetPropType<TypeTag, Properties::Discretization>;
407 bool useEcfv = std::is_same<D, Opm::EcfvDiscretization<TypeTag>>::value;
409 disc =
"element centered finite volume";
411 return std::string(
"")+
412 "Ground remediation problem where a dense oil infiltrates "+
413 "an aquifer with an embedded low-permability lens. " +
414 "This is the binary for the "+thermal+
" variant using "+deriv+
415 "and the "+disc+
" discretization";
426 template <
class Context>
428 unsigned timeIdx)
const
430 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
432 if (isInLens_(globalPos))
440 template <
class Context>
442 unsigned spaceIdx OPM_UNUSED,
443 unsigned timeIdx OPM_UNUSED)
const
449 template <
class Context>
451 unsigned spaceIdx,
unsigned timeIdx)
const
453 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
455 if (isInLens_(globalPos))
456 return lensMaterialParams_;
457 return outerMaterialParams_;
463 template <
class Context>
465 unsigned spaceIdx OPM_UNUSED,
466 unsigned timeIdx OPM_UNUSED)
const
467 {
return temperature_; }
481 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
483 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
485 std::ostringstream oss;
486 oss <<
"lens_" << Model::name()
487 <<
"_" << Model::discretizationName()
488 <<
"_" << (useAutoDiff?
"ad":
"fd");
510 this->model().checkConservativeness();
514 this->model().globalStorage(storage);
517 if (this->gridView().comm().rank() == 0) {
518 std::cout <<
"Storage: " << storage << std::endl << std::flush;
533 template <
class Context>
535 const Context& context,
537 unsigned timeIdx)
const
539 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
541 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
543 Scalar densityW = WettingPhase::density(temperature_, Scalar(1e5));
544 Scalar densityN = NonwettingPhase::density(temperature_, Scalar(1e5));
546 Scalar T =
temperature(context, spaceIdx, timeIdx);
550 if (onLeftBoundary_(pos)) {
551 Scalar height = this->boundingBoxMax()[1] - this->boundingBoxMin()[1];
552 Scalar depth = this->boundingBoxMax()[1] - pos[1];
553 Scalar alpha = (1 + 1.5 / height);
556 pw = 1e5 - alpha * densityW * this->gravity()[1] * depth;
560 Scalar depth = this->boundingBoxMax()[1] - pos[1];
563 pw = 1e5 - densityW * this->gravity()[1] * depth;
568 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
570 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
572 fs.setSaturation(wettingPhaseIdx, Sw);
573 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
574 fs.setTemperature(T);
576 Scalar pC[numPhases];
577 MaterialLaw::capillaryPressures(pC, matParams, fs);
578 fs.setPressure(wettingPhaseIdx, pw);
579 fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]);
581 fs.setDensity(wettingPhaseIdx, densityW);
582 fs.setDensity(nonWettingPhaseIdx, densityN);
584 fs.setViscosity(wettingPhaseIdx, WettingPhase::viscosity(temperature_, fs.pressure(wettingPhaseIdx)));
585 fs.setViscosity(nonWettingPhaseIdx, NonwettingPhase::viscosity(temperature_, fs.pressure(nonWettingPhaseIdx)));
588 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
590 else if (onInlet_(pos)) {
591 RateVector massRate(0.0);
593 massRate[contiNEqIdx] = -0.04;
596 values.setMassRate(massRate);
614 template <
class Context>
615 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
617 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
618 Scalar depth = this->boundingBoxMax()[1] - pos[1];
620 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
621 fs.setPressure(wettingPhaseIdx, 1e5);
624 fs.setSaturation(wettingPhaseIdx, Sw);
625 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
627 fs.setTemperature(temperature_);
629 typename FluidSystem::template ParameterCache<Scalar> paramCache;
630 paramCache.updatePhase(fs, wettingPhaseIdx);
631 Scalar densityW = FluidSystem::density(fs, paramCache, wettingPhaseIdx);
634 Scalar pw = 1e5 - densityW * this->gravity()[1] * depth;
637 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
638 Scalar pC[numPhases];
639 MaterialLaw::capillaryPressures(pC, matParams, fs);
642 fs.setPressure(wettingPhaseIdx, pw);
643 fs.setPressure(nonWettingPhaseIdx, pw + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]));
646 values.assignNaive(fs);
655 template <
class Context>
657 const Context& context OPM_UNUSED,
658 unsigned spaceIdx OPM_UNUSED,
659 unsigned timeIdx OPM_UNUSED)
const
660 { rate = Scalar(0.0); }
665 bool isInLens_(
const GlobalPosition& pos)
const
667 for (
unsigned i = 0; i < dim; ++i) {
668 if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i]
675 bool onLeftBoundary_(
const GlobalPosition& pos)
const
676 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
678 bool onRightBoundary_(
const GlobalPosition& pos)
const
679 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
681 bool onLowerBoundary_(
const GlobalPosition& pos)
const
682 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
684 bool onUpperBoundary_(
const GlobalPosition& pos)
const
685 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
687 bool onInlet_(
const GlobalPosition& pos)
const
689 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
690 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
691 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
694 GlobalPosition lensLowerLeft_;
695 GlobalPosition lensUpperRight_;
699 MaterialLawParams lensMaterialParams_;
700 MaterialLawParams outerMaterialParams_;
Soil contamination problem where DNAPL infiltrates a fully water saturated medium.
Definition: lensproblem.hh:266
static void registerParameters()
Definition: lensproblem.hh:362
void beginTimeStep()
Definition: lensproblem.hh:495
static std::string briefDescription()
Definition: lensproblem.hh:392
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:615
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:450
LensProblem(Simulator &simulator)
Definition: lensproblem.hh:310
void finishInit()
Definition: lensproblem.hh:317
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:534
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:427
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: lensproblem.hh:656
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: lensproblem.hh:441
std::string name() const
Definition: lensproblem.hh:479
void endTimeStep()
Definition: lensproblem.hh:507
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: lensproblem.hh:464
void beginIteration()
Definition: lensproblem.hh:501
Definition: groundwaterproblem.hh:61
Definition: groundwaterproblem.hh:63
Definition: groundwaterproblem.hh:65
Definition: groundwaterproblem.hh:67
Definition: groundwaterproblem.hh:69
Definition: groundwaterproblem.hh:71
Definition: lensproblem.hh:63