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>
162 class FingerProblem :
public GET_PROP_TYPE(TypeTag, BaseProblem)
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;
174 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
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)
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>
451 unsigned spaceIdx,
unsigned timeIdx)
const 453 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
455 if (onUpperBoundary_(pos) && !onInlet_(pos)) {
459 else if (onLowerBoundary_(pos)) {
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
Definition: baseauxiliarymodule.hh:37
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
Helper class for grid instantiation of the lens problem.
Definition: structuredgridmanager.hh:51
std::string name() const
The problem name.
Definition: fingerproblem.hh:241
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition: fingerproblem.hh:440
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:371
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
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:364
RestrictProlongOperator restrictProlongOperator()
return restriction and prolongation operator
Definition: fingerproblem.hh:233
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:390
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
#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 endTimeStep()
Called by the simulator after each time integration.
Definition: fingerproblem.hh:311
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: fingerproblem.hh:450
Defines the properties required for the immiscible multi-phase model.
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: fingerproblem.hh:409
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
FingerProblem(Simulator &simulator)
Definition: fingerproblem.hh:219
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394