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