All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
MSWellHelpers.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
23 #define OPM_MSWELLHELPERS_HEADER_INCLUDED
24 
25 #include <opm/common/ErrorMacros.hpp>
26 #include <dune/istl/solvers.hh>
27 #if HAVE_UMFPACK
28 #include <dune/istl/umfpack.hh>
29 #endif // HAVE_UMFPACK
30 #include <cmath>
31 
32 namespace Opm {
33 
34 namespace mswellhelpers
35 {
36 #if HAVE_UMFPACK
37  // obtain y = D^-1 * x with a direct solver
38  template <typename MatrixType, typename VectorType>
39  VectorType
40  invDXDirect(const MatrixType& D, VectorType x)
41  {
42  VectorType y(x.size());
43  y = 0.;
44 
45  Dune::UMFPack<MatrixType> linsolver(D, 0);
46 
47  // Object storing some statistics about the solving process
48  Dune::InverseOperatorResult res;
49 
50  // Solve
51  linsolver.apply(y, x, res);
52 
53  // Checking if there is any inf or nan in y
54  // it will be the solution before we find a way to catch the singularity of the matrix
55  for (size_t i_block = 0; i_block < y.size(); ++i_block) {
56  for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) {
57  if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) {
58  OPM_THROW(Opm::NumericalProblem, "nan or inf value found in invDXDirect due to singular matrix");
59  }
60  }
61  }
62 
63  return y;
64  }
65 #endif // HAVE_UMFPACK
66 
67 
68 
69 
70 
71  // obtain y = D^-1 * x with a BICSSTAB iterative solver
72  template <typename MatrixType, typename VectorType>
73  VectorType
74  invDX(const MatrixType& D, VectorType x)
75  {
76  // the function will change the value of x, so we should not use reference of x here.
77 
78  // TODO: store some of the following information to avoid to call it again and again for
79  // efficiency improvement.
80  // Bassically, only the solve / apply step is different.
81 
82  VectorType y(x.size());
83  y = 0.;
84 
85  Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
86 
87  // Sequential incomplete LU decomposition as the preconditioner
88  Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
89  // Dune::SeqILUn<MatrixType, VectorType, VectorType> preconditioner(D, 1, 0.92);
90  // Dune::SeqGS<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
91  // Dune::SeqJac<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
92 
93  // Preconditioned BICGSTAB solver
94  Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
95  preconditioner,
96  1.e-8, // desired residual reduction factor
97  250, // maximum number of iterations
98  0); // verbosity of the solver */
99 
100  // Object storing some statistics about the solving process
101  Dune::InverseOperatorResult res;
102 
103  // Solve
104  linsolver.apply(y, x, res);
105 
106  if ( !res.converged ) {
107  OPM_THROW(Opm::NumericalProblem, "the invDX does not get converged! ");
108  }
109 
110  return y;
111  }
112 
113 
114 
115 
116 
117  inline double haalandFormular(const double re, const double diameter, const double roughness)
118  {
119  const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
120 
121  // sqrt(1/f) should be non-positive
122  assert(value >= 0.0);
123 
124  return 1. / (value * value);
125  }
126 
127 
128 
129 
130 
131  inline double calculateFrictionFactor(const double area, const double diameter,
132  const double w, const double roughness, const double mu)
133  {
134 
135  double f = 0.;
136  // Reynolds number
137  const double re = std::abs(diameter * w / (area * mu));
138 
139  if ( re == 0.0 ) {
140  // make sure it is because the mass rate is zero
141  assert(w == 0.);
142  return 0.0;
143  }
144 
145  const double re_value1 = 200.;
146  const double re_value2 = 4000.;
147 
148  if (re < re_value1) {
149  f = 16. / re;
150  } else if (re > re_value2){
151  f = haalandFormular(re, diameter, roughness);
152  } else { // in between
153  const double f1 = 16. / re_value1;
154  const double f2 = haalandFormular(re_value2, diameter, roughness);
155 
156  f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
157  }
158  return f;
159  }
160 
161 
162 
163 
164 
165 
166  // calculating the friction pressure loss
167  // l is the segment length
168  // area is the segment cross area
169  // diameter is the segment inner diameter
170  // w is mass flow rate through the segment
171  // density is density
172  // roughness is the absolute roughness
173  // mu is the average phase viscosity
174  template <typename ValueType>
175  ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness,
176  const ValueType& density, const ValueType& w, const ValueType& mu)
177  {
178  const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value());
179  // \Note: a factor of 2 needs to be here based on the dimensional analysis
180  return 2. * f * l * w * w / (area * area * diameter * density);
181  }
182 
183 
184 
185 
186 
187  template <typename ValueType>
188  ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density)
189  {
190  return (0.5 * mass_rate * mass_rate / (area * area * density));
191  }
192 
193 
194 } // namespace mswellhelpers
195 
196 }
197 
198 #endif