My Project
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 
31 #include <opm/models/io/structuredgridvanguard.hh>
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 
44 #include <opm/models/immiscible/immiscibleproperties.hh>
45 #include <opm/models/discretization/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 Opm {
60 template <class TypeTag>
61 class FingerProblem;
62 
63 } // namespace Opm
64 
65 namespace Opm::Properties {
66 
67 // Create new type tags
68 namespace TTag {
69 struct FingerBaseProblem { using InheritsFrom = std::tuple<StructuredGridVanguard>; };
70 } // end namespace TTag
71 
72 #if HAVE_DUNE_ALUGRID
73 // use dune-alugrid if available
74 template<class TypeTag>
75 struct Grid<TypeTag, TTag::FingerBaseProblem>
76 { using type = Dune::ALUGrid</*dim=*/2,
77  /*dimWorld=*/2,
78  Dune::cube,
79  Dune::nonconforming>; };
80 #endif
81 
82 // declare the properties used by the finger problem
83 template<class TypeTag, class MyTypeTag>
84 struct InitialWaterSaturation { using type = UndefinedProperty; };
85 
86 // Set the problem property
87 template<class TypeTag>
88 struct Problem<TypeTag, TTag::FingerBaseProblem> { using type = Opm::FingerProblem<TypeTag>; };
89 
90 // Set the wetting phase
91 template<class TypeTag>
92 struct WettingPhase<TypeTag, TTag::FingerBaseProblem>
93 {
94 private:
95  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
96 
97 public:
98  using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
99 };
100 
101 // Set the non-wetting phase
102 template<class TypeTag>
103 struct NonwettingPhase<TypeTag, TTag::FingerBaseProblem>
104 {
105 private:
106  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
107 
108 public:
109  using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
110 };
111 
112 // Set the material Law
113 template<class TypeTag>
114 struct MaterialLaw<TypeTag, TTag::FingerBaseProblem>
115 {
116  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
117  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
118  using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
119  /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
120  /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
121 
122  // use the parker-lenhard hysteresis law
123  using ParkerLenhard = Opm::ParkerLenhard<Traits>;
124  using type = ParkerLenhard;
125 };
126 
127 // Write the solutions of individual newton iterations?
128 template<class TypeTag>
129 struct NewtonWriteConvergence<TypeTag, TTag::FingerBaseProblem> { static constexpr bool value = false; };
130 
131 // Use forward differences instead of central differences
132 template<class TypeTag>
133 struct NumericDifferenceMethod<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = +1; };
134 
135 // Enable constraints
136 template<class TypeTag>
137 struct EnableConstraints<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = true; };
138 
139 // Enable gravity
140 template<class TypeTag>
141 struct EnableGravity<TypeTag, TTag::FingerBaseProblem> { static constexpr bool value = true; };
142 
143 // define the properties specific for the finger problem
144 template<class TypeTag>
145 struct DomainSizeX<TypeTag, TTag::FingerBaseProblem>
146 {
147  using type = GetPropType<TypeTag, Scalar>;
148  static constexpr type value = 0.1;
149 };
150 template<class TypeTag>
151 struct DomainSizeY<TypeTag, TTag::FingerBaseProblem>
152 {
153  using type = GetPropType<TypeTag, Scalar>;
154  static constexpr type value = 0.3;
155 };
156 template<class TypeTag>
157 struct DomainSizeZ<TypeTag, TTag::FingerBaseProblem>
158 {
159  using type = GetPropType<TypeTag, Scalar>;
160  static constexpr type value = 0.1;
161 };
162 
163 template<class TypeTag>
164 struct InitialWaterSaturation<TypeTag, TTag::FingerBaseProblem>
165 {
166  using type = GetPropType<TypeTag, Scalar>;
167  static constexpr type value = 0.01;
168 };
169 
170 template<class TypeTag>
171 struct CellsX<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 20; };
172 template<class TypeTag>
173 struct CellsY<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 70; };
174 template<class TypeTag>
175 struct CellsZ<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = 1; };
176 
177 // The default for the end time of the simulation
178 template<class TypeTag>
179 struct EndTime<TypeTag, TTag::FingerBaseProblem>
180 {
181  using type = GetPropType<TypeTag, Scalar>;
182  static constexpr type value = 215;
183 };
184 
185 // The default for the initial time step size of the simulation
186 template<class TypeTag>
187 struct InitialTimeStepSize<TypeTag, TTag::FingerBaseProblem>
188 {
189  using type = GetPropType<TypeTag, Scalar>;
190  static constexpr type value = 10;
191 };
192 
193 } // namespace Opm::Properties
194 
195 namespace Opm {
196 
212 template <class TypeTag>
213 class FingerProblem : public GetPropType<TypeTag, Properties::BaseProblem>
214 {
216  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
217 
218  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
219  using GridView = GetPropType<TypeTag, Properties::GridView>;
220  using Indices = GetPropType<TypeTag, Properties::Indices>;
221  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
222  using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
223  using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
224  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
225  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
226  using Constraints = GetPropType<TypeTag, Properties::Constraints>;
227  using Model = GetPropType<TypeTag, Properties::Model>;
228 
229  enum {
230  // number of phases
231  numPhases = FluidSystem::numPhases,
232 
233  // phase indices
234  wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
235  nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
236 
237  // equation indices
238  contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
239 
240  // Grid and world dimension
241  dim = GridView::dimension,
242  dimWorld = GridView::dimensionworld
243  };
244 
245  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
246  using Stencil = GetPropType<TypeTag, Properties::Stencil> ;
247  enum { codim = Stencil::Entity::codimension };
248  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
249  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
250  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
251 
252  using ParkerLenhard = typename GetProp<TypeTag, Properties::MaterialLaw>::ParkerLenhard;
253  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
254  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
255 
256  using CoordScalar = typename GridView::ctype;
257  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
258  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
259 
260  using Grid = typename GridView :: Grid;
261 
262  using MaterialLawParamsContainer = Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > > ;
264 
265 public:
266  using RestrictProlongOperator = CopyRestrictProlong< Grid, MaterialLawParamsContainer >;
267 
271  FingerProblem(Simulator& simulator)
272  : ParentType(simulator),
273  materialParams_( simulator.vanguard().grid(), codim )
274  {
275  }
276 
281 
285  RestrictProlongOperator restrictProlongOperator()
286  {
287  return RestrictProlongOperator( materialParams_ );
288  }
289 
293  std::string name() const
294  { return
295  std::string("finger") +
296  "_" + Model::name() +
297  "_" + Model::discretizationName() +
298  (this->model().enableGridAdaptation()?"_adaptive":"");
299  }
300 
304  static void registerParameters()
305  {
306  ParentType::registerParameters();
307 
308  EWOMS_REGISTER_PARAM(TypeTag, Scalar, InitialWaterSaturation,
309  "The initial saturation in the domain [] of the "
310  "wetting phase");
311  }
312 
316  void finishInit()
317  {
318  ParentType::finishInit();
319 
320  eps_ = 3e-6;
321 
322  temperature_ = 273.15 + 20; // -> 20°C
323 
324  FluidSystem::init();
325 
326  // parameters for the Van Genuchten law of the main imbibition
327  // and the main drainage curves.
328  micParams_.setVgAlpha(0.0037);
329  micParams_.setVgN(4.7);
330  micParams_.finalize();
331 
332  mdcParams_.setVgAlpha(0.0037);
333  mdcParams_.setVgN(4.7);
334  mdcParams_.finalize();
335 
336  // initialize the material parameter objects of the individual
337  // finite volumes, resize will resize the container to the number of elements
338  materialParams_.resize();
339 
340  for (auto it = materialParams_.begin(),
341  end = materialParams_.end(); it != end; ++it ) {
342  std::shared_ptr< MaterialLawParams >& materialParams = *it ;
343  if( ! materialParams )
344  {
345  materialParams.reset( new MaterialLawParams() );
346  materialParams->setMicParams(&micParams_);
347  materialParams->setMdcParams(&mdcParams_);
348  materialParams->setSwr(0.0);
349  materialParams->setSnr(0.1);
350  materialParams->finalize();
351  ParkerLenhard::reset(*materialParams);
352  }
353  }
354 
355  K_ = this->toDimMatrix_(4.6e-10);
356 
357  setupInitialFluidState_();
358  }
359 
363  void endTimeStep()
364  {
365 #ifndef NDEBUG
366  // checkConservativeness() does not include the effect of constraints, so we
367  // disable it for this problem...
368  //this->model().checkConservativeness();
369 
370  // Calculate storage terms
371  EqVector storage;
372  this->model().globalStorage(storage);
373 
374  // Write mass balance information for rank 0
375  if (this->gridView().comm().rank() == 0) {
376  std::cout << "Storage: " << storage << std::endl << std::flush;
377  }
378 #endif // NDEBUG
379 
380  // update the history of the hysteresis law
381  ElementContext elemCtx(this->simulator());
382 
383  for (const auto& elem : elements(this->gridView())) {
384  elemCtx.updateAll(elem);
385  size_t numDofs = elemCtx.numDof(/*timeIdx=*/0);
386  for (unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
387  {
388  MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, /*timeIdx=*/0 );
389  const auto& fs = elemCtx.intensiveQuantities(scvIdx, /*timeIdx=*/0).fluidState();
390  ParkerLenhard::update(materialParam, fs);
391  }
392  }
393  }
394 
396 
401 
405  template <class Context>
406  Scalar temperature(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
407  { return temperature_; }
408 
412  template <class Context>
413  const DimMatrix& intrinsicPermeability(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
414  { return K_; }
415 
419  template <class Context>
420  Scalar porosity(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
421  { return 0.4; }
422 
426  template <class Context>
427  MaterialLawParams& materialLawParams(const Context& context,
428  unsigned spaceIdx, unsigned timeIdx)
429  {
430  const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
431  assert(materialParams_[entity]);
432  return *materialParams_[entity];
433  }
434 
438  template <class Context>
439  const MaterialLawParams& materialLawParams(const Context& context,
440  unsigned spaceIdx, unsigned timeIdx) const
441  {
442  const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
443  assert(materialParams_[entity]);
444  return *materialParams_[entity];
445  }
446 
448 
453 
457  template <class Context>
458  void boundary(BoundaryRateVector& values, const Context& context,
459  unsigned spaceIdx, unsigned timeIdx) const
460  {
461  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
462 
463  if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
464  values.setNoFlow();
465  else {
466  assert(onUpperBoundary_(pos));
467 
468  values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
469  }
470 
471  // override the value for the liquid phase by forced
472  // imbibition of water on inlet boundary segments
473  if (onInlet_(pos)) {
474  values[contiWettingEqIdx] = -0.001; // [kg/(m^2 s)]
475  }
476  }
477 
479 
484 
488  template <class Context>
489  void initial(PrimaryVariables& values, const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
490  {
491  // assign the primary variables
492  values.assignNaive(initialFluidState_);
493  }
494 
498  template <class Context>
499  void constraints(Constraints& constraints, const Context& context,
500  unsigned spaceIdx, unsigned timeIdx) const
501  {
502  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
503 
504  if (onUpperBoundary_(pos) && !onInlet_(pos)) {
505  constraints.setActive(true);
506  constraints.assignNaive(initialFluidState_);
507  }
508  else if (onLowerBoundary_(pos)) {
509  constraints.setActive(true);
510  constraints.assignNaive(initialFluidState_);
511  }
512  }
513 
520  template <class Context>
521  void source(RateVector& rate, const Context& /*context*/,
522  unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
523  { rate = Scalar(0.0); }
525 
526 private:
527  bool onLeftBoundary_(const GlobalPosition& pos) const
528  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
529 
530  bool onRightBoundary_(const GlobalPosition& pos) const
531  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
532 
533  bool onLowerBoundary_(const GlobalPosition& pos) const
534  { return pos[1] < this->boundingBoxMin()[1] + eps_; }
535 
536  bool onUpperBoundary_(const GlobalPosition& pos) const
537  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
538 
539  bool onInlet_(const GlobalPosition& pos) const
540  {
541  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
542  Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
543 
544  if (!onUpperBoundary_(pos))
545  return false;
546 
547  Scalar xInject[] = { 0.25, 0.75 };
548  Scalar injectLen[] = { 0.1, 0.1 };
549  for (unsigned i = 0; i < sizeof(xInject) / sizeof(Scalar); ++i) {
550  if (xInject[i] - injectLen[i] / 2 < lambda
551  && lambda < xInject[i] + injectLen[i] / 2)
552  return true;
553  }
554  return false;
555  }
556 
557  void setupInitialFluidState_()
558  {
559  auto& fs = initialFluidState_;
560  fs.setPressure(wettingPhaseIdx, /*pressure=*/1e5);
561 
562  Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
563  fs.setSaturation(wettingPhaseIdx, Sw);
564  fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
565 
566  fs.setTemperature(temperature_);
567 
568  // set the absolute pressures
569  Scalar pn = 1e5;
570  fs.setPressure(nonWettingPhaseIdx, pn);
571  fs.setPressure(wettingPhaseIdx, pn);
572 
573  typename FluidSystem::template ParameterCache<Scalar> paramCache;
574  paramCache.updateAll(fs);
575  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
576  fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
577  fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
578  }
579 
580  }
581 
582  DimMatrix K_;
583 
584  typename MaterialLawParams::VanGenuchtenParams micParams_;
585  typename MaterialLawParams::VanGenuchtenParams mdcParams_;
586 
587  MaterialLawParamsContainer materialParams_;
588 
589  Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
590 
591  Scalar temperature_;
592  Scalar eps_;
593 };
594 
595 } // namespace Opm
596 
597 #endif
Two-phase problem featuring some gravity-driven saturation fingers.
Definition: fingerproblem.hh:214
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:413
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:458
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:420
MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx)
Definition: fingerproblem.hh:427
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:521
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:499
void finishInit()
Definition: fingerproblem.hh:316
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:406
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:489
void endTimeStep()
Definition: fingerproblem.hh:363
RestrictProlongOperator restrictProlongOperator()
Definition: fingerproblem.hh:285
static void registerParameters()
Definition: fingerproblem.hh:304
std::string name() const
Definition: fingerproblem.hh:293
FingerProblem(Simulator &simulator)
Definition: fingerproblem.hh:271
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:439
Definition: fingerproblem.hh:84
Definition: fingerproblem.hh:69