All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
fingerproblem.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_FINGER_PROBLEM_HH
29 #define EWOMS_FINGER_PROBLEM_HH
30 
32 
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>
38 
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>
43 
45 #include <ewoms/disc/common/restrictprolong.hh>
46 
47 #if HAVE_DUNE_ALUGRID
48 #include <dune/alugrid/grid.hh>
49 #endif
50 
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>
55 
56 #include <vector>
57 #include <string>
58 
59 namespace Ewoms {
60 template <class TypeTag>
62 
63 namespace Properties {
65 
66 #if HAVE_DUNE_ALUGRID
67 // use dune-alugrid if available
68 SET_TYPE_PROP(FingerBaseProblem,
69  Grid,
70  Dune::ALUGrid</*dim=*/2,
71  /*dimWorld=*/2,
72  Dune::cube,
73  Dune::nonconforming>);
74 #endif
75 
76 // declare the properties used by the finger problem
77 NEW_PROP_TAG(InitialWaterSaturation);
78 
79 // Set the problem property
80 SET_TYPE_PROP(FingerBaseProblem, Problem, Ewoms::FingerProblem<TypeTag>);
81 
82 // Set the wetting phase
83 SET_PROP(FingerBaseProblem, WettingPhase)
84 {
85 private:
86  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
87 
88 public:
89  typedef Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> > type;
90 };
91 
92 // Set the non-wetting phase
93 SET_PROP(FingerBaseProblem, NonwettingPhase)
94 {
95 private:
96  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
97 
98 public:
99  typedef Opm::GasPhase<Scalar, Opm::Air<Scalar> > type;
100 };
101 
102 // Set the material Law
103 SET_PROP(FingerBaseProblem, MaterialLaw)
104 {
105  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
106  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
107  typedef Opm::TwoPhaseMaterialTraits<Scalar,
108  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
109  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx> Traits;
110 
111  // use the parker-lenhard hysteresis law
112  typedef Opm::ParkerLenhard<Traits> ParkerLenhard;
113  typedef ParkerLenhard type;
114 };
115 
116 // Write the solutions of individual newton iterations?
117 SET_BOOL_PROP(FingerBaseProblem, NewtonWriteConvergence, false);
118 
119 // Use forward differences instead of central differences
120 SET_INT_PROP(FingerBaseProblem, NumericDifferenceMethod, +1);
121 
122 // Enable constraints
123 SET_INT_PROP(FingerBaseProblem, EnableConstraints, true);
124 
125 // Enable gravity
126 SET_BOOL_PROP(FingerBaseProblem, EnableGravity, true);
127 
128 // define the properties specific for the finger problem
129 SET_SCALAR_PROP(FingerBaseProblem, DomainSizeX, 0.1);
130 SET_SCALAR_PROP(FingerBaseProblem, DomainSizeY, 0.3);
131 SET_SCALAR_PROP(FingerBaseProblem, DomainSizeZ, 0.1);
132 
133 SET_SCALAR_PROP(FingerBaseProblem, InitialWaterSaturation, 0.01);
134 
135 SET_INT_PROP(FingerBaseProblem, CellsX, 20);
136 SET_INT_PROP(FingerBaseProblem, CellsY, 70);
137 SET_INT_PROP(FingerBaseProblem, CellsZ, 1);
138 
139 // The default for the end time of the simulation
140 SET_SCALAR_PROP(FingerBaseProblem, EndTime, 215);
141 
142 // The default for the initial time step size of the simulation
143 SET_SCALAR_PROP(FingerBaseProblem, InitialTimeStepSize, 10);
144 } // namespace Properties
145 
161 template <class TypeTag>
162 class FingerProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
163 {
165  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
166 
167  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
168  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
169  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
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;
176  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
177 
178  enum {
179  // number of phases
180 
181  // phase indices
182  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
183  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
184 
185  // equation indices
186  contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
187 
188  // Grid and world dimension
189  dim = GridView::dimension,
190  dimWorld = GridView::dimensionworld
191  };
192 
193  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
194  typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
195  enum { codim = Stencil::Entity::codimension };
196  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
197  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
198  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
199 
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;
203 
204  typedef typename GridView::ctype CoordScalar;
205  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
206  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
207 
208  typedef typename GridView :: Grid Grid;
209 
210  typedef Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > > MaterialLawParamsContainer;
212 
213 public:
214  typedef CopyRestrictProlong< Grid, MaterialLawParamsContainer > RestrictProlongOperator;
215 
220  : ParentType(simulator),
221  materialParams_( simulator.gridManager().grid(), codim )
222  {
223  }
224 
228 
234  {
235  return RestrictProlongOperator( materialParams_ );
236  }
237 
241  std::string name() const
242  { return
243  std::string("finger") +
244  "_" + Model::name() +
245  "_" + Model::discretizationName() +
246  (this->model().enableGridAdaptation()?"_adaptive":"");
247  }
248 
252  static void registerParameters()
253  {
254  ParentType::registerParameters();
255 
256  EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation,
257  "The initial saturation in the domain [] of the "
258  "wetting phase");
259  }
260 
264  void finishInit()
265  {
266  ParentType::finishInit();
267 
268  eps_ = 3e-6;
269 
270  temperature_ = 273.15 + 20; // -> 20°C
271 
272  FluidSystem::init();
273 
274  // parameters for the Van Genuchten law of the main imbibition
275  // and the main drainage curves.
276  micParams_.setVgAlpha(0.0037);
277  micParams_.setVgN(4.7);
278  micParams_.finalize();
279 
280  mdcParams_.setVgAlpha(0.0037);
281  mdcParams_.setVgN(4.7);
282  mdcParams_.finalize();
283 
284  // initialize the material parameter objects of the individual
285  // finite volumes, resize will resize the container to the number of elements
286  materialParams_.resize();
287 
288  for (auto it = materialParams_.begin(),
289  end = materialParams_.end(); it != end; ++it ) {
290  std::shared_ptr< MaterialLawParams >& materialParams = *it ;
291  if( ! materialParams )
292  {
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);
300  }
301  }
302 
303  K_ = this->toDimMatrix_(4.6e-10);
304 
305  setupInitialFluidState_();
306  }
307 
311  void endTimeStep()
312  {
313 #ifndef NDEBUG
314  // checkConservativeness() does not include the effect of constraints, so we
315  // disable it for this problem...
316  //this->model().checkConservativeness();
317 
318  // Calculate storage terms
319  EqVector storage;
320  this->model().globalStorage(storage);
321 
322  // Write mass balance information for rank 0
323  if (this->gridView().comm().rank() == 0) {
324  std::cout << "Storage: " << storage << std::endl << std::flush;
325  }
326 #endif // NDEBUG
327 
328  // update the history of the hysteresis law
329  ElementContext elemCtx(this->simulator());
330 
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(/*timeIdx=*/0);
337  for (unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
338  {
339  MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, /*timeIdx=*/0 );
340  const auto& fs = elemCtx.intensiveQuantities(scvIdx, /*timeIdx=*/0).fluidState();
341  ParkerLenhard::update(materialParam, fs);
342  }
343  }
344  }
345 
347 
351 
356  template <class Context>
357  Scalar temperature(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
358  { return temperature_; }
359 
363  template <class Context>
364  const DimMatrix& intrinsicPermeability(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
365  { return K_; }
366 
370  template <class Context>
371  Scalar porosity(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
372  { return 0.4; }
373 
377  template <class Context>
378  MaterialLawParams& materialLawParams(const Context& context,
379  unsigned spaceIdx, unsigned timeIdx)
380  {
381  const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
382  assert(materialParams_[entity]);
383  return *materialParams_[entity];
384  }
385 
389  template <class Context>
390  const MaterialLawParams& materialLawParams(const Context& context,
391  unsigned spaceIdx, unsigned timeIdx) const
392  {
393  const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
394  assert(materialParams_[entity]);
395  return *materialParams_[entity];
396  }
397 
399 
403 
408  template <class Context>
409  void boundary(BoundaryRateVector& values, const Context& context,
410  unsigned spaceIdx, unsigned timeIdx) const
411  {
412  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
413 
414  if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
415  values.setNoFlow();
416  else {
417  assert(onUpperBoundary_(pos));
418 
419  values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
420  }
421 
422  // override the value for the liquid phase by forced
423  // imbibition of water on inlet boundary segments
424  if (onInlet_(pos)) {
425  values[contiWettingEqIdx] = -0.001; // [kg/(m^2 s)]
426  }
427  }
428 
430 
434 
439  template <class Context>
440  void initial(PrimaryVariables& values, const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
441  {
442  // assign the primary variables
443  values.assignNaive(initialFluidState_);
444  }
445 
449  template <class Context>
450  void constraints(Constraints& constraints, const Context& context,
451  unsigned spaceIdx, unsigned timeIdx) const
452  {
453  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
454 
455  if (onUpperBoundary_(pos) && !onInlet_(pos)) {
456  constraints.setActive(true);
457  constraints.assignNaive(initialFluidState_);
458  }
459  else if (onLowerBoundary_(pos)) {
460  constraints.setActive(true);
461  constraints.assignNaive(initialFluidState_);
462  }
463  }
464 
471  template <class Context>
472  void source(RateVector& rate, const Context& /*context*/,
473  unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
474  { rate = Scalar(0.0); }
476 
477 private:
478  bool onLeftBoundary_(const GlobalPosition& pos) const
479  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
480 
481  bool onRightBoundary_(const GlobalPosition& pos) const
482  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
483 
484  bool onLowerBoundary_(const GlobalPosition& pos) const
485  { return pos[1] < this->boundingBoxMin()[1] + eps_; }
486 
487  bool onUpperBoundary_(const GlobalPosition& pos) const
488  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
489 
490  bool onInlet_(const GlobalPosition& pos) const
491  {
492  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
493  Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
494 
495  if (!onUpperBoundary_(pos))
496  return false;
497 
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)
503  return true;
504  }
505  return false;
506  }
507 
508  void setupInitialFluidState_()
509  {
510  auto& fs = initialFluidState_;
511  fs.setPressure(wettingPhaseIdx, /*pressure=*/1e5);
512 
513  Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
514  fs.setSaturation(wettingPhaseIdx, Sw);
515  fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
516 
517  fs.setTemperature(temperature_);
518 
519  // set the absolute pressures
520  Scalar pn = 1e5;
521  fs.setPressure(nonWettingPhaseIdx, pn);
522  fs.setPressure(wettingPhaseIdx, pn);
523  }
524 
525  DimMatrix K_;
526 
527  typename MaterialLawParams::VanGenuchtenParams micParams_;
528  typename MaterialLawParams::VanGenuchtenParams mdcParams_;
529 
530  MaterialLawParamsContainer materialParams_;
531 
532  Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
533 
534  Scalar temperature_;
535  Scalar eps_;
536 };
537 
538 } // namespace Ewoms
539 
540 #endif
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