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 <dune/grid/common/GridAdapter.hpp>
18 
19 #include <opm/core/grid.h>
20 #include <opm/core/linalg/LinearSolverIstl.hpp>
21 #include <opm/core/linalg/sparse_sys.h>
22 #include <opm/core/pressure/tpfa/ifs_tpfa.h>
23 #include <opm/core/pressure/tpfa/trans_tpfa.h>
24 #include <opm/core/transport/implicit/NormSupport.hpp>
25 #include <opm/core/transport/implicit/ImplicitTransport.hpp>
26 #include <opm/core/utility/parameters/ParameterGroup.hpp>
27 
28 #include <opm/porsol/common/ReservoirPropertyCapillary.hpp>
29 
30 #include <algorithm>
31 #include <cstddef>
32 #include <memory>
33 #include <string>
34 #include <vector>
35 
36 #include <boost/lexical_cast.hpp>
37 
38 namespace Opm {
39  template <class Vector>
40  class MaxNormStl {
41  public:
42  static double
43  norm(const Vector& v) {
44  return ImplicitTransportDefault::AccumulationNorm <Vector, ImplicitTransportDefault::MaxAbs>::norm(v);
45  }
46  };
47 
48  template <class Vector>
49  class MaxNormDune {
50  public:
51  static double
52  norm(const Vector& v) {
53  return v.infinity_norm();
54  }
55  };
56 
57  template <int np = 2>
59  public:
60  ReservoirState(const UnstructuredGrid* g)
61  : press_ (g->number_of_cells),
62  fpress_(g->number_of_faces),
63  flux_ (g->number_of_faces),
64  sat_ (np * g->number_of_cells)
65  {}
66 
67  ::std::vector<double>& pressure () { return press_ ; }
68  ::std::vector<double>& facepressure() { return fpress_; }
69  ::std::vector<double>& faceflux () { return flux_ ; }
70  ::std::vector<double>& saturation () { return sat_ ; }
71 
72  const ::std::vector<double>& faceflux () const { return flux_; }
73  const ::std::vector<double>& saturation () const { return sat_ ; }
74 
75  private:
76  ::std::vector<double> press_ ;
77  ::std::vector<double> fpress_;
78  ::std::vector<double> flux_ ;
79  ::std::vector<double> sat_ ;
80  };
81 
82  class Rock {
83  public:
84  Rock(::std::size_t nc, ::std::size_t dim)
85  : dim_ (dim ),
86  perm_(nc * dim * dim),
87  poro_(nc ) {}
88 
89  const ::std::vector<double>& perm() const { return perm_; }
90  const ::std::vector<double>& poro() const { return poro_; }
91 
92  void
93  perm_homogeneous(double k) {
94  setVector(0.0, perm_);
95 
96  const ::std::size_t d2 = dim_ * dim_;
97 
98  for (::std::size_t c = 0, nc = poro_.size(); c < nc; ++c) {
99  for (::std::size_t i = 0; i < dim_; ++i) {
100  perm_[c*d2 + i*(dim_ + 1)] = k;
101  }
102  }
103  }
104 
105  void
106  poro_homogeneous(double phi) {
107  setVector(phi, poro_);
108  }
109 
110  private:
111  void
112  setVector(double x, ::std::vector<double>& v) {
113  ::std::fill(v.begin(), v.end(), x);
114  }
115 
116  ::std::size_t dim_ ;
117  ::std::vector<double> perm_;
118  ::std::vector<double> poro_;
119  };
120 
122  public:
124  : r_(r)
125  {}
126 
127  template <class Sat ,
128  class Mob ,
129  class DMob>
130  void
131  mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const {
132  const double s1 = s[0];
133 
134  r_.phaseMobilities (c, s1, mob );
135  r_.phaseMobilitiesDeriv(c, s1, dmob);
136  }
137 
138  template <class Sat ,
139  class PC ,
140  class DPC>
141  void
142  pc(int c, const Sat& s, PC& pc, DPC& dpc) const {
143  const double s1 = s[0];
144  pc = r_.capillaryPressure(c, s1);
145  dpc = r_.capillaryPressureDeriv(c, s1);
146  }
147 
148  double density(int p) const {
149  if (p == 0) {
150  return r_.densityFirstPhase();
151  } else {
152  return r_.densitySecondPhase();
153  }
154  }
155 
156  double s_min(int c) const {
157  return r_.s_min(c);
158  }
159 
160  double s_max(int c) const {
161  return r_.s_max(c);
162  }
163 
164  private:
166  };
167 
169  public:
170  typedef Dune::FieldVector<double, 1> ScalarVectorBlockType;
171  typedef Dune::FieldMatrix<double, 1, 1> ScalarMatrixBlockType;
172 
173  typedef Dune::BlockVector<ScalarVectorBlockType> ScalarBlockVector;
174  typedef Dune::BCRSMatrix <ScalarMatrixBlockType> ScalarBCRSMatrix;
175 
176  void
177  solve(const ScalarBCRSMatrix& A,
178  const ScalarBlockVector& b,
179  ScalarBlockVector& x)
180  {
181  Dune::MatrixAdapter<ScalarBCRSMatrix,
182  ScalarBlockVector,
183  ScalarBlockVector> opA(A);
184 
185  Dune::SeqILU0<ScalarBCRSMatrix,
186  ScalarBlockVector,
187  ScalarBlockVector> precond(A, 1.0);
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:
205  {
206  Opm::ParameterGroup params;
207 
208  params.insertParameter("linsolver_tolerance",
209  boost::lexical_cast<std::string>(5.0e-9));
210 
211  params.insertParameter("linsolver_verbosity",
212  boost::lexical_cast<std::string>(1));
213 
214  params.insertParameter("linsolver_type",
215  boost::lexical_cast<std::string>(1));
216 
217  ls_.reset(new LinearSolverIstl(params));
218  }
219 
220  void
221  solve(struct CSRMatrix* A,
222  const double* b,
223  double* x)
224  {
225  ls_->solve(A->m, A->nnz, A->ia, A->ja, A->sa, b, x);
226  }
227 
228  private:
229  std::unique_ptr<Opm::LinearSolverIstl> ls_;
230  };
231 
233  public:
234  TransportSource() : nsrc(0), pf(0) {}
235 
236  int nsrc;
237  int pf;
238  ::std::vector< int > cell;
239  ::std::vector<double> pressure;
240  ::std::vector<double> flux;
241  ::std::vector<double> saturation;
242  ::std::vector<double> periodic_cells;
243  ::std::vector<double> periodic_faces;
244  };
245 
246  template <class Arr>
247  void
248  append_transport_source(int c, double p, double v, const Arr& s,
249  TransportSource& src)
250  {
251  src.cell .push_back(c);
252  src.pressure .push_back(p);
253  src.flux .push_back(v);
254  src.saturation.insert(src.saturation.end(),
255  s.begin(), s.end());
256  ++src.nsrc;
257  }
258 
259  void
260  compute_porevolume(const UnstructuredGrid* g,
261  const Rock& rock,
262  std::vector<double>& porevol);
263 } // namespace Opm
264 
265 #endif // OPENRS_IMPLICITTRANSPORTDEFS_HEADER
Definition: ImplicitTransportDefs.hpp:202
Definition: ImplicitTransportDefs.hpp:168
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 densitySecondPhase() const
Density of second (oil) phase.
Definition: ReservoirPropertyCommon_impl.hpp:373
Definition: ImplicitTransportDefs.hpp:40
Definition: ImplicitTransportDefs.hpp:121
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
Rock()
Default constructor.
Definition: Rock_impl.hpp:37
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:82
double s_max(int c) const
Definition: ReservoirPropertyCommon_impl.hpp:506
Definition: ImplicitTransportDefs.hpp:58
Definition: ImplicitTransportDefs.hpp:232
Definition: ImplicitTransportDefs.hpp:49
double capillaryPressure(int cell_index, double saturation) const
Capillary pressure.
Definition: ReservoirPropertyCommon_impl.hpp:467
void phaseMobilities(int cell_index, double saturation, Vector &mobility) const
Mobilities for both phases.
Definition: ReservoirPropertyCapillary_impl.hpp:93
double densityFirstPhase() const
Density of first (water) phase.
Definition: ReservoirPropertyCommon_impl.hpp:366