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>
146 class InfiltrationProblem :
public GET_PROP_TYPE(TypeTag, BaseProblem)
148 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
150 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
151 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
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;
155 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
156 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
157 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
158 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
159 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
160 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
163 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
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);
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));
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_;
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
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void endTimeStep()
Called by the simulator after each time integration.
Definition: infiltrationproblem.hh:262
std::string name() const
The problem name.
Definition: infiltrationproblem.hh:252
Definition: baseauxiliarymodule.hh:37
InfiltrationProblem(Simulator &simulator)
Definition: infiltrationproblem.hh:194
Scalar heatCapacitySolid(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:327
bool shouldWriteRestartFile() const
Returns true if a restart file should be written to disk.
Definition: infiltrationproblem.hh:246
Isothermal NAPL infiltration problem where LNAPL contaminates the unsaturated and the saturated groun...
Definition: infiltrationproblem.hh:53
#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
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: infiltrationproblem.hh:382
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:306
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:282
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: infiltrationproblem.hh:292
Declares the properties required for the compositional multi-phase primary variable switching model...
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: infiltrationproblem.hh:316
#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
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_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394