My Project
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 
31 #include <opm/models/ncp/ncpproperties.hh>
32 
33 #include <opm/models/io/cubegridvanguard.hh>
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 
41 #include <dune/grid/yaspgrid.hh>
42 #include <dune/common/version.hh>
43 #include <dune/common/fvector.hh>
44 #include <dune/common/fmatrix.hh>
45 
46 #include <sstream>
47 #include <string>
48 
49 namespace Opm {
50 template <class TypeTag>
51 class DiffusionProblem;
52 }
53 
54 namespace Opm::Properties {
55 
56 namespace TTag {
57 
59 
60 } // namespace TTag
61 
62 // Set the grid implementation to be used
63 template<class TypeTag>
64 struct Grid<TypeTag, TTag::DiffusionBaseProblem> { using type = Dune::YaspGrid</*dim=*/1>; };
65 
66 // set the Vanguard property
67 template<class TypeTag>
68 struct Vanguard<TypeTag, TTag::DiffusionBaseProblem> { using type = Opm::CubeGridVanguard<TypeTag>; };
69 
70 // Set the problem property
71 template<class TypeTag>
72 struct Problem<TypeTag, TTag::DiffusionBaseProblem> { using type = Opm::DiffusionProblem<TypeTag>; };
73 
74 // Set the fluid system
75 template<class TypeTag>
76 struct FluidSystem<TypeTag, TTag::DiffusionBaseProblem>
77 {
78 private:
79  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
80 
81 public:
82  using type = Opm::H2ON2FluidSystem<Scalar>;
83 };
84 
85 // Set the material Law
86 template<class TypeTag>
87 struct MaterialLaw<TypeTag, TTag::DiffusionBaseProblem>
88 {
89 private:
90  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
92 
93  static_assert(FluidSystem::numPhases == 2,
94  "A fluid system with two phases is required "
95  "for this problem!");
96 
97  using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
98  /*wettingPhaseIdx=*/FluidSystem::liquidPhaseIdx,
99  /*nonWettingPhaseIdx=*/FluidSystem::gasPhaseIdx>;
100 
101 public:
102  using type = Opm::LinearMaterial<Traits>;
103 };
104 
105 // Enable molecular diffusion for this problem
106 template<class TypeTag>
107 struct EnableDiffusion<TypeTag, TTag::DiffusionBaseProblem> { static constexpr bool value = true; };
108 
109 // Disable gravity
110 template<class TypeTag>
111 struct EnableGravity<TypeTag, TTag::DiffusionBaseProblem> { static constexpr bool value = false; };
112 
113 // define the properties specific for the diffusion problem
114 template<class TypeTag>
115 struct DomainSizeX<TypeTag, TTag::DiffusionBaseProblem>
116 {
117  using type = GetPropType<TypeTag, Scalar>;
118  static constexpr type value = 1.0;
119 };
120 template<class TypeTag>
121 struct DomainSizeY<TypeTag, TTag::DiffusionBaseProblem>
122 {
123  using type = GetPropType<TypeTag, Scalar>;
124  static constexpr type value = 1.0;
125 };
126 template<class TypeTag>
127 struct DomainSizeZ<TypeTag, TTag::DiffusionBaseProblem>
128 {
129  using type = GetPropType<TypeTag, Scalar>;
130  static constexpr type value = 1.0;
131 };
132 
133 template<class TypeTag>
134 struct CellsX<TypeTag, TTag::DiffusionBaseProblem> { static constexpr int value = 250; };
135 template<class TypeTag>
136 struct CellsY<TypeTag, TTag::DiffusionBaseProblem> { static constexpr int value = 1; };
137 template<class TypeTag>
138 struct CellsZ<TypeTag, TTag::DiffusionBaseProblem> { static constexpr int value = 1; };
139 
140 // The default for the end time of the simulation
141 template<class TypeTag>
142 struct EndTime<TypeTag, TTag::DiffusionBaseProblem>
143 {
144  using type = GetPropType<TypeTag, Scalar>;
145  static constexpr type value = 1e6;
146 };
147 
148 // The default for the initial time step size of the simulation
149 template<class TypeTag>
150 struct InitialTimeStepSize<TypeTag, TTag::DiffusionBaseProblem>
151 {
152  using type = GetPropType<TypeTag, Scalar>;
153  static constexpr type value = 1000;
154 };
155 
156 } // namespace Opm::Properties
157 
158 namespace Opm {
169 template <class TypeTag>
170 class DiffusionProblem : public GetPropType<TypeTag, Properties::BaseProblem>
171 {
172  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
173 
174  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
175  using GridView = GetPropType<TypeTag, Properties::GridView>;
176  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
177  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
178  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
179  using Model = GetPropType<TypeTag, Properties::Model>;
180 
181  enum {
182  // number of phases
183  numPhases = FluidSystem::numPhases,
184 
185  // phase indices
186  liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
187  gasPhaseIdx = FluidSystem::gasPhaseIdx,
188 
189  // component indices
190  H2OIdx = FluidSystem::H2OIdx,
191  N2Idx = FluidSystem::N2Idx,
192 
193  // Grid and world dimension
194  dim = GridView::dimension,
195  dimWorld = GridView::dimensionworld
196  };
197 
198  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
199  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
200  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
201 
202  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
203  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
204 
205  using CoordScalar = typename GridView::ctype;
206  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
207 
208  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
209 
210 public:
214  DiffusionProblem(Simulator& simulator)
215  : ParentType(simulator)
216  { }
217 
221  void finishInit()
222  {
223  ParentType::finishInit();
224 
225  FluidSystem::init();
226 
227  temperature_ = 273.15 + 20.0;
228 
229  materialParams_.finalize();
230 
231  K_ = this->toDimMatrix_(1e-12); // [m^2]
232 
233  setupInitialFluidStates_();
234  }
235 
240 
244  std::string name() const
245  { return std::string("diffusion_") + Model::name(); }
246 
250  void endTimeStep()
251  {
252 #ifndef NDEBUG
253  this->model().checkConservativeness();
254 
255  // Calculate storage terms
256  EqVector storage;
257  this->model().globalStorage(storage);
258 
259  // Write mass balance information for rank 0
260  if (this->gridView().comm().rank() == 0) {
261  std::cout << "Storage: " << storage << std::endl << std::flush;
262  }
263 #endif // NDEBUG
264  }
265 
267 
272 
276  template <class Context>
277  const DimMatrix& intrinsicPermeability(const Context& /*context*/,
278  unsigned /*spaceIdx*/,
279  unsigned /*timeIdx*/) const
280  { return K_; }
281 
285  template <class Context>
286  Scalar porosity(const Context& /*context*/,
287  unsigned /*spaceIdx*/,
288  unsigned /*timeIdx*/) const
289  { return 0.35; }
290 
294  template <class Context>
295  const MaterialLawParams&
296  materialLawParams(const Context& /*context*/,
297  unsigned /*spaceIdx*/,
298  unsigned /*timeIdx*/) const
299  { return materialParams_; }
300 
304  template <class Context>
305  Scalar temperature(const Context& /*context*/,
306  unsigned /*spaceIdx*/,
307  unsigned /*timeIdx*/) const
308  { return temperature_; }
309 
311 
316 
322  template <class Context>
323  void boundary(BoundaryRateVector& values,
324  const Context& /*context*/,
325  unsigned /*spaceIdx*/,
326  unsigned /*timeIdx*/) const
327  { values.setNoFlow(); }
328 
330 
335 
339  template <class Context>
340  void initial(PrimaryVariables& values,
341  const Context& context,
342  unsigned spaceIdx,
343  unsigned timeIdx) const
344  {
345  const auto& pos = context.pos(spaceIdx, timeIdx);
346  if (onLeftSide_(pos))
347  values.assignNaive(leftInitialFluidState_);
348  else
349  values.assignNaive(rightInitialFluidState_);
350  }
351 
358  template <class Context>
359  void source(RateVector& rate,
360  const Context& /*context*/,
361  unsigned /*spaceIdx*/,
362  unsigned /*timeIdx*/) const
363  { rate = Scalar(0.0); }
364 
366 
367 private:
368  bool onLeftSide_(const GlobalPosition& pos) const
369  { return pos[0] < (this->boundingBoxMin()[0] + this->boundingBoxMax()[0]) / 2; }
370 
371  void setupInitialFluidStates_()
372  {
373  // create the initial fluid state for the left half of the domain
374  leftInitialFluidState_.setTemperature(temperature_);
375 
376  Scalar Sl = 0.0;
377  leftInitialFluidState_.setSaturation(liquidPhaseIdx, Sl);
378  leftInitialFluidState_.setSaturation(gasPhaseIdx, 1 - Sl);
379 
380  Scalar p = 1e5;
381  leftInitialFluidState_.setPressure(liquidPhaseIdx, p);
382  leftInitialFluidState_.setPressure(gasPhaseIdx, p);
383 
384  Scalar xH2O = 0.01;
385  leftInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
386  leftInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
387 
388  using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
389  typename FluidSystem::template ParameterCache<Scalar> paramCache;
390  CFRP::solve(leftInitialFluidState_, paramCache, gasPhaseIdx,
391  /*setViscosity=*/false, /*setEnthalpy=*/false);
392 
393  // create the initial fluid state for the right half of the domain
394  rightInitialFluidState_.assign(leftInitialFluidState_);
395  xH2O = 0.0;
396  rightInitialFluidState_.setMoleFraction(gasPhaseIdx, H2OIdx, xH2O);
397  rightInitialFluidState_.setMoleFraction(gasPhaseIdx, N2Idx, 1 - xH2O);
398  CFRP::solve(rightInitialFluidState_, paramCache, gasPhaseIdx,
399  /*setViscosity=*/false, /*setEnthalpy=*/false);
400  }
401 
402  DimMatrix K_;
403  MaterialLawParams materialParams_;
404 
405  Opm::CompositionalFluidState<Scalar, FluidSystem> leftInitialFluidState_;
406  Opm::CompositionalFluidState<Scalar, FluidSystem> rightInitialFluidState_;
407  Scalar temperature_;
408 };
409 
410 } // namespace Opm
411 
412 #endif
1D problem which is driven by molecular diffusion.
Definition: diffusionproblem.hh:171
const MaterialLawParams & materialLawParams(const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:296
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:305
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: diffusionproblem.hh:340
void finishInit()
Definition: diffusionproblem.hh:221
Scalar porosity(const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:286
void boundary(BoundaryRateVector &values, const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:323
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:359
DiffusionProblem(Simulator &simulator)
Definition: diffusionproblem.hh:214
void endTimeStep()
Definition: diffusionproblem.hh:250
std::string name() const
Definition: diffusionproblem.hh:244
const DimMatrix & intrinsicPermeability(const Context &, unsigned, unsigned) const
Definition: diffusionproblem.hh:277
Definition: diffusionproblem.hh:58