All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
diffusionproblem.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_POWER_INJECTION_PROBLEM_HH
29 #define EWOMS_POWER_INJECTION_PROBLEM_HH
30 
32 
34 
35 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
36 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
37 #include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>
38 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
39 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
40 #include <opm/common/Unused.hpp>
41 
42 #include <dune/grid/yaspgrid.hh>
43 #include <dune/common/version.hh>
44 #include <dune/common/fvector.hh>
45 #include <dune/common/fmatrix.hh>
46 
47 #include <sstream>
48 #include <string>
49 
50 namespace Ewoms {
51 template <class TypeTag>
53 }
54 
55 namespace Ewoms {
56 namespace Properties {
57 
58 NEW_TYPE_TAG(DiffusionBaseProblem);
59 
60 // Set the grid implementation to be used
61 SET_TYPE_PROP(DiffusionBaseProblem, Grid, Dune::YaspGrid</*dim=*/1>);
62 
63 // set the GridManager property
64 SET_TYPE_PROP(DiffusionBaseProblem, GridManager, Ewoms::CubeGridManager<TypeTag>);
65 
66 // Set the problem property
67 SET_TYPE_PROP(DiffusionBaseProblem, Problem, Ewoms::DiffusionProblem<TypeTag>);
68 
69 // Set the fluid system
70 SET_PROP(DiffusionBaseProblem, FluidSystem)
71 {
72 private:
73  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
74 
75 public:
76  typedef Opm::FluidSystems::H2ON2<Scalar> type;
77 };
78 
79 // Set the material Law
80 SET_PROP(DiffusionBaseProblem, MaterialLaw)
81 {
82 private:
83  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
84  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
85 
86  static_assert(FluidSystem::numPhases == 2,
87  "A fluid system with two phases is required "
88  "for this problem!");
89 
90  typedef Opm::TwoPhaseMaterialTraits<Scalar,
91  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
92  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>
93  Traits;
94 
95 public:
96  typedef Opm::LinearMaterial<Traits> type;
97 };
98 
99 // Enable molecular diffusion for this problem
100 SET_BOOL_PROP(DiffusionBaseProblem, EnableDiffusion, true);
101 
102 // Disable gravity
103 SET_BOOL_PROP(DiffusionBaseProblem, EnableGravity, false);
104 
105 // define the properties specific for the diffusion problem
106 SET_SCALAR_PROP(DiffusionBaseProblem, DomainSizeX, 1.0);
107 SET_SCALAR_PROP(DiffusionBaseProblem, DomainSizeY, 1.0);
108 SET_SCALAR_PROP(DiffusionBaseProblem, DomainSizeZ, 1.0);
109 
110 SET_INT_PROP(DiffusionBaseProblem, CellsX, 250);
111 SET_INT_PROP(DiffusionBaseProblem, CellsY, 1);
112 SET_INT_PROP(DiffusionBaseProblem, CellsZ, 1);
113 
114 // The default for the end time of the simulation
115 SET_SCALAR_PROP(DiffusionBaseProblem, EndTime, 1e6);
116 
117 // The default for the initial time step size of the simulation
118 SET_SCALAR_PROP(DiffusionBaseProblem, InitialTimeStepSize, 1000);
119 }
120 } // namespace Ewoms, Properties
121 
122 namespace Ewoms {
133 template <class TypeTag>
134 class DiffusionProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
135 {
136  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
137 
138  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
139  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
140  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
141  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
142  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
143  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
144 
145  enum {
146  // number of phases
147  numPhases = FluidSystem::numPhases,
148 
149  // phase indices
150  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
151  gasPhaseIdx = FluidSystem::gasPhaseIdx,
152 
153  // component indices
154  H2OIdx = FluidSystem::H2OIdx,
155  N2Idx = FluidSystem::N2Idx,
156 
157  // Grid and world dimension
158  dim = GridView::dimension,
159  dimWorld = GridView::dimensionworld
160  };
161 
162  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
163  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
164  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
165 
166  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
167  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
168 
169  typedef typename GridView::ctype CoordScalar;
170  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
171 
172  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
173 
174 public:
179  : ParentType(simulator)
180  { }
181 
185  void finishInit()
186  {
187  ParentType::finishInit();
188 
189  FluidSystem::init();
190 
191  temperature_ = 273.15 + 20.0;
192 
193  materialParams_.finalize();
194 
195  K_ = this->toDimMatrix_(1e-12); // [m^2]
196 
197  setupInitialFluidStates_();
198  }
199 
203 
208  std::string name() const
209  { return std::string("diffusion_") + Model::name(); }
210 
214  void endTimeStep()
215  {
216 #ifndef NDEBUG
217  this->model().checkConservativeness();
218 
219  // Calculate storage terms
220  EqVector storage;
221  this->model().globalStorage(storage);
222 
223  // Write mass balance information for rank 0
224  if (this->gridView().comm().rank() == 0) {
225  std::cout << "Storage: " << storage << std::endl << std::flush;
226  }
227 #endif // NDEBUG
228  }
229 
231 
235 
240  template <class Context>
241  const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED,
242  unsigned spaceIdx OPM_UNUSED,
243  unsigned timeIdx OPM_UNUSED) const
244  { return K_; }
245 
249  template <class Context>
250  Scalar porosity(const Context& context OPM_UNUSED,
251  unsigned spaceIdx OPM_UNUSED,
252  unsigned timeIdx OPM_UNUSED) const
253  { return 0.35; }
254 
258  template <class Context>
259  const MaterialLawParams&
260  materialLawParams(const Context& context OPM_UNUSED,
261  unsigned spaceIdx OPM_UNUSED,
262  unsigned timeIdx OPM_UNUSED) const
263  { return materialParams_; }
264 
268  template <class Context>
269  Scalar temperature(const Context& context OPM_UNUSED,
270  unsigned spaceIdx OPM_UNUSED,
271  unsigned timeIdx OPM_UNUSED) const
272  { return temperature_; }
273 
275 
279 
286  template <class Context>
287  void boundary(BoundaryRateVector& values,
288  const Context& context OPM_UNUSED,
289  unsigned spaceIdx OPM_UNUSED,
290  unsigned timeIdx OPM_UNUSED) const
291  { values.setNoFlow(); }
292 
294 
298 
303  template <class Context>
304  void initial(PrimaryVariables& values,
305  const Context& context,
306  unsigned spaceIdx,
307  unsigned timeIdx) const
308  {
309  const auto& pos = context.pos(spaceIdx, timeIdx);
310  if (onLeftSide_(pos))
311  values.assignNaive(leftInitialFluidState_);
312  else
313  values.assignNaive(rightInitialFluidState_);
314  }
315 
322  template <class Context>
323  void source(RateVector& rate,
324  const Context& context OPM_UNUSED,
325  unsigned spaceIdx OPM_UNUSED,
326  unsigned timeIdx OPM_UNUSED) const
327  { rate = Scalar(0.0); }
328 
330 
331 private:
332  bool onLeftSide_(const GlobalPosition& pos) const
333  { return pos[0] < (this->boundingBoxMin()[0] + this->boundingBoxMax()[0]) / 2; }
334 
335  void setupInitialFluidStates_()
336  {
337  // create the initial fluid state for the left half of the domain
338  leftInitialFluidState_.setTemperature(temperature_);
339 
340  Scalar Sl = 0.0;
341  leftInitialFluidState_.setSaturation(liquidPhaseIdx, Sl);
342  leftInitialFluidState_.setSaturation(gasPhaseIdx, 1 - Sl);
343 
344  Scalar p = 1e5;
345  leftInitialFluidState_.setPressure(liquidPhaseIdx, p);
346  leftInitialFluidState_.setPressure(gasPhaseIdx, p);
347 
348  Scalar xH2O = 0.01;
349  leftInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
350  leftInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
351 
352  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
353  typename FluidSystem::template ParameterCache<Scalar> paramCache;
354  CFRP::solve(leftInitialFluidState_, paramCache, gasPhaseIdx,
355  /*setViscosity=*/false, /*setEnthalpy=*/false);
356 
357  // create the initial fluid state for the right half of the domain
358  rightInitialFluidState_.assign(leftInitialFluidState_);
359  xH2O = 0.0;
360  rightInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
361  rightInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
362  CFRP::solve(rightInitialFluidState_, paramCache, gasPhaseIdx,
363  /*setViscosity=*/false, /*setEnthalpy=*/false);
364  }
365 
366  DimMatrix K_;
367  MaterialLawParams materialParams_;
368 
369  Opm::CompositionalFluidState<Scalar, FluidSystem> leftInitialFluidState_;
370  Opm::CompositionalFluidState<Scalar, FluidSystem> rightInitialFluidState_;
371  Scalar temperature_;
372 };
373 
374 } // namespace Ewoms
375 
376 #endif
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
const MaterialLawParams & materialLawParams(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: diffusionproblem.hh:260
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: diffusionproblem.hh:323
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: diffusionproblem.hh:185
void endTimeStep()
Called by the simulator after each time integration.
Definition: diffusionproblem.hh:214
std::string name() const
The problem name.
Definition: diffusionproblem.hh:208
1D problem which is driven by molecular diffusion.
Definition: diffusionproblem.hh:52
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
Declares the properties required for the NCP compositional multi-phase model.
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: diffusionproblem.hh:304
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: diffusionproblem.hh:269
Provides a grid manager which a regular grid made of quadrilaterals.
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar porosity(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: diffusionproblem.hh:250
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
const DimMatrix & intrinsicPermeability(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: diffusionproblem.hh:241
Provides a grid manager which a regular grid made of quadrilaterals.
Definition: cubegridmanager.hh:65
void boundary(BoundaryRateVector &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the boundary conditions for a boundary segment.
Definition: diffusionproblem.hh:287
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394