28 #ifndef EWOMS_FINGER_PROBLEM_HH
29 #define EWOMS_FINGER_PROBLEM_HH
33 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
34 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
35 #include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
36 #include <opm/material/fluidmatrixinteractions/ParkerLenhard.hpp>
37 #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/Air.hpp>
45 #include <ewoms/disc/common/restrictprolong.hh>
48 #include <dune/alugrid/grid.hh>
51 #include <dune/common/version.hh>
52 #include <dune/common/fvector.hh>
53 #include <dune/common/fmatrix.hh>
54 #include <dune/grid/utility/persistentcontainer.hh>
60 template <
class TypeTag>
63 namespace Properties {
73 Dune::nonconforming>);
83 SET_PROP(FingerBaseProblem, WettingPhase)
89 typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
93 SET_PROP(FingerBaseProblem, NonwettingPhase)
99 typedef Opm::GasPhase<Scalar, Opm::Air<Scalar> > type;
103 SET_PROP(FingerBaseProblem, MaterialLaw)
106 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
107 typedef Opm::TwoPhaseMaterialTraits<Scalar,
108 FluidSystem::wettingPhaseIdx,
109 FluidSystem::nonWettingPhaseIdx> Traits;
112 typedef Opm::ParkerLenhard<Traits> ParkerLenhard;
113 typedef ParkerLenhard type;
117 SET_BOOL_PROP(FingerBaseProblem, NewtonWriteConvergence, false);
120 SET_INT_PROP(FingerBaseProblem, NumericDifferenceMethod, +1);
123 SET_INT_PROP(FingerBaseProblem, EnableConstraints, true);
161 template <class TypeTag>
165 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
170 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
171 typedef typename
GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
172 typedef typename
GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
173 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
175 typedef typename
GET_PROP_TYPE(TypeTag, Constraints) Constraints;
182 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
183 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
186 contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
189 dim = GridView::dimension,
190 dimWorld = GridView::dimensionworld
193 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
195 enum { codim = Stencil::Entity::codimension };
197 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
198 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
200 typedef typename
GET_PROP(TypeTag, MaterialLaw)::ParkerLenhard ParkerLenhard;
201 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
202 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
204 typedef typename GridView::ctype CoordScalar;
205 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
206 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
208 typedef typename GridView :: Grid Grid;
210 typedef Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > > MaterialLawParamsContainer;
214 typedef CopyRestrictProlong< Grid, MaterialLawParamsContainer > RestrictProlongOperator;
220 : ParentType(simulator),
221 materialParams_( simulator.gridManager().grid(), codim )
243 std::string(
"finger") +
244 "_" + Model::name() +
245 "_" + Model::discretizationName() +
246 (this->model().enableGridAdaptation()?
"_adaptive":
"");
254 ParentType::registerParameters();
257 "The initial saturation in the domain [] of the "
266 ParentType::finishInit();
270 temperature_ = 273.15 + 20;
276 micParams_.setVgAlpha(0.0037);
277 micParams_.setVgN(4.7);
278 micParams_.finalize();
280 mdcParams_.setVgAlpha(0.0037);
281 mdcParams_.setVgN(4.7);
282 mdcParams_.finalize();
286 materialParams_.resize();
288 for (
auto it = materialParams_.begin(),
289 end = materialParams_.end(); it != end; ++it ) {
290 std::shared_ptr< MaterialLawParams >& materialParams = *it ;
291 if( ! materialParams )
293 materialParams.reset(
new MaterialLawParams() );
294 materialParams->setMicParams(&micParams_);
295 materialParams->setMdcParams(&mdcParams_);
296 materialParams->setSwr(0.0);
297 materialParams->setSnr(0.1);
298 materialParams->finalize();
299 ParkerLenhard::reset(*materialParams);
303 K_ = this->toDimMatrix_(4.6e-10);
305 setupInitialFluidState_();
320 this->model().globalStorage(storage);
323 if (this->gridView().comm().rank() == 0) {
324 std::cout <<
"Storage: " << storage << std::endl << std::flush;
329 ElementContext elemCtx(this->simulator());
331 auto elemIt = this->gridView().template begin<0>();
332 const auto& elemEndIt = this->gridView().template end<0>();
333 for (; elemIt != elemEndIt; ++elemIt) {
334 const auto& elem = *elemIt;
335 elemCtx.updateAll( elem );
336 size_t numDofs = elemCtx.numDof(0);
337 for (
unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
339 MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, 0 );
340 const auto& fs = elemCtx.intensiveQuantities(scvIdx, 0).fluidState();
341 ParkerLenhard::update(materialParam, fs);
356 template <
class Context>
358 {
return temperature_; }
363 template <
class Context>
370 template <
class Context>
371 Scalar
porosity(
const Context& ,
unsigned ,
unsigned )
const
377 template <
class Context>
379 unsigned spaceIdx,
unsigned timeIdx)
381 const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
382 assert(materialParams_[entity]);
383 return *materialParams_[entity];
389 template <
class Context>
391 unsigned spaceIdx,
unsigned timeIdx)
const
393 const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
394 assert(materialParams_[entity]);
395 return *materialParams_[entity];
408 template <
class Context>
409 void boundary(BoundaryRateVector& values,
const Context& context,
410 unsigned spaceIdx,
unsigned timeIdx)
const
412 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
414 if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
417 assert(onUpperBoundary_(pos));
419 values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
425 values[contiWettingEqIdx] = -0.001;
439 template <
class Context>
440 void initial(PrimaryVariables& values,
const Context& ,
unsigned ,
unsigned )
const
443 values.assignNaive(initialFluidState_);
449 template <
class Context>
450 void constraints(Constraints& constraints,
const Context& context,
451 unsigned spaceIdx,
unsigned timeIdx)
const
453 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
455 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
456 constraints.setActive(
true);
457 constraints.assignNaive(initialFluidState_);
459 else if (onLowerBoundary_(pos)) {
460 constraints.setActive(
true);
461 constraints.assignNaive(initialFluidState_);
471 template <
class Context>
472 void source(RateVector& rate,
const Context& ,
473 unsigned ,
unsigned )
const
474 { rate = Scalar(0.0); }
478 bool onLeftBoundary_(
const GlobalPosition& pos)
const
479 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
481 bool onRightBoundary_(
const GlobalPosition& pos)
const
482 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
484 bool onLowerBoundary_(
const GlobalPosition& pos)
const
485 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
487 bool onUpperBoundary_(
const GlobalPosition& pos)
const
488 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
490 bool onInlet_(
const GlobalPosition& pos)
const
492 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
493 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
495 if (!onUpperBoundary_(pos))
498 Scalar xInject[] = { 0.25, 0.75 };
499 Scalar injectLen[] = { 0.1, 0.1 };
500 for (
unsigned i = 0; i <
sizeof(xInject) /
sizeof(Scalar); ++i) {
501 if (xInject[i] - injectLen[i] / 2 < lambda
502 && lambda < xInject[i] + injectLen[i] / 2)
508 void setupInitialFluidState_()
510 auto& fs = initialFluidState_;
511 fs.setPressure(wettingPhaseIdx, 1e5);
514 fs.setSaturation(wettingPhaseIdx, Sw);
515 fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
517 fs.setTemperature(temperature_);
521 fs.setPressure(nonWettingPhaseIdx, pn);
522 fs.setPressure(wettingPhaseIdx, pn);
527 typename MaterialLawParams::VanGenuchtenParams micParams_;
528 typename MaterialLawParams::VanGenuchtenParams mdcParams_;
530 MaterialLawParamsContainer materialParams_;
532 Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
static void registerParameters()
Definition: fingerproblem.hh:252
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
Two-phase problem featuring some gravity-driven saturation fingers.
Definition: fingerproblem.hh:61
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:364
std::string name() const
The problem name.
Definition: fingerproblem.hh:241
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: fingerproblem.hh:472
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: fingerproblem.hh:450
Helper class for grid instantiation of the lens problem.
Definition: structuredgridmanager.hh:51
Definition: restrictprolong.hh:35
#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
#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
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
RestrictProlongOperator restrictProlongOperator()
return restriction and prolongation operator
Definition: fingerproblem.hh:233
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition: fingerproblem.hh:440
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: fingerproblem.hh:264
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:357
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:371
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx)
Definition: fingerproblem.hh:378
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: fingerproblem.hh:409
void endTimeStep()
Called by the simulator after each time integration.
Definition: fingerproblem.hh:311
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 GET_PROP(TypeTag, PropTagName)
Retrieve a property for a type tag.
Definition: propertysystem.hh:454
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:390
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394