My Project
ImplicitTransportDefs.hpp
1 #ifndef OPENRS_IMPLICITTRANSPORTDEFS_HEADER
2 #define OPENRS_IMPLICITTRANSPORTDEFS_HEADER
3 
4 #include <opm/common/utility/platform_dependent/disable_warnings.h>
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 
9 #include <dune/istl/operators.hh>
10 #include <dune/istl/solvers.hh>
11 #include <dune/istl/bvector.hh>
12 #include <dune/istl/bcrsmatrix.hh>
13 #include <dune/istl/preconditioners.hh>
14 
15 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
16 
17 #include <opm/grid/common/GridAdapter.hpp>
18 
19 #include <opm/grid/UnstructuredGrid.h>
20 #include <opm/core/transport/implicit/NormSupport.hpp>
21 #include <opm/core/transport/implicit/ImplicitTransport.hpp>
22 #include <opm/common/utility/parameters/ParameterGroup.hpp>
23 
24 #include <opm/porsol/common/ReservoirPropertyCapillary.hpp>
25 
26 #include <algorithm>
27 #include <cstddef>
28 #include <memory>
29 #include <string>
30 #include <vector>
31 
32 namespace Opm {
33  template <class Vector>
34  class MaxNormStl {
35  public:
36  static double
37  norm(const Vector& v) {
39  }
40  };
41 
42  template <class Vector>
43  class MaxNormDune {
44  public:
45  static double
46  norm(const Vector& v) {
47  return v.infinity_norm();
48  }
49  };
50 
51  template <int np = 2>
53  public:
54  ReservoirState(const UnstructuredGrid* g)
55  : press_ (g->number_of_cells),
56  fpress_(g->number_of_faces),
57  flux_ (g->number_of_faces),
58  sat_ (np * g->number_of_cells)
59  {}
60 
61  ::std::vector<double>& pressure () { return press_ ; }
62  ::std::vector<double>& facepressure() { return fpress_; }
63  ::std::vector<double>& faceflux () { return flux_ ; }
64  ::std::vector<double>& saturation () { return sat_ ; }
65 
66  const ::std::vector<double>& faceflux () const { return flux_; }
67  const ::std::vector<double>& saturation () const { return sat_ ; }
68 
69  private:
70  ::std::vector<double> press_ ;
71  ::std::vector<double> fpress_;
72  ::std::vector<double> flux_ ;
73  ::std::vector<double> sat_ ;
74  };
75 
76  class Rock {
77  public:
78  Rock(::std::size_t nc, ::std::size_t dim)
79  : dim_ (dim ),
80  perm_(nc * dim * dim),
81  poro_(nc ) {}
82 
83  const ::std::vector<double>& perm() const { return perm_; }
84  const ::std::vector<double>& poro() const { return poro_; }
85 
86  void
87  perm_homogeneous(double k) {
88  setVector(0.0, perm_);
89 
90  const ::std::size_t d2 = dim_ * dim_;
91 
92  for (::std::size_t c = 0, nc = poro_.size(); c < nc; ++c) {
93  for (::std::size_t i = 0; i < dim_; ++i) {
94  perm_[c*d2 + i*(dim_ + 1)] = k;
95  }
96  }
97  }
98 
99  void
100  poro_homogeneous(double phi) {
101  setVector(phi, poro_);
102  }
103 
104  private:
105  void
106  setVector(double x, ::std::vector<double>& v) {
107  ::std::fill(v.begin(), v.end(), x);
108  }
109 
110  ::std::size_t dim_ ;
111  ::std::vector<double> perm_;
112  ::std::vector<double> poro_;
113  };
114 
116  public:
118  : r_(r)
119  {}
120 
121  template <class Sat ,
122  class Mob ,
123  class DMob>
124  void
125  mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const {
126  const double s1 = s[0];
127 
128  r_.phaseMobilities (c, s1, mob );
129  r_.phaseMobilitiesDeriv(c, s1, dmob);
130  }
131 
132  template <class Sat ,
133  class PC ,
134  class DPC>
135  void
136  pc(int c, const Sat& s, PC& pc_arg, DPC& dpc) const {
137  const double s1 = s[0];
138  pc_arg = r_.capillaryPressure(c, s1);
139  dpc = r_.capillaryPressureDeriv(c, s1);
140  }
141 
142  double density(int p) const {
143  if (p == 0) {
144  return r_.densityFirstPhase();
145  } else {
146  return r_.densitySecondPhase();
147  }
148  }
149 
150  double s_min(int c) const {
151  return r_.s_min(c);
152  }
153 
154  double s_max(int c) const {
155  return r_.s_max(c);
156  }
157 
158  private:
160  };
161 
163  public:
164  typedef Dune::FieldVector<double, 1> ScalarVectorBlockType;
165  typedef Dune::FieldMatrix<double, 1, 1> ScalarMatrixBlockType;
166 
167  typedef Dune::BlockVector<ScalarVectorBlockType> ScalarBlockVector;
168  typedef Dune::BCRSMatrix <ScalarMatrixBlockType> ScalarBCRSMatrix;
169 
170  void
171  solve(const ScalarBCRSMatrix& A,
172  const ScalarBlockVector& b,
173  ScalarBlockVector& x)
174  {
175  Dune::MatrixAdapter<ScalarBCRSMatrix,
176  ScalarBlockVector,
177  ScalarBlockVector> opA(A);
178 
179 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
180  Dune::SeqILU<ScalarBCRSMatrix,
181  ScalarBlockVector,
182  ScalarBlockVector> precond(A, 1.0);
183 #else
184  Dune::SeqILU0<ScalarBCRSMatrix,
185  ScalarBlockVector,
186  ScalarBlockVector> precond(A, 1.0);
187 #endif
188 
189  int maxit = A.N();
190  double tol = 5.0e-7;
191  int verb = 0;
192 
193  Dune::BiCGSTABSolver<ScalarBlockVector>
194  solver(opA, precond, tol, maxit, verb);
195 
196  ScalarBlockVector bcpy(b);
197  Dune::InverseOperatorResult res;
198  solver.apply(x, bcpy, res);
199  }
200  };
201 
203  public:
204  TransportSource() : nsrc(0), pf(0) {}
205 
206  int nsrc;
207  int pf;
208  ::std::vector< int > cell;
209  ::std::vector<double> pressure;
210  ::std::vector<double> flux;
211  ::std::vector<double> saturation;
212  ::std::vector<double> periodic_cells;
213  ::std::vector<double> periodic_faces;
214  };
215 
216  template <class Arr>
217  void
218  append_transport_source(int c, double p, double v, const Arr& s,
219  TransportSource& src)
220  {
221  src.cell .push_back(c);
222  src.pressure .push_back(p);
223  src.flux .push_back(v);
224  src.saturation.insert(src.saturation.end(),
225  s.begin(), s.end());
226  ++src.nsrc;
227  }
228 
229  void
230  compute_porevolume(const UnstructuredGrid* g,
231  const Rock& rock,
232  std::vector<double>& porevol);
233 } // namespace Opm
234 
235 #endif // OPENRS_IMPLICITTRANSPORTDEFS_HEADER
Definition: ImplicitTransportDefs.hpp:162
Definition: ImplicitTransportDefs.hpp:43
Definition: ImplicitTransportDefs.hpp:34
void phaseMobilities(int cell_index, double saturation, Vector &mobility) const
Mobilities for both phases.
Definition: ReservoirPropertyCapillary_impl.hpp:93
double densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:373
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:366
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:467
double capillaryPressureDeriv(int c, double s) const
Derivative of Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:480
double s_min(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:493
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:506
Definition: ImplicitTransportDefs.hpp:52
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:76
Rock()
Default constructor.
Definition: Rock_impl.hpp:37
Definition: ImplicitTransportDefs.hpp:202
Definition: ImplicitTransportDefs.hpp:115
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43