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>
54template <
class TypeTag>
58namespace Opm::Properties {
62struct LensBaseProblem {
using InheritsFrom = std::tuple<StructuredGridVanguard>; };
66template<
class TypeTag,
class MyTypeTag>
68template<
class TypeTag,
class MyTypeTag>
69struct LensLowerLeftY {
using type = UndefinedProperty; };
70template<
class TypeTag,
class MyTypeTag>
71struct LensLowerLeftZ {
using type = UndefinedProperty; };
72template<
class TypeTag,
class MyTypeTag>
73struct LensUpperRightX {
using type = UndefinedProperty; };
74template<
class TypeTag,
class MyTypeTag>
75struct LensUpperRightY {
using type = UndefinedProperty; };
76template<
class TypeTag,
class MyTypeTag>
77struct LensUpperRightZ {
using type = UndefinedProperty; };
80template<
class TypeTag>
84template<
class TypeTag>
85struct Grid<TypeTag, TTag::LensBaseProblem> {
using type = Dune::YaspGrid<2>; };
88template<
class TypeTag>
89struct WettingPhase<TypeTag, TTag::LensBaseProblem>
92 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
95 using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
99template<
class TypeTag>
100struct NonwettingPhase<TypeTag, TTag::LensBaseProblem>
103 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
106 using type = Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> >;
110template<
class TypeTag>
111struct 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>;
133template<
class TypeTag>
134struct NewtonWriteConvergence<TypeTag, TTag::LensBaseProblem> {
static constexpr bool value =
false; };
137template<
class TypeTag>
138struct NumericDifferenceMethod<TypeTag, TTag::LensBaseProblem> {
static constexpr int value = +1; };
141template<
class TypeTag>
142struct EnableGravity<TypeTag, TTag::LensBaseProblem> {
static constexpr bool value =
true; };
145template<
class TypeTag>
148 using type = GetPropType<TypeTag, Scalar>;
149 static constexpr type value = 1.0;
151template<
class TypeTag>
154 using type = GetPropType<TypeTag, Scalar>;
155 static constexpr type value = 2.0;
157template<
class TypeTag>
160 using type = GetPropType<TypeTag, Scalar>;
161 static constexpr type value = 0.0;
163template<
class TypeTag>
166 using type = GetPropType<TypeTag, Scalar>;
167 static constexpr type value = 4.0;
169template<
class TypeTag>
172 using type = GetPropType<TypeTag, Scalar>;
173 static constexpr type value = 3.0;
175template<
class TypeTag>
178 using type = GetPropType<TypeTag, Scalar>;
179 static constexpr type value = 1.0;
182template<
class TypeTag>
183struct DomainSizeX<TypeTag, TTag::LensBaseProblem>
185 using type = GetPropType<TypeTag, Scalar>;
186 static constexpr type value = 6.0;
188template<
class TypeTag>
189struct DomainSizeY<TypeTag, TTag::LensBaseProblem>
191 using type = GetPropType<TypeTag, Scalar>;
192 static constexpr type value = 4.0;
194template<
class TypeTag>
195struct DomainSizeZ<TypeTag, TTag::LensBaseProblem>
197 using type = GetPropType<TypeTag, Scalar>;
198 static constexpr type value = 1.0;
201template<
class TypeTag>
202struct CellsX<TypeTag, TTag::LensBaseProblem> {
static constexpr unsigned value = 48; };
203template<
class TypeTag>
204struct CellsY<TypeTag, TTag::LensBaseProblem> {
static constexpr unsigned value = 32; };
205template<
class TypeTag>
206struct CellsZ<TypeTag, TTag::LensBaseProblem> {
static constexpr unsigned value = 16; };
209template<
class TypeTag>
210struct EndTime<TypeTag, TTag::LensBaseProblem>
212 using type = GetPropType<TypeTag, Scalar>;
213 static constexpr type value = 30e3;
217template<
class TypeTag>
218struct InitialTimeStepSize<TypeTag, TTag::LensBaseProblem>
220 using type = GetPropType<TypeTag, Scalar>;
221 static constexpr type value = 250;
225template<
class TypeTag>
226struct VtkWriteIntrinsicPermeabilities<TypeTag, TTag::LensBaseProblem> {
static constexpr bool value =
true; };
229template<
class TypeTag>
230struct EnableStorageCache<TypeTag, TTag::LensBaseProblem> {
static constexpr bool value =
true; };
233template<
class TypeTag>
234struct EnableIntensiveQuantityCache<TypeTag, TTag::LensBaseProblem> {
static constexpr bool value =
true; };
263template <
class TypeTag>
264class 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] = Parameters::get<TypeTag, Properties::LensLowerLeftX>();
325 lensLowerLeft_[1] = Parameters::get<TypeTag, Properties::LensLowerLeftY>();
326 lensUpperRight_[0] = Parameters::get<TypeTag, Properties::LensUpperRightX>();
327 lensUpperRight_[1] = Parameters::get<TypeTag, Properties::LensUpperRightY>();
330 lensLowerLeft_[2] = Parameters::get<TypeTag, Properties::LensLowerLeftZ>();
331 lensUpperRight_[2] = Parameters::get<TypeTag, Properties::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 Parameters::registerParam<TypeTag, Properties::LensLowerLeftX>
366 (
"The x-coordinate of the lens' lower-left corner [m].");
367 Parameters::registerParam<TypeTag, Properties::LensLowerLeftY>
368 (
"The y-coordinate of the lens' lower-left corner [m].");
369 Parameters::registerParam<TypeTag, Properties::LensUpperRightX>
370 (
"The x-coordinate of the lens' upper-right corner [m].");
371 Parameters::registerParam<TypeTag, Properties::LensUpperRightY>
372 (
"The y-coordinate of the lens' upper-right corner [m].");
375 Parameters::registerParam<TypeTag, Properties::LensLowerLeftZ>
376 (
"The z-coordinate of the lens' lower-left corner [m].");
377 Parameters::registerParam<TypeTag, Properties::LensUpperRightZ>
378 (
"The z-coordinate of the lens' upper-right corner [m].");
387 std::string thermal =
"isothermal";
388 bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
390 thermal =
"non-isothermal";
392 std::string deriv =
"finite difference";
393 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
394 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
396 deriv =
"automatic differentiation";
398 std::string disc =
"vertex centered finite volume";
399 using D = GetPropType<TypeTag, Properties::Discretization>;
400 bool useEcfv = std::is_same<D, Opm::EcfvDiscretization<TypeTag>>::value;
402 disc =
"element centered finite volume";
404 return std::string(
"")+
405 "Ground remediation problem where a dense oil infiltrates "+
406 "an aquifer with an embedded low-permability lens. " +
407 "This is the binary for the "+thermal+
" variant using "+deriv+
408 "and the "+disc+
" discretization";
419 template <
class Context>
421 unsigned timeIdx)
const
423 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
425 if (isInLens_(globalPos))
433 template <
class Context>
442 template <
class Context>
444 unsigned spaceIdx,
unsigned timeIdx)
const
446 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
448 if (isInLens_(globalPos))
449 return lensMaterialParams_;
450 return outerMaterialParams_;
456 template <
class Context>
460 {
return temperature_; }
474 using LLS = GetPropType<TypeTag, Properties::LocalLinearizerSplice>;
476 bool useAutoDiff = std::is_same<LLS, Properties::TTag::AutoDiffLocalLinearizer>::value;
478 using FM = GetPropType<TypeTag, Properties::FluxModule>;
479 bool useTrans = std::is_same<FM, Opm::TransFluxModule<TypeTag>>::value;
481 std::ostringstream oss;
482 oss <<
"lens_" << Model::name()
483 <<
"_" << Model::discretizationName()
484 <<
"_" << (useAutoDiff?
"ad":
"fd");
513 this->model().globalStorage(storage);
516 if (this->gridView().comm().rank() == 0) {
517 std::cout <<
"Storage: " << storage << std::endl << std::flush;
532 template <
class Context>
534 const Context& context,
536 unsigned timeIdx)
const
538 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
540 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
542 Scalar densityW = WettingPhase::density(temperature_, Scalar(1e5));
543 Scalar densityN = NonwettingPhase::density(temperature_, Scalar(1e5));
545 Scalar T =
temperature(context, spaceIdx, timeIdx);
549 if (onLeftBoundary_(pos)) {
550 Scalar height = this->boundingBoxMax()[1] - this->boundingBoxMin()[1];
551 Scalar depth = this->boundingBoxMax()[1] - pos[1];
552 Scalar alpha = (1 + 1.5 / height);
555 pw = 1e5 - alpha * densityW * this->gravity()[1] * depth;
559 Scalar depth = this->boundingBoxMax()[1] - pos[1];
562 pw = 1e5 - densityW * this->gravity()[1] * depth;
567 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
569 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
571 fs.setSaturation(wettingPhaseIdx, Sw);
572 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
573 fs.setTemperature(T);
575 Scalar pC[numPhases];
576 MaterialLaw::capillaryPressures(pC, matParams, fs);
577 fs.setPressure(wettingPhaseIdx, pw);
578 fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]);
580 fs.setDensity(wettingPhaseIdx, densityW);
581 fs.setDensity(nonWettingPhaseIdx, densityN);
583 fs.setViscosity(wettingPhaseIdx, WettingPhase::viscosity(temperature_, fs.pressure(wettingPhaseIdx)));
584 fs.setViscosity(nonWettingPhaseIdx, NonwettingPhase::viscosity(temperature_, fs.pressure(nonWettingPhaseIdx)));
587 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
589 else if (onInlet_(pos)) {
590 RateVector massRate(0.0);
592 massRate[contiNEqIdx] = -0.04;
595 values.setMassRate(massRate);
613 template <
class Context>
614 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
616 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
617 Scalar depth = this->boundingBoxMax()[1] - pos[1];
619 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
620 fs.setPressure(wettingPhaseIdx, 1e5);
623 fs.setSaturation(wettingPhaseIdx, Sw);
624 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
626 fs.setTemperature(temperature_);
628 typename FluidSystem::template ParameterCache<Scalar> paramCache;
629 paramCache.updatePhase(fs, wettingPhaseIdx);
630 Scalar densityW = FluidSystem::density(fs, paramCache, wettingPhaseIdx);
633 Scalar pw = 1e5 - densityW * this->gravity()[1] * depth;
636 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
637 Scalar pC[numPhases];
638 MaterialLaw::capillaryPressures(pC, matParams, fs);
641 fs.setPressure(wettingPhaseIdx, pw);
642 fs.setPressure(nonWettingPhaseIdx, pw + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]));
645 values.assignNaive(fs);
654 template <
class Context>
659 { rate = Scalar(0.0); }
664 bool isInLens_(
const GlobalPosition& pos)
const
666 for (
unsigned i = 0; i < dim; ++i) {
667 if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i]
674 bool onLeftBoundary_(
const GlobalPosition& pos)
const
675 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
677 bool onRightBoundary_(
const GlobalPosition& pos)
const
678 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
680 bool onLowerBoundary_(
const GlobalPosition& pos)
const
681 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
683 bool onUpperBoundary_(
const GlobalPosition& pos)
const
684 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
686 bool onInlet_(
const GlobalPosition& pos)
const
688 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
689 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
690 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
693 GlobalPosition lensLowerLeft_;
694 GlobalPosition lensUpperRight_;
698 MaterialLawParams lensMaterialParams_;
699 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:457
static void registerParameters()
Definition lensproblem.hh:361
void beginTimeStep()
Definition lensproblem.hh:494
static std::string briefDescription()
Definition lensproblem.hh:385
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition lensproblem.hh:443
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition lensproblem.hh:614
Scalar porosity(const Context &, unsigned, unsigned) const
Definition lensproblem.hh:434
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:533
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition lensproblem.hh:420
std::string name() const
Definition lensproblem.hh:472
void endTimeStep()
Definition lensproblem.hh:506
void beginIteration()
Definition lensproblem.hh:500
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition lensproblem.hh:655
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