27#ifndef EWOMS_INFILTRATION_PROBLEM_HH
28#define EWOMS_INFILTRATION_PROBLEM_HH
30#include <dune/common/fmatrix.hh>
31#include <dune/common/fvector.hh>
32#include <dune/common/version.hh>
34#include <dune/grid/yaspgrid.hh>
35#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
37#include <opm/material/common/Valgrind.hpp>
38#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
39#include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
40#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
41#include <opm/material/fluidstates/CompositionalFluidState.hpp>
42#include <opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp>
54template <
class TypeTag>
55class InfiltrationProblem;
65template<
class TypeTag>
66struct Grid<TypeTag, TTag::InfiltrationBaseProblem> {
using type = Dune::YaspGrid<2>; };
69template<
class TypeTag>
73template<
class TypeTag>
75{
using type = Opm::H2OAirMesityleneFluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
78template<
class TypeTag>
85 using Traits= Opm::ThreePhaseMaterialTraits<
87 FluidSystem::waterPhaseIdx,
88 FluidSystem::naplPhaseIdx,
89 FluidSystem::gasPhaseIdx>;
92 using type = Opm::ThreePhaseParkerVanGenuchten<Traits>;
120template <
class TypeTag>
141 conti0EqIdx = Indices::conti0EqIdx,
144 numPhases = FluidSystem::numPhases,
147 NAPLIdx = FluidSystem::NAPLIdx,
148 H2OIdx = FluidSystem::H2OIdx,
149 airIdx = FluidSystem::airIdx,
152 waterPhaseIdx = FluidSystem::waterPhaseIdx,
153 gasPhaseIdx = FluidSystem::gasPhaseIdx,
154 naplPhaseIdx = FluidSystem::naplPhaseIdx,
157 dim = GridView::dimension,
158 dimWorld = GridView::dimensionworld
161 using CoordScalar =
typename GridView::ctype;
162 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
163 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
170 : ParentType(simulator)
179 ParentType::finishInit();
181 temperature_ = 273.15 + 10.0;
182 FluidSystem::init(temperature_ - 1,
190 fineK_ = this->toDimMatrix_(1e-11);
191 coarseK_ = this->toDimMatrix_(1e-11);
197 materialParams_.setSwr(0.12);
198 materialParams_.setSwrx(0.12);
199 materialParams_.setSnr(0.07);
200 materialParams_.setSgr(0.03);
203 materialParams_.setVgAlpha(0.0005);
204 materialParams_.setVgN(4.);
205 materialParams_.setkrRegardsSnr(
false);
207 materialParams_.finalize();
208 materialParams_.checkDefined();
216 ParentType::registerParameters();
218 Parameters::SetDefault<Parameters::GridFile>(
"./data/infiltration_50x3.dgf");
219 Parameters::SetDefault<Parameters::NumericDifferenceMethod>(1);
220 Parameters::SetDefault<Parameters::EndTime<Scalar>>(6e3);
221 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(60.0);
222 Parameters::SetDefault<Parameters::EnableGravity>(
true);
243 std::ostringstream oss;
244 oss <<
"infiltration_" << Model::name();
254 this->model().checkConservativeness();
258 this->model().globalStorage(storage);
261 if (this->gridView().comm().rank() == 0) {
262 std::cout <<
"Storage: " << storage << std::endl << std::flush;
270 template <
class Context>
274 {
return temperature_; }
279 template <
class Context>
283 unsigned timeIdx)
const
285 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
286 if (isFineMaterial_(pos))
294 template <
class Context>
298 {
return porosity_; }
303 template <
class Context>
304 const MaterialLawParams&
308 {
return materialParams_; }
320 template <
class Context>
322 const Context& context,
324 unsigned timeIdx)
const
326 const auto& pos = context.pos(spaceIdx, timeIdx);
328 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
329 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
331 initialFluidState_(fs, context, spaceIdx, timeIdx);
333 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
335 else if (onInlet_(pos)) {
336 RateVector molarRate(0.0);
337 molarRate[conti0EqIdx + NAPLIdx] = -0.001;
339 values.setMolarRate(molarRate);
340 Opm::Valgrind::CheckDefined(values);
356 template <
class Context>
358 const Context& context,
360 unsigned timeIdx)
const
362 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
364 initialFluidState_(fs, context, spaceIdx, timeIdx);
367 values.assignMassConservative(fs, matParams,
true);
368 Opm::Valgrind::CheckDefined(values);
377 template <
class Context>
382 { rate = Scalar(0.0); }
387 bool onLeftBoundary_(
const GlobalPosition& pos)
const
388 {
return pos[0] < eps_; }
390 bool onRightBoundary_(
const GlobalPosition& pos)
const
391 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
393 bool onLowerBoundary_(
const GlobalPosition& pos)
const
394 {
return pos[1] < eps_; }
396 bool onUpperBoundary_(
const GlobalPosition& pos)
const
397 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
399 bool onInlet_(
const GlobalPosition& pos)
const
400 {
return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
402 template <
class Flu
idState,
class Context>
403 void initialFluidState_(FluidState& fs,
const Context& context,
404 unsigned spaceIdx,
unsigned timeIdx)
const
406 const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
410 Scalar densityW = 1000.0;
411 Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
417 Scalar Sw = matParams.Swr();
418 Scalar Swr = matParams.Swr();
419 Scalar Sgr = matParams.Sgr();
426 Opm::Valgrind::CheckDefined(Sw);
427 Opm::Valgrind::CheckDefined(Sg);
429 fs.setSaturation(waterPhaseIdx, Sw);
430 fs.setSaturation(gasPhaseIdx, Sg);
431 fs.setSaturation(naplPhaseIdx, 0);
434 fs.setTemperature(temperature_);
437 Scalar pcAll[numPhases];
439 if (onLeftBoundary_(pos))
441 MaterialLaw::capillaryPressures(pcAll, matParams, fs);
442 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
443 fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
446 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
447 fs.setMoleFraction(gasPhaseIdx, airIdx,
448 1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
449 fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
451 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
452 typename FluidSystem::template ParameterCache<Scalar> paramCache;
453 CFRP::solve(fs, paramCache, gasPhaseIdx,
457 fs.setMoleFraction(waterPhaseIdx, H2OIdx,
458 1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
461 bool isFineMaterial_(
const GlobalPosition& pos)
const
462 {
return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
469 MaterialLawParams materialParams_;
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:122
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: infiltrationproblem.hh:378
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:271
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:281
std::string name() const
The problem name.
Definition: infiltrationproblem.hh:241
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:295
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: infiltrationproblem.hh:321
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition: infiltrationproblem.hh:305
void endTimeStep()
Called by the simulator after each time integration.
Definition: infiltrationproblem.hh:251
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:169
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: infiltrationproblem.hh:235
static void registerParameters()
Definition: infiltrationproblem.hh:214
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: infiltrationproblem.hh:357
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: infiltrationproblem.hh:177
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
Declares the properties required for the compositional multi-phase primary variable switching model.
Opm::H2OAirMesityleneFluidSystem< GetPropType< TypeTag, Properties::Scalar > > type
Definition: infiltrationproblem.hh:75
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:69
Dune::YaspGrid< 2 > type
Definition: infiltrationproblem.hh:66
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::ThreePhaseParkerVanGenuchten< Traits > type
Definition: infiltrationproblem.hh:92
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Definition: infiltrationproblem.hh:61