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;
140 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
141 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
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;
145 typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
146 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
147 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
148 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
149 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
151 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
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);
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); }
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_;
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
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:307
std::string name() const
The problem name.
Definition: richardslensproblem.hh:243
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: richardslensproblem.hh:298
Definition: baseauxiliarymodule.hh:37
RichardsLensProblem(Simulator &simulator)
Definition: richardslensproblem.hh:178
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: richardslensproblem.hh:188
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: richardslensproblem.hh:396
#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 INHERITS_FROM(...)
Syntactic sugar for NEW_TYPE_TAG.
Definition: propertysystem.hh:230
Scalar referencePressure(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Return the reference pressure [Pa] of the wetting phase.
Definition: richardslensproblem.hh:329
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
This model implements a variant of the Richards equation for quasi-twophase flow. ...
Scalar temperature(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:274
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: richardslensproblem.hh:351
A water infiltration problem with a low-permeability lens embedded into a high-permeability domain...
Definition: richardslensproblem.hh:50
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: richardslensproblem.hh:284
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
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394