27 #ifndef EWOMS_INFILTRATION_PROBLEM_HH
28 #define EWOMS_INFILTRATION_PROBLEM_HH
30 #include <opm/models/pvs/pvsproperties.hh>
32 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
33 #include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>
34 #include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
35 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
36 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
37 #include <opm/material/common/Valgrind.hpp>
39 #include <dune/grid/yaspgrid.hh>
40 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
42 #include <dune/common/version.hh>
43 #include <dune/common/fvector.hh>
44 #include <dune/common/fmatrix.hh>
50 template <
class TypeTag>
51 class InfiltrationProblem;
54 namespace Opm::Properties {
61 template<
class TypeTag>
62 struct Grid<TypeTag, TTag::InfiltrationBaseProblem> {
using type = Dune::YaspGrid<2>; };
65 template<
class TypeTag>
69 template<
class TypeTag>
70 struct FluidSystem<TypeTag, TTag::InfiltrationBaseProblem>
71 {
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
74 template<
class TypeTag>
75 struct EnableGravity<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
bool value =
true; };
78 template<
class TypeTag>
79 struct NewtonWriteConvergence<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
bool value =
false; };
82 template<
class TypeTag>
83 struct NumericDifferenceMethod<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
int value = 1; };
86 template<
class TypeTag>
87 struct MaterialLaw<TypeTag, TTag::InfiltrationBaseProblem>
90 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
93 using Traits= Opm::ThreePhaseMaterialTraits<
95 FluidSystem::waterPhaseIdx,
96 FluidSystem::naplPhaseIdx,
97 FluidSystem::gasPhaseIdx>;
100 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
104 template<
class TypeTag>
105 struct EndTime<TypeTag, TTag::InfiltrationBaseProblem>
107 using type = GetPropType<TypeTag, Scalar>;
108 static constexpr type value = 6e3;
112 template<
class TypeTag>
113 struct InitialTimeStepSize<TypeTag, TTag::InfiltrationBaseProblem>
115 using type = GetPropType<TypeTag, Scalar>;
116 static constexpr type value = 60;
120 template<
class TypeTag>
121 struct GridFile<TypeTag, TTag::InfiltrationBaseProblem>
122 {
static constexpr
auto value =
"./data/infiltration_50x3.dgf"; };
149 template <
class TypeTag>
152 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
154 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
155 using GridView = GetPropType<TypeTag, Properties::GridView>;
156 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
157 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
158 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
159 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
160 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
161 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
162 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
163 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
164 using Model = GetPropType<TypeTag, Properties::Model>;
167 using Indices = GetPropType<TypeTag, Properties::Indices>;
170 conti0EqIdx = Indices::conti0EqIdx,
173 numPhases = FluidSystem::numPhases,
176 NAPLIdx = FluidSystem::NAPLIdx,
177 H2OIdx = FluidSystem::H2OIdx,
178 airIdx = FluidSystem::airIdx,
181 waterPhaseIdx = FluidSystem::waterPhaseIdx,
182 gasPhaseIdx = FluidSystem::gasPhaseIdx,
183 naplPhaseIdx = FluidSystem::naplPhaseIdx,
186 dim = GridView::dimension,
187 dimWorld = GridView::dimensionworld
190 using CoordScalar =
typename GridView::ctype;
191 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
192 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
199 : ParentType(simulator)
208 ParentType::finishInit();
210 temperature_ = 273.15 + 10.0;
211 FluidSystem::init(temperature_ - 1,
219 fineK_ = this->toDimMatrix_(1e-11);
220 coarseK_ = this->toDimMatrix_(1e-11);
226 materialParams_.setSwr(0.12);
227 materialParams_.setSwrx(0.12);
228 materialParams_.setSnr(0.07);
229 materialParams_.setSgr(0.03);
232 materialParams_.setVgAlpha(0.0005);
233 materialParams_.setVgN(4.);
234 materialParams_.setkrRegardsSnr(
false);
236 materialParams_.finalize();
237 materialParams_.checkDefined();
258 std::ostringstream oss;
259 oss <<
"infiltration_" << Model::name();
269 this->model().checkConservativeness();
273 this->model().globalStorage(storage);
276 if (this->gridView().comm().rank() == 0) {
277 std::cout <<
"Storage: " << storage << std::endl << std::flush;
285 template <
class Context>
289 {
return temperature_; }
294 template <
class Context>
298 unsigned timeIdx)
const
300 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
301 if (isFineMaterial_(pos))
309 template <
class Context>
313 {
return porosity_; }
318 template <
class Context>
319 const MaterialLawParams&
323 {
return materialParams_; }
335 template <
class Context>
337 const Context& context,
339 unsigned timeIdx)
const
341 const auto& pos = context.pos(spaceIdx, timeIdx);
343 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
344 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
346 initialFluidState_(fs, context, spaceIdx, timeIdx);
348 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
350 else if (onInlet_(pos)) {
351 RateVector molarRate(0.0);
352 molarRate[conti0EqIdx + NAPLIdx] = -0.001;
354 values.setMolarRate(molarRate);
355 Opm::Valgrind::CheckDefined(values);
371 template <
class Context>
373 const Context& context,
375 unsigned timeIdx)
const
377 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
379 initialFluidState_(fs, context, spaceIdx, timeIdx);
382 values.assignMassConservative(fs, matParams,
true);
383 Opm::Valgrind::CheckDefined(values);
392 template <
class Context>
397 { rate = Scalar(0.0); }
402 bool onLeftBoundary_(
const GlobalPosition& pos)
const
403 {
return pos[0] < eps_; }
405 bool onRightBoundary_(
const GlobalPosition& pos)
const
406 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
408 bool onLowerBoundary_(
const GlobalPosition& pos)
const
409 {
return pos[1] < eps_; }
411 bool onUpperBoundary_(
const GlobalPosition& pos)
const
412 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
414 bool onInlet_(
const GlobalPosition& pos)
const
415 {
return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
417 template <
class Flu
idState,
class Context>
418 void initialFluidState_(FluidState& fs,
const Context& context,
419 unsigned spaceIdx,
unsigned timeIdx)
const
421 const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
425 Scalar densityW = 1000.0;
426 Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
432 Scalar Sw = matParams.Swr();
433 Scalar Swr = matParams.Swr();
434 Scalar Sgr = matParams.Sgr();
441 Opm::Valgrind::CheckDefined(Sw);
442 Opm::Valgrind::CheckDefined(Sg);
444 fs.setSaturation(waterPhaseIdx, Sw);
445 fs.setSaturation(gasPhaseIdx, Sg);
446 fs.setSaturation(naplPhaseIdx, 0);
449 fs.setTemperature(temperature_);
452 Scalar pcAll[numPhases];
454 if (onLeftBoundary_(pos))
456 MaterialLaw::capillaryPressures(pcAll, matParams, fs);
457 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
458 fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
461 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
462 fs.setMoleFraction(gasPhaseIdx, airIdx,
463 1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
464 fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
466 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
467 typename FluidSystem::template ParameterCache<Scalar> paramCache;
468 CFRP::solve(fs, paramCache, gasPhaseIdx,
472 fs.setMoleFraction(waterPhaseIdx, H2OIdx,
473 1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
476 bool isFineMaterial_(
const GlobalPosition& pos)
const
477 {
return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
484 MaterialLawParams materialParams_;
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:151
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:393
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:296
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:286
std::string name() const
Definition: infiltrationproblem.hh:256
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:310
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:336
void endTimeStep()
Definition: infiltrationproblem.hh:266
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:198
bool shouldWriteRestartFile() const
Definition: infiltrationproblem.hh:250
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:320
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:372
void finishInit()
Definition: infiltrationproblem.hh:206
Definition: infiltrationproblem.hh:57