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>
38 #include <opm/material/common/Unused.hpp>
40 #include <dune/grid/yaspgrid.hh>
41 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
43 #include <dune/common/version.hh>
44 #include <dune/common/fvector.hh>
45 #include <dune/common/fmatrix.hh>
51 template <
class TypeTag>
52 class InfiltrationProblem;
55 namespace Opm::Properties {
62 template<
class TypeTag>
63 struct Grid<TypeTag, TTag::InfiltrationBaseProblem> {
using type = Dune::YaspGrid<2>; };
66 template<
class TypeTag>
70 template<
class TypeTag>
71 struct FluidSystem<TypeTag, TTag::InfiltrationBaseProblem>
72 {
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
75 template<
class TypeTag>
76 struct EnableGravity<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
bool value =
true; };
79 template<
class TypeTag>
80 struct NewtonWriteConvergence<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
bool value =
false; };
83 template<
class TypeTag>
84 struct NumericDifferenceMethod<TypeTag, TTag::InfiltrationBaseProblem> {
static constexpr
int value = 1; };
87 template<
class TypeTag>
88 struct MaterialLaw<TypeTag, TTag::InfiltrationBaseProblem>
91 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
92 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
94 using Traits= Opm::ThreePhaseMaterialTraits<
96 FluidSystem::waterPhaseIdx,
97 FluidSystem::naplPhaseIdx,
98 FluidSystem::gasPhaseIdx>;
101 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
105 template<
class TypeTag>
106 struct EndTime<TypeTag, TTag::InfiltrationBaseProblem>
108 using type = GetPropType<TypeTag, Scalar>;
109 static constexpr type value = 6e3;
113 template<
class TypeTag>
114 struct InitialTimeStepSize<TypeTag, TTag::InfiltrationBaseProblem>
116 using type = GetPropType<TypeTag, Scalar>;
117 static constexpr type value = 60;
121 template<
class TypeTag>
122 struct GridFile<TypeTag, TTag::InfiltrationBaseProblem>
123 {
static constexpr
auto value =
"./data/infiltration_50x3.dgf"; };
150 template <
class TypeTag>
153 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
155 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
156 using GridView = GetPropType<TypeTag, Properties::GridView>;
157 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
158 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
159 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
160 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
161 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
162 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
163 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
164 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
165 using Model = GetPropType<TypeTag, Properties::Model>;
168 using Indices = GetPropType<TypeTag, Properties::Indices>;
171 conti0EqIdx = Indices::conti0EqIdx,
174 numPhases = FluidSystem::numPhases,
177 NAPLIdx = FluidSystem::NAPLIdx,
178 H2OIdx = FluidSystem::H2OIdx,
179 airIdx = FluidSystem::airIdx,
182 waterPhaseIdx = FluidSystem::waterPhaseIdx,
183 gasPhaseIdx = FluidSystem::gasPhaseIdx,
184 naplPhaseIdx = FluidSystem::naplPhaseIdx,
187 dim = GridView::dimension,
188 dimWorld = GridView::dimensionworld
191 using CoordScalar =
typename GridView::ctype;
192 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
193 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
200 : ParentType(simulator)
209 ParentType::finishInit();
211 temperature_ = 273.15 + 10.0;
212 FluidSystem::init(temperature_ - 1,
220 fineK_ = this->toDimMatrix_(1e-11);
221 coarseK_ = this->toDimMatrix_(1e-11);
227 materialParams_.setSwr(0.12);
228 materialParams_.setSwrx(0.12);
229 materialParams_.setSnr(0.07);
230 materialParams_.setSgr(0.03);
233 materialParams_.setVgAlpha(0.0005);
234 materialParams_.setVgN(4.);
235 materialParams_.setkrRegardsSnr(
false);
237 materialParams_.finalize();
238 materialParams_.checkDefined();
259 std::ostringstream oss;
260 oss <<
"infiltration_" << Model::name();
270 this->model().checkConservativeness();
274 this->model().globalStorage(storage);
277 if (this->gridView().comm().rank() == 0) {
278 std::cout <<
"Storage: " << storage << std::endl << std::flush;
286 template <
class Context>
288 unsigned spaceIdx OPM_UNUSED,
289 unsigned timeIdx OPM_UNUSED)
const
290 {
return temperature_; }
295 template <
class Context>
299 unsigned timeIdx)
const
301 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
302 if (isFineMaterial_(pos))
310 template <
class Context>
312 unsigned spaceIdx OPM_UNUSED,
313 unsigned timeIdx OPM_UNUSED)
const
314 {
return porosity_; }
319 template <
class Context>
320 const MaterialLawParams&
322 unsigned spaceIdx OPM_UNUSED,
323 unsigned timeIdx OPM_UNUSED)
const
324 {
return materialParams_; }
336 template <
class Context>
338 const Context& context,
340 unsigned timeIdx)
const
342 const auto& pos = context.pos(spaceIdx, timeIdx);
344 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
345 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
347 initialFluidState_(fs, context, spaceIdx, timeIdx);
349 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
351 else if (onInlet_(pos)) {
352 RateVector molarRate(0.0);
353 molarRate[conti0EqIdx + NAPLIdx] = -0.001;
355 values.setMolarRate(molarRate);
356 Opm::Valgrind::CheckDefined(values);
372 template <
class Context>
374 const Context& context,
376 unsigned timeIdx)
const
378 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
380 initialFluidState_(fs, context, spaceIdx, timeIdx);
383 values.assignMassConservative(fs, matParams,
true);
384 Opm::Valgrind::CheckDefined(values);
393 template <
class Context>
395 const Context& context OPM_UNUSED,
396 unsigned spaceIdx OPM_UNUSED,
397 unsigned timeIdx OPM_UNUSED)
const
398 { rate = Scalar(0.0); }
403 bool onLeftBoundary_(
const GlobalPosition& pos)
const
404 {
return pos[0] < eps_; }
406 bool onRightBoundary_(
const GlobalPosition& pos)
const
407 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
409 bool onLowerBoundary_(
const GlobalPosition& pos)
const
410 {
return pos[1] < eps_; }
412 bool onUpperBoundary_(
const GlobalPosition& pos)
const
413 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
415 bool onInlet_(
const GlobalPosition& pos)
const
416 {
return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
418 template <
class Flu
idState,
class Context>
419 void initialFluidState_(FluidState& fs,
const Context& context,
420 unsigned spaceIdx,
unsigned timeIdx)
const
422 const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
426 Scalar densityW = 1000.0;
427 Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
433 Scalar Sw = matParams.Swr();
434 Scalar Swr = matParams.Swr();
435 Scalar Sgr = matParams.Sgr();
442 Opm::Valgrind::CheckDefined(Sw);
443 Opm::Valgrind::CheckDefined(Sg);
445 fs.setSaturation(waterPhaseIdx, Sw);
446 fs.setSaturation(gasPhaseIdx, Sg);
447 fs.setSaturation(naplPhaseIdx, 0);
450 fs.setTemperature(temperature_);
453 Scalar pcAll[numPhases];
455 if (onLeftBoundary_(pos))
457 MaterialLaw::capillaryPressures(pcAll, matParams, fs);
458 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
459 fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
462 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
463 fs.setMoleFraction(gasPhaseIdx, airIdx,
464 1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
465 fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
467 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
468 typename FluidSystem::template ParameterCache<Scalar> paramCache;
469 CFRP::solve(fs, paramCache, gasPhaseIdx,
473 fs.setMoleFraction(waterPhaseIdx, H2OIdx,
474 1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
477 bool isFineMaterial_(
const GlobalPosition& pos)
const
478 {
return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
485 MaterialLawParams materialParams_;
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:152
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:297
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:321
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:311
std::string name() const
Definition: infiltrationproblem.hh:257
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:337
void endTimeStep()
Definition: infiltrationproblem.hh:267
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:199
bool shouldWriteRestartFile() const
Definition: infiltrationproblem.hh:251
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:394
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:373
void finishInit()
Definition: infiltrationproblem.hh:207
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:287
Definition: infiltrationproblem.hh:58