My Project
outflowproblem.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 */
27 #ifndef EWOMS_OUTFLOW_PROBLEM_HH
28 #define EWOMS_OUTFLOW_PROBLEM_HH
29 
30 #include <opm/models/pvs/pvsproperties.hh>
31 
32 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
33 #include <opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp>
34 #include <opm/material/common/Unused.hpp>
35 
36 #include <dune/grid/yaspgrid.hh>
37 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
38 
39 #include <dune/common/version.hh>
40 #include <dune/common/fvector.hh>
41 #include <dune/common/fmatrix.hh>
42 
43 namespace Opm {
44 template <class TypeTag>
45 class OutflowProblem;
46 }
47 
48 namespace Opm::Properties {
49 
50 namespace TTag {
51 
53 
54 } // namespace TTag
55 
56 // Set the grid type
57 template<class TypeTag>
58 struct Grid<TypeTag, TTag::OutflowBaseProblem> { using type = Dune::YaspGrid<2>; };
59 
60 // Set the problem property
61 template<class TypeTag>
62 struct Problem<TypeTag, TTag::OutflowBaseProblem> { using type = Opm::OutflowProblem<TypeTag>; };
63 
64 // Set fluid system
65 template<class TypeTag>
66 struct FluidSystem<TypeTag, TTag::OutflowBaseProblem>
67 {
68 private:
69  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
70 
71 public:
72  // Two-component single phase fluid system
73  using type = Opm::H2ON2LiquidPhaseFluidSystem<Scalar>;
74 };
75 
76 // Disable gravity
77 template<class TypeTag>
78 struct EnableGravity<TypeTag, TTag::OutflowBaseProblem> { static constexpr bool value = false; };
79 
80 // Also write mass fractions to the output
81 template<class TypeTag>
82 struct VtkWriteMassFractions<TypeTag, TTag::OutflowBaseProblem> { static constexpr bool value = true; };
83 
84 // The default for the end time of the simulation
85 template<class TypeTag>
86 struct EndTime<TypeTag, TTag::OutflowBaseProblem>
87 {
88  using type = GetPropType<TypeTag, Scalar>;
89  static constexpr type value = 100;
90 };
91 
92 // The default for the initial time step size of the simulation
93 template<class TypeTag>
94 struct InitialTimeStepSize<TypeTag, TTag::OutflowBaseProblem>
95 {
96  using type = GetPropType<TypeTag, Scalar>;
97  static constexpr type value = 1;
98 };
99 
100 // The default DGF file to load
101 template<class TypeTag>
102 struct GridFile<TypeTag, TTag::OutflowBaseProblem> { static constexpr auto value = "./data/outflow.dgf"; };
103 
104 } // namespace Opm::Properties
105 
106 namespace Opm {
124 template <class TypeTag>
125 class OutflowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
126 {
127  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
128 
129  using GridView = GetPropType<TypeTag, Properties::GridView>;
130  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
131  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
132  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
133  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
134  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
135  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
136  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
137  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
138 
139  // copy some indices for convenience
140  enum {
141  // Grid and world dimension
142  dim = GridView::dimension,
143  dimWorld = GridView::dimensionworld,
144 
145  numPhases = FluidSystem::numPhases,
146 
147  // component indices
148  H2OIdx = FluidSystem::H2OIdx,
149  N2Idx = FluidSystem::N2Idx
150  };
151 
152  using CoordScalar = typename GridView::ctype;
153  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
154 
155  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
156 
157 public:
161  OutflowProblem(Simulator& simulator)
162  : ParentType(simulator)
163  , eps_(1e-6)
164  { }
165 
169  void finishInit()
170  {
171  ParentType::finishInit();
172 
173  temperature_ = 273.15 + 20;
174  FluidSystem::init(/*minT=*/temperature_ - 1, /*maxT=*/temperature_ + 2,
175  /*numT=*/3,
176  /*minp=*/0.8e5, /*maxp=*/2.5e5, /*nump=*/500);
177 
178  // set parameters of porous medium
179  perm_ = this->toDimMatrix_(1e-10);
180  porosity_ = 0.4;
181  tortuosity_ = 0.28;
182  }
183 
188 
192  std::string name() const
193  { return "outflow"; }
194 
198  void endTimeStep()
199  {
200 #ifndef NDEBUG
201  this->model().checkConservativeness();
202 
203  // Calculate storage terms
204  EqVector storage;
205  this->model().globalStorage(storage);
206 
207  // Write mass balance information for rank 0
208  if (this->gridView().comm().rank() == 0) {
209  std::cout << "Storage: " << storage << std::endl << std::flush;
210  }
211 #endif // NDEBUG
212  }
213 
219  template <class Context>
220  Scalar temperature(const Context& context OPM_UNUSED,
221  unsigned spaceIdx OPM_UNUSED,
222  unsigned timeIdx OPM_UNUSED) const
223  { return temperature_; } // in [K]
224 
230  template <class Context>
231  const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED,
232  unsigned spaceIdx OPM_UNUSED,
233  unsigned timeIdx OPM_UNUSED) const
234  { return perm_; }
235 
241  template <class Context>
242  Scalar porosity(const Context& context OPM_UNUSED,
243  unsigned spaceIdx OPM_UNUSED,
244  unsigned timeIdx OPM_UNUSED) const
245  { return porosity_; }
246 
247 #if 0
252  template <class Context>
253  Scalar tortuosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
254  { return tortuosity_; }
255 
260  template <class Context>
261  Scalar dispersivity(const Context& context,
262  unsigned spaceIdx, unsigned timeIdx) const
263  { return 0; }
264 #endif
265 
267 
272 
276  template <class Context>
277  void boundary(BoundaryRateVector& values, const Context& context,
278  unsigned spaceIdx, unsigned timeIdx) const
279  {
280  const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
281 
282  if (onLeftBoundary_(globalPos)) {
283  Opm::CompositionalFluidState<Scalar, FluidSystem,
284  /*storeEnthalpy=*/false> fs;
285  initialFluidState_(fs, context, spaceIdx, timeIdx);
286  fs.setPressure(/*phaseIdx=*/0, fs.pressure(/*phaseIdx=*/0) + 1e5);
287 
288  Scalar xlN2 = 2e-4;
289  fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, xlN2);
290  fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1 - xlN2);
291 
292  typename FluidSystem::template ParameterCache<Scalar> paramCache;
293  paramCache.updateAll(fs);
294  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
295  fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
296  fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
297  }
298 
299  // impose an freeflow boundary condition
300  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
301  }
302  else if (onRightBoundary_(globalPos)) {
303  Opm::CompositionalFluidState<Scalar, FluidSystem,
304  /*storeEnthalpy=*/false> fs;
305  initialFluidState_(fs, context, spaceIdx, timeIdx);
306 
307  // impose an outflow boundary condition
308  values.setOutFlow(context, spaceIdx, timeIdx, fs);
309  }
310  else
311  // no flow on top and bottom
312  values.setNoFlow();
313  }
314 
316 
321 
325  template <class Context>
326  void initial(PrimaryVariables& values,
327  const Context& context,
328  unsigned spaceIdx,
329  unsigned timeIdx) const
330  {
331  Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/false> fs;
332  initialFluidState_(fs, context, spaceIdx, timeIdx);
333 
334  values.assignNaive(fs);
335  }
336 
343  template <class Context>
344  void source(RateVector& rate,
345  const Context& context OPM_UNUSED,
346  unsigned spaceIdx OPM_UNUSED,
347  unsigned timeIdx OPM_UNUSED) const
348  { rate = Scalar(0.0); }
349 
351 
352 private:
353  bool onLeftBoundary_(const GlobalPosition& pos) const
354  { return pos[0] < eps_; }
355 
356  bool onRightBoundary_(const GlobalPosition& pos) const
357  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
358 
359  template <class FluidState, class Context>
360  void initialFluidState_(FluidState& fs, const Context& context,
361  unsigned spaceIdx, unsigned timeIdx) const
362  {
363  Scalar T = temperature(context, spaceIdx, timeIdx);
364  // Scalar rho = FluidSystem::H2O::liquidDensity(T, /*pressure=*/1.5e5);
365  // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
366  // this->boundingBoxMax()[dim - 1];
367  // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
368  // this->boundingBoxMax()[dim - 1];
369 
370  fs.setSaturation(/*phaseIdx=*/0, 1.0);
371  fs.setPressure(/*phaseIdx=*/0, 1e5 /* + rho*z */);
372  fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1.0);
373  fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, 0);
374  fs.setTemperature(T);
375 
376  typename FluidSystem::template ParameterCache<Scalar> paramCache;
377  paramCache.updateAll(fs);
378  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
379  fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
380  fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
381  }
382  }
383 
384  const Scalar eps_;
385 
386  MaterialLawParams materialParams_;
387  DimMatrix perm_;
388  Scalar temperature_;
389  Scalar porosity_;
390  Scalar tortuosity_;
391 };
392 } // namespace Opm
393 
394 #endif
Problem where dissolved nitrogen is transported with the water phase from the left side to the right.
Definition: outflowproblem.hh:126
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: outflowproblem.hh:326
void finishInit()
Definition: outflowproblem.hh:169
void endTimeStep()
Definition: outflowproblem.hh:198
OutflowProblem(Simulator &simulator)
Definition: outflowproblem.hh:161
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: outflowproblem.hh:277
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:220
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:231
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:344
std::string name() const
Definition: outflowproblem.hh:192
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:242
Definition: outflowproblem.hh:52