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
32namespace 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>
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