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