28 #ifndef EWOMS_LENS_PROBLEM_HH 29 #define EWOMS_LENS_PROBLEM_HH 35 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp> 36 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp> 37 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp> 38 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp> 39 #include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp> 40 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp> 41 #include <opm/material/components/SimpleH2O.hpp> 42 #include <opm/material/components/Dnapl.hpp> 43 #include <opm/common/Unused.hpp> 45 #include <dune/common/version.hh> 46 #include <dune/common/fvector.hh> 47 #include <dune/common/fmatrix.hh> 54 template <
class TypeTag>
57 namespace Properties {
75 SET_PROP(LensBaseProblem, WettingPhase)
81 typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
85 SET_PROP(LensBaseProblem, NonwettingPhase)
91 typedef Opm::LiquidPhase<Scalar, Opm::DNAPL<Scalar> > type;
95 SET_PROP(LensBaseProblem, MaterialLaw)
98 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
99 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
100 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
103 typedef Opm::TwoPhaseMaterialTraits<Scalar,
104 FluidSystem::wettingPhaseIdx,
105 FluidSystem::nonWettingPhaseIdx> Traits;
109 typedef Opm::RegularizedVanGenuchten<Traits> EffectiveLaw;
113 typedef Opm::EffToAbsLaw<EffectiveLaw> type;
117 SET_BOOL_PROP(LensBaseProblem, NewtonWriteConvergence,
false);
120 SET_INT_PROP(LensBaseProblem, NumericDifferenceMethod, +1);
148 SET_BOOL_PROP(LensBaseProblem, VtkWriteIntrinsicPermeabilities,
true);
154 SET_BOOL_PROP(LensBaseProblem, EnableIntensiveQuantityCache,
true);
180 template <
class TypeTag>
181 class LensProblem :
public GET_PROP_TYPE(TypeTag, BaseProblem)
183 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
185 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
186 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
187 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
188 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
189 typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
190 typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
191 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
192 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
193 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
197 numPhases = FluidSystem::numPhases,
200 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
201 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
204 contiNEqIdx = Indices::conti0EqIdx + nonWettingPhaseIdx,
207 dim = GridView::dimension,
208 dimWorld = GridView::dimensionworld
211 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
212 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
213 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
214 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
215 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
217 typedef typename GridView::ctype CoordScalar;
218 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
220 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
227 : ParentType(simulator)
235 ParentType::finishInit();
240 temperature_ = 273.15 + 20;
243 lensUpperRight_[0] =
EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightX);
244 lensUpperRight_[1] =
EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightY);
248 lensUpperRight_[2] =
EWOMS_GET_PARAM(TypeTag, Scalar, LensUpperRightZ);
252 lensMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.18);
253 lensMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
254 outerMaterialParams_.setResidualSaturation(wettingPhaseIdx, 0.05);
255 outerMaterialParams_.setResidualSaturation(nonWettingPhaseIdx, 0.0);
258 lensMaterialParams_.setVgAlpha(0.00045);
259 lensMaterialParams_.setVgN(7.3);
260 outerMaterialParams_.setVgAlpha(0.0037);
261 outerMaterialParams_.setVgN(4.7);
263 lensMaterialParams_.finalize();
264 outerMaterialParams_.finalize();
266 lensK_ = this->toDimMatrix_(9.05e-12);
267 outerK_ = this->toDimMatrix_(4.6e-10);
271 this->gravity_[1] = -9.81;
280 ParentType::registerParameters();
283 "The x-coordinate of the lens' lower-left corner " 286 "The y-coordinate of the lens' lower-left corner " 289 "The x-coordinate of the lens' upper-right corner " 292 "The y-coordinate of the lens' upper-right corner " 297 "The z-coordinate of the lens' lower-left " 300 "The z-coordinate of the lens' upper-right " 313 template <
class Context>
315 unsigned timeIdx)
const 317 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
319 if (isInLens_(globalPos))
327 template <
class Context>
329 unsigned spaceIdx OPM_UNUSED,
330 unsigned timeIdx OPM_UNUSED)
const 336 template <
class Context>
338 unsigned spaceIdx,
unsigned timeIdx)
const 340 const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
342 if (isInLens_(globalPos))
343 return lensMaterialParams_;
344 return outerMaterialParams_;
350 template <
class Context>
352 unsigned spaceIdx OPM_UNUSED,
353 unsigned timeIdx OPM_UNUSED)
const 354 {
return temperature_; }
368 typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizerSplice) LLS;
370 bool useAutoDiff = std::is_same<LLS, TTAG(AutoDiffLocalLinearizer)>::value;
372 std::ostringstream oss;
373 oss <<
"lens_" << Model::name()
374 <<
"_" << Model::discretizationName()
375 <<
"_" << (useAutoDiff?
"ad":
"fd");
397 this->model().checkConservativeness();
401 this->model().globalStorage(storage);
404 if (this->gridView().comm().rank() == 0) {
405 std::cout <<
"Storage: " << storage << std::endl << std::flush;
420 template <
class Context>
422 const Context& context,
424 unsigned timeIdx)
const 426 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
428 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
430 Scalar densityW = WettingPhase::density(temperature_,
433 Scalar T =
temperature(context, spaceIdx, timeIdx);
437 if (onLeftBoundary_(pos)) {
438 Scalar height = this->boundingBoxMax()[1] - this->boundingBoxMin()[1];
439 Scalar depth = this->boundingBoxMax()[1] - pos[1];
440 Scalar alpha = (1 + 1.5 / height);
443 pw = 1e5 - alpha * densityW * this->gravity()[1] * depth;
447 Scalar depth = this->boundingBoxMax()[1] - pos[1];
450 pw = 1e5 - densityW * this->gravity()[1] * depth;
455 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
457 Opm::ImmiscibleFluidState<Scalar, FluidSystem,
459 fs.setSaturation(wettingPhaseIdx, Sw);
460 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
461 fs.setTemperature(T);
463 Scalar pC[numPhases];
464 MaterialLaw::capillaryPressures(pC, matParams, fs);
465 fs.setPressure(wettingPhaseIdx, pw);
466 fs.setPressure(nonWettingPhaseIdx, pw + pC[nonWettingPhaseIdx] - pC[wettingPhaseIdx]);
469 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
471 else if (onInlet_(pos)) {
472 RateVector massRate(0.0);
474 massRate[contiNEqIdx] = -0.04;
477 values.setMassRate(massRate);
495 template <
class Context>
496 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const 498 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
499 Scalar depth = this->boundingBoxMax()[1] - pos[1];
501 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
502 fs.setPressure(wettingPhaseIdx, 1e5);
505 fs.setSaturation(wettingPhaseIdx, Sw);
506 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
508 fs.setTemperature(temperature_);
510 typename FluidSystem::template ParameterCache<Scalar> paramCache;
511 paramCache.updatePhase(fs, wettingPhaseIdx);
512 Scalar densityW = FluidSystem::density(fs, paramCache, wettingPhaseIdx);
515 Scalar pw = 1e5 - densityW * this->gravity()[1] * depth;
518 const MaterialLawParams& matParams = this->
materialLawParams(context, spaceIdx, timeIdx);
519 Scalar pC[numPhases];
520 MaterialLaw::capillaryPressures(pC, matParams, fs);
523 fs.setPressure(wettingPhaseIdx, pw);
524 fs.setPressure(nonWettingPhaseIdx, pw + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]));
527 values.assignNaive(fs);
536 template <
class Context>
538 const Context& context OPM_UNUSED,
539 unsigned spaceIdx OPM_UNUSED,
540 unsigned timeIdx OPM_UNUSED)
const 541 { rate = Scalar(0.0); }
546 bool isInLens_(
const GlobalPosition& pos)
const 548 for (
unsigned i = 0; i < dim; ++i) {
549 if (pos[i] < lensLowerLeft_[i] - eps_ || pos[i] > lensUpperRight_[i]
556 bool onLeftBoundary_(
const GlobalPosition& pos)
const 557 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
559 bool onRightBoundary_(
const GlobalPosition& pos)
const 560 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
562 bool onLowerBoundary_(
const GlobalPosition& pos)
const 563 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
565 bool onUpperBoundary_(
const GlobalPosition& pos)
const 566 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
568 bool onInlet_(
const GlobalPosition& pos)
const 570 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
571 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
572 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
575 GlobalPosition lensLowerLeft_;
576 GlobalPosition lensUpperRight_;
580 MaterialLawParams lensMaterialParams_;
581 MaterialLawParams outerMaterialParams_;
static void registerParameters()
Definition: lensproblem.hh:278
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:314
Definition: baseauxiliarymodule.hh:37
Helper class for grid instantiation of the lens problem.
Definition: structuredgridmanager.hh:51
LensProblem(Simulator &simulator)
Definition: lensproblem.hh:226
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: lensproblem.hh:537
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
Helper class for grid instantiation of the lens problem.
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: lensproblem.hh:421
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
void beginTimeStep()
Called by the simulator before each time integration.
Definition: lensproblem.hh:382
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
void beginIteration()
Called by the simulator before each Newton-Raphson iteration.
Definition: lensproblem.hh:388
void endTimeStep()
Called by the simulator after each time integration.
Definition: lensproblem.hh:394
std::string name() const
The problem name.
Definition: lensproblem.hh:366
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: lensproblem.hh:351
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
Calculates the local residual and its Jacobian for a single element of the grid.
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: lensproblem.hh:233
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: lensproblem.hh:328
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: lensproblem.hh:337
Defines the properties required for the immiscible multi-phase model.
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: lensproblem.hh:496
Soil contamination problem where DNAPL infiltrates a fully water saturated medium.
Definition: lensproblem.hh:55
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394