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  auto elemIt = this->gridView().template begin<0>();
384  const auto& elemEndIt = this->gridView().template end<0>();
385  for (; elemIt != elemEndIt; ++elemIt) {
386  const auto& elem = *elemIt;
387  elemCtx.updateAll( elem );
388  size_t numDofs = elemCtx.numDof(/*timeIdx=*/0);
389  for (unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
390  {
391  MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, /*timeIdx=*/0 );
392  const auto& fs = elemCtx.intensiveQuantities(scvIdx, /*timeIdx=*/0).fluidState();
393  ParkerLenhard::update(materialParam, fs);
394  }
395  }
396  }
397 
399 
404 
408  template <class Context>
409  Scalar temperature(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
410  { return temperature_; }
411 
415  template <class Context>
416  const DimMatrix& intrinsicPermeability(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
417  { return K_; }
418 
422  template <class Context>
423  Scalar porosity(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
424  { return 0.4; }
425 
429  template <class Context>
430  MaterialLawParams& materialLawParams(const Context& context,
431  unsigned spaceIdx, unsigned timeIdx)
432  {
433  const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
434  assert(materialParams_[entity]);
435  return *materialParams_[entity];
436  }
437 
441  template <class Context>
442  const MaterialLawParams& materialLawParams(const Context& context,
443  unsigned spaceIdx, unsigned timeIdx) const
444  {
445  const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
446  assert(materialParams_[entity]);
447  return *materialParams_[entity];
448  }
449 
451 
456 
460  template <class Context>
461  void boundary(BoundaryRateVector& values, const Context& context,
462  unsigned spaceIdx, unsigned timeIdx) const
463  {
464  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
465 
466  if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
467  values.setNoFlow();
468  else {
469  assert(onUpperBoundary_(pos));
470 
471  values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
472  }
473 
474  // override the value for the liquid phase by forced
475  // imbibition of water on inlet boundary segments
476  if (onInlet_(pos)) {
477  values[contiWettingEqIdx] = -0.001; // [kg/(m^2 s)]
478  }
479  }
480 
482 
487 
491  template <class Context>
492  void initial(PrimaryVariables& values, const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
493  {
494  // assign the primary variables
495  values.assignNaive(initialFluidState_);
496  }
497 
501  template <class Context>
502  void constraints(Constraints& constraints, const Context& context,
503  unsigned spaceIdx, unsigned timeIdx) const
504  {
505  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
506 
507  if (onUpperBoundary_(pos) && !onInlet_(pos)) {
508  constraints.setActive(true);
509  constraints.assignNaive(initialFluidState_);
510  }
511  else if (onLowerBoundary_(pos)) {
512  constraints.setActive(true);
513  constraints.assignNaive(initialFluidState_);
514  }
515  }
516 
523  template <class Context>
524  void source(RateVector& rate, const Context& /*context*/,
525  unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
526  { rate = Scalar(0.0); }
528 
529 private:
530  bool onLeftBoundary_(const GlobalPosition& pos) const
531  { return pos[0] < this->boundingBoxMin()[0] + eps_; }
532 
533  bool onRightBoundary_(const GlobalPosition& pos) const
534  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
535 
536  bool onLowerBoundary_(const GlobalPosition& pos) const
537  { return pos[1] < this->boundingBoxMin()[1] + eps_; }
538 
539  bool onUpperBoundary_(const GlobalPosition& pos) const
540  { return pos[1] > this->boundingBoxMax()[1] - eps_; }
541 
542  bool onInlet_(const GlobalPosition& pos) const
543  {
544  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
545  Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
546 
547  if (!onUpperBoundary_(pos))
548  return false;
549 
550  Scalar xInject[] = { 0.25, 0.75 };
551  Scalar injectLen[] = { 0.1, 0.1 };
552  for (unsigned i = 0; i < sizeof(xInject) / sizeof(Scalar); ++i) {
553  if (xInject[i] - injectLen[i] / 2 < lambda
554  && lambda < xInject[i] + injectLen[i] / 2)
555  return true;
556  }
557  return false;
558  }
559 
560  void setupInitialFluidState_()
561  {
562  auto& fs = initialFluidState_;
563  fs.setPressure(wettingPhaseIdx, /*pressure=*/1e5);
564 
565  Scalar Sw = EWOMS_GET_PARAM(TypeTag, Scalar, InitialWaterSaturation);
566  fs.setSaturation(wettingPhaseIdx, Sw);
567  fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
568 
569  fs.setTemperature(temperature_);
570 
571  // set the absolute pressures
572  Scalar pn = 1e5;
573  fs.setPressure(nonWettingPhaseIdx, pn);
574  fs.setPressure(wettingPhaseIdx, pn);
575 
576  typename FluidSystem::template ParameterCache<Scalar> paramCache;
577  paramCache.updateAll(fs);
578  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
579  fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
580  fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
581  }
582 
583  }
584 
585  DimMatrix K_;
586 
587  typename MaterialLawParams::VanGenuchtenParams micParams_;
588  typename MaterialLawParams::VanGenuchtenParams mdcParams_;
589 
590  MaterialLawParamsContainer materialParams_;
591 
592  Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
593 
594  Scalar temperature_;
595  Scalar eps_;
596 };
597 
598 } // namespace Opm
599 
600 #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:416
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:461
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:423
MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx)
Definition: fingerproblem.hh:430
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:524
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: fingerproblem.hh:502
void finishInit()
Definition: fingerproblem.hh:316
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:409
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition: fingerproblem.hh:492
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:442
Definition: fingerproblem.hh:84
Definition: fingerproblem.hh:69