27 #ifndef EWOMS_INFILTRATION_PROBLEM_HH
28 #define EWOMS_INFILTRATION_PROBLEM_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/heatconduction/Somerton.hpp>
38 #include <opm/common/Valgrind.hpp>
39 #include <opm/common/Unused.hpp>
41 #include <dune/grid/yaspgrid.hh>
42 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
44 #include <dune/common/version.hh>
45 #include <dune/common/fvector.hh>
46 #include <dune/common/fmatrix.hh>
52 template <
class TypeTag>
57 namespace Properties {
61 SET_TYPE_PROP(InfiltrationBaseProblem, Grid, Dune::YaspGrid<2>);
69 InfiltrationBaseProblem, FluidSystem,
70 Opm::FluidSystems::H2OAirMesitylene<
typename GET_PROP_TYPE(TypeTag, Scalar)>);
76 SET_BOOL_PROP(InfiltrationBaseProblem, NewtonWriteConvergence,
false);
79 SET_INT_PROP(InfiltrationBaseProblem, NumericDifferenceMethod, 1);
82 SET_PROP(InfiltrationBaseProblem, MaterialLaw)
86 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
88 typedef Opm::ThreePhaseMaterialTraits<
90 FluidSystem::waterPhaseIdx,
91 FluidSystem::naplPhaseIdx,
92 FluidSystem::gasPhaseIdx> Traits;
95 typedef Opm::ThreePhaseParkerVanGenuchten<Traits> type;
99 SET_PROP(InfiltrationBaseProblem, HeatConductionLaw)
103 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
107 typedef Opm::Somerton<FluidSystem, Scalar> type;
118 "./data/infiltration_50x3.dgf");
145 template <
class TypeTag>
148 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
152 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
153 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
154 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
156 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
157 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
159 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
166 conti0EqIdx = Indices::conti0EqIdx,
169 numPhases = FluidSystem::numPhases,
172 NAPLIdx = FluidSystem::NAPLIdx,
173 H2OIdx = FluidSystem::H2OIdx,
174 airIdx = FluidSystem::airIdx,
177 waterPhaseIdx = FluidSystem::waterPhaseIdx,
178 gasPhaseIdx = FluidSystem::gasPhaseIdx,
179 naplPhaseIdx = FluidSystem::naplPhaseIdx,
182 dim = GridView::dimension,
183 dimWorld = GridView::dimensionworld
186 typedef typename GridView::ctype CoordScalar;
187 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
188 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
195 : ParentType(simulator)
204 ParentType::finishInit();
206 temperature_ = 273.15 + 10.0;
207 FluidSystem::init(temperature_ - 1,
215 fineK_ = this->toDimMatrix_(1e-11);
216 coarseK_ = this->toDimMatrix_(1e-11);
222 materialParams_.setSwr(0.12);
223 materialParams_.setSwrx(0.12);
224 materialParams_.setSnr(0.07);
225 materialParams_.setSgr(0.03);
228 materialParams_.setVgAlpha(0.0005);
229 materialParams_.setVgN(4.);
230 materialParams_.setkrRegardsSnr(
false);
232 materialParams_.finalize();
233 materialParams_.checkDefined();
254 std::ostringstream oss;
255 oss <<
"infiltration_" << Model::name();
265 this->model().checkConservativeness();
269 this->model().globalStorage(storage);
272 if (this->gridView().comm().rank() == 0) {
273 std::cout <<
"Storage: " << storage << std::endl << std::flush;
281 template <
class Context>
283 unsigned spaceIdx OPM_UNUSED,
284 unsigned timeIdx OPM_UNUSED)
const
285 {
return temperature_; }
290 template <
class Context>
294 unsigned timeIdx)
const
296 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
297 if (isFineMaterial_(pos))
305 template <
class Context>
307 unsigned spaceIdx OPM_UNUSED,
308 unsigned timeIdx OPM_UNUSED)
const
309 {
return porosity_; }
314 template <
class Context>
315 const MaterialLawParams&
317 unsigned spaceIdx OPM_UNUSED,
318 unsigned timeIdx OPM_UNUSED)
const
319 {
return materialParams_; }
326 template <
class Context>
328 unsigned spaceIdx OPM_UNUSED,
329 unsigned timeIdx OPM_UNUSED)
const
345 template <
class Context>
347 const Context& context,
349 unsigned timeIdx)
const
351 const auto& pos = context.pos(spaceIdx, timeIdx);
353 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
354 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
356 initialFluidState_(fs, context, spaceIdx, timeIdx);
358 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
360 else if (onInlet_(pos)) {
361 RateVector molarRate(0.0);
362 molarRate[conti0EqIdx + NAPLIdx] = -0.001;
364 values.setMolarRate(molarRate);
365 Opm::Valgrind::CheckDefined(values);
381 template <
class Context>
383 const Context& context,
385 unsigned timeIdx)
const
387 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
389 initialFluidState_(fs, context, spaceIdx, timeIdx);
391 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
392 values.assignMassConservative(fs, matParams,
true);
393 Opm::Valgrind::CheckDefined(values);
402 template <
class Context>
404 const Context& context OPM_UNUSED,
405 unsigned spaceIdx OPM_UNUSED,
406 unsigned timeIdx OPM_UNUSED)
const
407 { rate = Scalar(0.0); }
412 bool onLeftBoundary_(
const GlobalPosition& pos)
const
413 {
return pos[0] < eps_; }
415 bool onRightBoundary_(
const GlobalPosition& pos)
const
416 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
418 bool onLowerBoundary_(
const GlobalPosition& pos)
const
419 {
return pos[1] < eps_; }
421 bool onUpperBoundary_(
const GlobalPosition& pos)
const
422 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
424 bool onInlet_(
const GlobalPosition& pos)
const
425 {
return onUpperBoundary_(pos) && 50 < pos[0] && pos[0] < 75; }
427 template <
class Flu
idState,
class Context>
428 void initialFluidState_(FluidState& fs,
const Context& context,
429 unsigned spaceIdx,
unsigned timeIdx)
const
431 const GlobalPosition pos = context.pos(spaceIdx, timeIdx);
435 Scalar densityW = 1000.0;
436 Scalar pc = 9.81 * densityW * (y - (5 - 5e-4 * x));
441 const auto& matParams = materialLawParams(context, spaceIdx, timeIdx);
442 Scalar Sw = matParams.Swr();
443 Scalar Swr = matParams.Swr();
444 Scalar Sgr = matParams.Sgr();
451 Opm::Valgrind::CheckDefined(Sw);
452 Opm::Valgrind::CheckDefined(Sg);
454 fs.setSaturation(waterPhaseIdx, Sw);
455 fs.setSaturation(gasPhaseIdx, Sg);
456 fs.setSaturation(naplPhaseIdx, 0);
459 fs.setTemperature(temperature_);
462 Scalar pcAll[numPhases];
464 if (onLeftBoundary_(pos))
466 MaterialLaw::capillaryPressures(pcAll, matParams, fs);
467 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
468 fs.setPressure(phaseIdx, pg + (pcAll[phaseIdx] - pcAll[gasPhaseIdx]));
471 fs.setMoleFraction(gasPhaseIdx, H2OIdx, 1e-6);
472 fs.setMoleFraction(gasPhaseIdx, airIdx,
473 1 - fs.moleFraction(gasPhaseIdx, H2OIdx));
474 fs.setMoleFraction(gasPhaseIdx, NAPLIdx, 0);
476 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
477 typename FluidSystem::template ParameterCache<Scalar> paramCache;
478 CFRP::solve(fs, paramCache, gasPhaseIdx,
482 fs.setMoleFraction(waterPhaseIdx, H2OIdx,
483 1 - fs.moleFraction(waterPhaseIdx, H2OIdx));
486 bool isFineMaterial_(
const GlobalPosition& pos)
const
487 {
return 70. <= pos[0] && pos[0] <= 85. && 7.0 <= pos[1] && pos[1] <= 7.50; }
494 MaterialLawParams materialParams_;
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:306
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:327
void endTimeStep()
Called by the simulator after each time integration.
Definition: infiltrationproblem.hh:262
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:194
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:53
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:282
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: infiltrationproblem.hh:202
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: infiltrationproblem.hh:246
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: infiltrationproblem.hh:346
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:316
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: infiltrationproblem.hh:403
Declares the properties required for the compositional multi-phase primary variable switching model...
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: infiltrationproblem.hh:382
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:292
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
std::string name() const
The problem name.
Definition: infiltrationproblem.hh:252
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394