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 
31 
32 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
33 #include <opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp>
34 #include <opm/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 Ewoms {
44 template <class TypeTag>
46 }
47 
48 namespace Ewoms {
49 namespace Properties {
50 NEW_TYPE_TAG(OutflowBaseProblem);
51 
52 // Set the grid type
53 SET_TYPE_PROP(OutflowBaseProblem, Grid, Dune::YaspGrid<2>);
54 
55 // Set the problem property
56 SET_TYPE_PROP(OutflowBaseProblem, Problem, Ewoms::OutflowProblem<TypeTag>);
57 
58 // Set fluid system
59 SET_PROP(OutflowBaseProblem, FluidSystem)
60 {
61 private:
62  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
63 
64 public:
65  // Two-component single phase fluid system
66  typedef Opm::FluidSystems::H2ON2LiquidPhase<Scalar> type;
67 };
68 
69 // Disable gravity
70 SET_BOOL_PROP(OutflowBaseProblem, EnableGravity, false);
71 
72 // Also write mass fractions to the output
73 SET_BOOL_PROP(OutflowBaseProblem, VtkWriteMassFractions, true);
74 
75 // The default for the end time of the simulation
76 SET_SCALAR_PROP(OutflowBaseProblem, EndTime, 100);
77 
78 // The default for the initial time step size of the simulation
79 SET_SCALAR_PROP(OutflowBaseProblem, InitialTimeStepSize, 1);
80 
81 // The default DGF file to load
82 SET_STRING_PROP(OutflowBaseProblem, GridFile, "./data/outflow.dgf");
83 } // namespace Properties
84 } // namespace Ewoms
85 
86 namespace Ewoms {
104 template <class TypeTag>
105 class OutflowProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
106 {
107  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
108 
109  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
110  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
111  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
112  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
113  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
114  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
115  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
116  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
117  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
118 
119  // copy some indices for convenience
120  enum {
121  // Grid and world dimension
122  dim = GridView::dimension,
123  dimWorld = GridView::dimensionworld,
124 
125  // component indices
126  H2OIdx = FluidSystem::H2OIdx,
127  N2Idx = FluidSystem::N2Idx
128  };
129 
130  typedef typename GridView::ctype CoordScalar;
131  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
132 
133  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
134 
135 public:
140  : ParentType(simulator)
141  , eps_(1e-6)
142  { }
143 
147  void finishInit()
148  {
149  ParentType::finishInit();
150 
151  temperature_ = 273.15 + 20;
152  FluidSystem::init(/*minT=*/temperature_ - 1, /*maxT=*/temperature_ + 2,
153  /*numT=*/3,
154  /*minp=*/0.8e5, /*maxp=*/2.5e5, /*nump=*/500);
155 
156  // set parameters of porous medium
157  perm_ = this->toDimMatrix_(1e-10);
158  porosity_ = 0.4;
159  tortuosity_ = 0.28;
160  }
161 
165 
170  std::string name() const
171  { return "outflow"; }
172 
176  void endTimeStep()
177  {
178 #ifndef NDEBUG
179  this->model().checkConservativeness();
180 
181  // Calculate storage terms
182  EqVector storage;
183  this->model().globalStorage(storage);
184 
185  // Write mass balance information for rank 0
186  if (this->gridView().comm().rank() == 0) {
187  std::cout << "Storage: " << storage << std::endl << std::flush;
188  }
189 #endif // NDEBUG
190  }
191 
197  template <class Context>
198  Scalar temperature(const Context& context OPM_UNUSED,
199  unsigned spaceIdx OPM_UNUSED,
200  unsigned timeIdx OPM_UNUSED) const
201  { return temperature_; } // in [K]
202 
208  template <class Context>
209  const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED,
210  unsigned spaceIdx OPM_UNUSED,
211  unsigned timeIdx OPM_UNUSED) const
212  { return perm_; }
213 
219  template <class Context>
220  Scalar porosity(const Context& context OPM_UNUSED,
221  unsigned spaceIdx OPM_UNUSED,
222  unsigned timeIdx OPM_UNUSED) const
223  { return porosity_; }
224 
225 #if 0
226 
230  template <class Context>
231  Scalar tortuosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
232  { return tortuosity_; }
233 
238  template <class Context>
239  Scalar dispersivity(const Context& context,
240  unsigned spaceIdx, unsigned timeIdx) const
241  { return 0; }
242 #endif
243 
245 
249 
254  template <class Context>
255  void boundary(BoundaryRateVector& values, const Context& context,
256  unsigned spaceIdx, unsigned timeIdx) const
257  {
258  const GlobalPosition& globalPos = context.pos(spaceIdx, timeIdx);
259 
260  if (onLeftBoundary_(globalPos)) {
261  Opm::CompositionalFluidState<Scalar, FluidSystem,
262  /*storeEnthalpy=*/false> fs;
263  initialFluidState_(fs, context, spaceIdx, timeIdx);
264  fs.setPressure(/*phaseIdx=*/0, fs.pressure(/*phaseIdx=*/0) + 1e5);
265 
266  Scalar xlN2 = 2e-4;
267  fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, xlN2);
268  fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1 - xlN2);
269 
270  // impose an freeflow boundary condition
271  values.setFreeFlow(context, spaceIdx, timeIdx, fs);
272  }
273  else if (onRightBoundary_(globalPos)) {
274  Opm::CompositionalFluidState<Scalar, FluidSystem,
275  /*storeEnthalpy=*/false> fs;
276  initialFluidState_(fs, context, spaceIdx, timeIdx);
277 
278  // impose an outflow boundary condition
279  values.setOutFlow(context, spaceIdx, timeIdx, fs);
280  }
281  else
282  // no flow on top and bottom
283  values.setNoFlow();
284  }
285 
287 
291 
296  template <class Context>
297  void initial(PrimaryVariables& values,
298  const Context& context,
299  unsigned spaceIdx,
300  unsigned timeIdx) const
301  {
302  Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/false> fs;
303  initialFluidState_(fs, context, spaceIdx, timeIdx);
304 
305  values.assignNaive(fs);
306  }
307 
314  template <class Context>
315  void source(RateVector& rate,
316  const Context& context OPM_UNUSED,
317  unsigned spaceIdx OPM_UNUSED,
318  unsigned timeIdx OPM_UNUSED) const
319  { rate = Scalar(0.0); }
320 
322 
323 private:
324  bool onLeftBoundary_(const GlobalPosition& pos) const
325  { return pos[0] < eps_; }
326 
327  bool onRightBoundary_(const GlobalPosition& pos) const
328  { return pos[0] > this->boundingBoxMax()[0] - eps_; }
329 
330  template <class FluidState, class Context>
331  void initialFluidState_(FluidState& fs, const Context& context,
332  unsigned spaceIdx, unsigned timeIdx) const
333  {
334  Scalar T = temperature(context, spaceIdx, timeIdx);
335  // Scalar rho = FluidSystem::H2O::liquidDensity(T, /*pressure=*/1.5e5);
336  // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
337  // this->boundingBoxMax()[dim - 1];
338  // Scalar z = context.pos(spaceIdx, timeIdx)[dim - 1] -
339  // this->boundingBoxMax()[dim - 1];
340 
341  fs.setSaturation(/*phaseIdx=*/0, 1.0);
342  fs.setPressure(/*phaseIdx=*/0, 1e5 /* + rho*z */);
343  fs.setMoleFraction(/*phaseIdx=*/0, H2OIdx, 1.0);
344  fs.setMoleFraction(/*phaseIdx=*/0, N2Idx, 0);
345  fs.setTemperature(T);
346  }
347 
348  const Scalar eps_;
349 
350  MaterialLawParams materialParams_;
351  DimMatrix perm_;
352  Scalar temperature_;
353  Scalar porosity_;
354  Scalar tortuosity_;
355 };
356 } // namespace Ewoms
357 
358 #endif
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:220
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:209
Definition: baseauxiliarymodule.hh:37
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: outflowproblem.hh:315
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: outflowproblem.hh:198
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: outflowproblem.hh:147
std::string name() const
The problem name.
Definition: outflowproblem.hh:170
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
void endTimeStep()
Called by the simulator after each time integration.
Definition: outflowproblem.hh:176
Declares the properties required for the compositional multi-phase primary variable switching model...
OutflowProblem(Simulator &simulator)
Definition: outflowproblem.hh:139
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: outflowproblem.hh:297
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the boundary conditions for a boundary segment.
Definition: outflowproblem.hh:255
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
Problem where dissolved nitrogen is transported with the water phase from the left side to the right...
Definition: outflowproblem.hh:45
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394