28 #ifndef EWOMS_RICHARDS_LENS_PROBLEM_HH
29 #define EWOMS_RICHARDS_LENS_PROBLEM_HH
33 #include <opm/material/components/SimpleH2O.hpp>
34 #include <opm/material/fluidsystems/LiquidPhase.hpp>
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/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>
49 template <
class TypeTag>
52 namespace Properties {
68 typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
75 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76 enum { wettingPhaseIdx = FluidSystem::wettingPhaseIdx };
77 enum { nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx };
80 typedef Opm::TwoPhaseMaterialTraits<Scalar,
81 FluidSystem::wettingPhaseIdx,
82 FluidSystem::nonWettingPhaseIdx>
87 typedef Opm::RegularizedVanGenuchten<Traits> EffectiveLaw;
91 typedef Opm::EffToAbsLaw<EffectiveLaw> type;
98 SET_INT_PROP(RichardsLensProblem, NumericDifferenceMethod, 0);
101 SET_INT_PROP(RichardsLensProblem, NewtonMaxIterations, 28);
104 SET_INT_PROP(RichardsLensProblem, NewtonTargetIterations, 18);
107 SET_BOOL_PROP(RichardsLensProblem, NewtonWriteConvergence, false);
116 SET_STRING_PROP(RichardsLensProblem, GridFile, "./data/richardslens_24x16.dgf");
135 template <class TypeTag>
136 class RichardsLensProblem : public
GET_PROP_TYPE(TypeTag, BaseProblem)
138 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
142 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
143 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
144 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
146 typedef typename
GET_PROP_TYPE(TypeTag, Simulator) Simulator;
148 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
154 pressureWIdx = Indices::pressureWIdx,
155 contiEqIdx = Indices::contiEqIdx,
156 wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
157 nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
158 numPhases = FluidSystem::numPhases,
161 dimWorld = GridView::dimensionworld
165 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
167 typedef typename MaterialLaw::Params MaterialLawParams;
169 typedef typename GridView::ctype CoordScalar;
170 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
171 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
172 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
179 : ParentType(simulator)
182 dofIsInLens_.resize(simulator.model().numGridDof());
190 ParentType::finishInit();
195 lensLowerLeft_[0] = 1.0;
196 lensLowerLeft_[1] = 2.0;
198 lensUpperRight_[0] = 4.0;
199 lensUpperRight_[1] = 3.0;
203 lensMaterialParams_.setVgAlpha(0.00045);
204 lensMaterialParams_.setVgN(7.3);
205 lensMaterialParams_.finalize();
207 outerMaterialParams_.setVgAlpha(0.0037);
208 outerMaterialParams_.setVgN(4.7);
209 outerMaterialParams_.finalize();
218 lensK_ = this->toDimMatrix_(1e-12);
219 outerK_ = this->toDimMatrix_(5e-12);
222 Stencil stencil(this->gridView(), this->simulator().model().dofMapper() );
223 auto elemIt = this->gridView().template begin<0>();
224 auto elemEndIt = this->gridView().template end<0>();
225 for (; elemIt != elemEndIt; ++elemIt) {
226 stencil.update(*elemIt);
227 for (
unsigned dofIdx = 0; dofIdx < stencil.numPrimaryDof(); ++ dofIdx) {
228 unsigned globalDofIdx = stencil.globalSpaceIndex(dofIdx);
229 const auto& dofPos = stencil.subControlVolume(dofIdx).center();
230 dofIsInLens_[globalDofIdx] = isInLens_(dofPos);
245 std::ostringstream oss;
246 oss <<
"lens_richards_"
247 << Model::discretizationName();
257 this->model().checkConservativeness();
261 this->model().globalStorage(storage);
264 if (this->gridView().comm().rank() == 0) {
265 std::cout <<
"Storage: " << storage << std::endl << std::flush;
273 template <
class Context>
274 Scalar
temperature(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
275 {
return temperature(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
277 Scalar temperature(
unsigned globalSpaceIdx OPM_UNUSED,
unsigned timeIdx OPM_UNUSED)
const
278 {
return 273.15 + 10; }
283 template <
class Context>
286 unsigned timeIdx)
const
288 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
297 template <
class Context>
299 unsigned spaceIdx OPM_UNUSED,
300 unsigned timeIdx OPM_UNUSED)
const
306 template <
class Context>
309 unsigned timeIdx)
const
311 unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
312 return materialLawParams(globalSpaceIdx, timeIdx);
315 const MaterialLawParams& materialLawParams(
unsigned globalSpaceIdx,
316 unsigned timeIdx OPM_UNUSED)
const
318 if (dofIsInLens_[globalSpaceIdx])
319 return lensMaterialParams_;
320 return outerMaterialParams_;
328 template <
class Context>
331 unsigned timeIdx)
const
332 {
return referencePressure(context.globalSpaceIndex(spaceIdx, timeIdx), timeIdx); }
336 Scalar referencePressure(
unsigned globalSpaceIdx OPM_UNUSED,
337 unsigned timeIdx OPM_UNUSED)
const
350 template <
class Context>
352 const Context& context,
354 unsigned timeIdx)
const
356 const auto& pos = context.pos(spaceIdx, timeIdx);
358 if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
359 const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
362 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
363 fs.setSaturation(wettingPhaseIdx, Sw);
364 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
367 MaterialLaw::capillaryPressures(pC, materialParams, fs);
368 fs.setPressure(wettingPhaseIdx, pnRef_ + pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
369 fs.setPressure(nonWettingPhaseIdx, pnRef_);
371 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
373 else if (onInlet_(pos)) {
374 RateVector massRate(0.0);
377 massRate[contiEqIdx] = -0.04;
379 values.setMassRate(massRate);
395 template <
class Context>
397 const Context& context,
399 unsigned timeIdx)
const
401 const auto& materialParams = this->materialLawParams(context, spaceIdx, timeIdx);
404 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
405 fs.setSaturation(wettingPhaseIdx, Sw);
406 fs.setSaturation(nonWettingPhaseIdx, 1.0 - Sw);
409 MaterialLaw::capillaryPressures(pC, materialParams, fs);
410 values[pressureWIdx] = pnRef_ + (pC[wettingPhaseIdx] - pC[nonWettingPhaseIdx]);
419 template <
class Context>
421 const Context& context OPM_UNUSED,
422 unsigned spaceIdx OPM_UNUSED,
423 unsigned timeIdx OPM_UNUSED)
const
424 { rate = Scalar(0.0); }
429 bool onLeftBoundary_(
const GlobalPosition& pos)
const
430 {
return pos[0] < this->boundingBoxMin()[0] + eps_; }
432 bool onRightBoundary_(
const GlobalPosition& pos)
const
433 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
435 bool onLowerBoundary_(
const GlobalPosition& pos)
const
436 {
return pos[1] < this->boundingBoxMin()[1] + eps_; }
438 bool onUpperBoundary_(
const GlobalPosition& pos)
const
439 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
441 bool onInlet_(
const GlobalPosition& pos)
const
443 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
444 Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
445 return onUpperBoundary_(pos) && 0.5 < lambda && lambda < 2.0 / 3.0;
448 bool isInLens_(
const GlobalPosition& pos)
const
450 for (
unsigned i = 0; i < dimWorld; ++i) {
451 if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
457 GlobalPosition lensLowerLeft_;
458 GlobalPosition lensUpperRight_;
462 MaterialLawParams lensMaterialParams_;
463 MaterialLawParams outerMaterialParams_;
465 std::vector<bool> dofIsInLens_;
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
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: richardslensproblem.hh:420
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:298
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: richardslensproblem.hh:188
#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
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:274
#define INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
void endTimeStep()
Called by the simulator after each time integration.
Definition: richardslensproblem.hh:254
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:329
std::string name() const
The problem name.
Definition: richardslensproblem.hh:243
This model implements a variant of the Richards equation for quasi-twophase flow. ...
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:284
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:307
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain...
Definition: richardslensproblem.hh:50
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: richardslensproblem.hh:396
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#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: richardslensproblem.hh:351
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394