My Project
MSWellHelpers.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4  Copyright 2020 Equinor ASA.
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 
23 #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
24 #define OPM_MSWELLHELPERS_HEADER_INCLUDED
25 
26 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
27 #include <opm/simulators/utils/DeferredLogger.hpp>
28 #include <opm/common/ErrorMacros.hpp>
29 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/SICD.hpp>
30 #include <opm/common/OpmLog/OpmLog.hpp>
31 #include <string>
32 #include<dune/istl/matrix.hh>
33 #include <dune/istl/preconditioners.hh>
34 #include <dune/istl/solvers.hh>
35 
36 #if HAVE_UMFPACK
37 #include <dune/istl/umfpack.hh>
38 #endif // HAVE_UMFPACK
39 #include <cmath>
40 
41 namespace Opm {
42 
43 namespace mswellhelpers
44 {
45 
47  template <typename MatrixType, typename VectorType>
48  VectorType
49  applyUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver, VectorType x)
50  {
51 #if HAVE_UMFPACK
52  if (!linsolver)
53  {
54  linsolver.reset(new Dune::UMFPack<MatrixType>(D, 0));
55  }
56 
57  // The copy of x seems mandatory for calling UMFPack!
58  VectorType y(x.size());
59  y = 0.;
60 
61  // Object storing some statistics about the solving process
62  Dune::InverseOperatorResult res;
63 
64  // Solve
65  linsolver->apply(y, x, res);
66 
67  // Checking if there is any inf or nan in y
68  // it will be the solution before we find a way to catch the singularity of the matrix
69  for (size_t i_block = 0; i_block < y.size(); ++i_block) {
70  for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) {
71  if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) {
72  const std::string msg{"nan or inf value found after UMFPack solve due to singular matrix"};
73  OpmLog::debug(msg);
74  OPM_THROW_NOLOG(NumericalIssue, msg);
75  }
76  }
77  }
78  return y;
79 #else
80  // this is not thread safe
81  OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. "
82  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
83 #endif // HAVE_UMFPACK
84  }
85 
86 
87 
89  template <typename MatrixType, typename VectorType>
90  Dune::Matrix<typename MatrixType::block_type>
91  invertWithUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver)
92  {
93 #if HAVE_UMFPACK
94  const int sz = D.M();
95  const int bsz = D[0][0].M();
96  VectorType e(sz);
97  e = 0.0;
98 
99  // Make a full block matrix.
100  Dune::Matrix<typename MatrixType::block_type> inv(sz, sz);
101 
102  // Create inverse by passing basis vectors to the solver.
103  for (int ii = 0; ii < sz; ++ii) {
104  for (int jj = 0; jj < bsz; ++jj) {
105  e[ii][jj] = 1.0;
106  auto col = applyUMFPack(D, linsolver, e);
107  for (int cc = 0; cc < sz; ++cc) {
108  for (int dd = 0; dd < bsz; ++dd) {
109  inv[cc][ii][dd][jj] = col[cc][dd];
110  }
111  }
112  e[ii][jj] = 0.0;
113  }
114  }
115 
116  return inv;
117 #else
118  // this is not thread safe
119  OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. "
120  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
121 #endif // HAVE_UMFPACK
122  }
123 
124 
125 
126  // obtain y = D^-1 * x with a BICSSTAB iterative solver
127  template <typename MatrixType, typename VectorType>
128  VectorType
129  invDX(const MatrixType& D, VectorType x, DeferredLogger& deferred_logger)
130  {
131  // the function will change the value of x, so we should not use reference of x here.
132 
133  // TODO: store some of the following information to avoid to call it again and again for
134  // efficiency improvement.
135  // Bassically, only the solve / apply step is different.
136 
137  VectorType y(x.size());
138  y = 0.;
139 
140  Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
141 
142  // Sequential incomplete LU decomposition as the preconditioner
143 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
144  Dune::SeqILU<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
145 #else
146  Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
147 #endif
148  // Dune::SeqILUn<MatrixType, VectorType, VectorType> preconditioner(D, 1, 0.92);
149  // Dune::SeqGS<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
150  // Dune::SeqJac<MatrixType, VectorType, VectorType> preconditioner(D, 1, 1);
151 
152  // Preconditioned BICGSTAB solver
153  Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
154  preconditioner,
155  1.e-8, // desired residual reduction factor
156  250, // maximum number of iterations
157  0); // verbosity of the solver */
158 
159  // Object storing some statistics about the solving process
160  Dune::InverseOperatorResult res;
161 
162  // Solve
163  linsolver.apply(y, x, res);
164 
165  if ( !res.converged ) {
166  OPM_DEFLOG_THROW(NumericalIssue, "the invDX did not converge ", deferred_logger);
167  }
168 
169  return y;
170  }
171 
172 
173 
174 
175  template <typename ValueType>
176  inline ValueType haalandFormular(const ValueType& re, const double diameter, const double roughness)
177  {
178  const ValueType value = -3.6 * log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
179 
180  // sqrt(1/f) should be non-positive
181  assert(value >= 0.0);
182 
183  return 1. / (value * value);
184  }
185 
186 
187 
188 
189  template <typename ValueType>
190  inline ValueType calculateFrictionFactor(const double area, const double diameter,
191  const ValueType& w, const double roughness, const ValueType& mu)
192  {
193 
194  ValueType f = 0.;
195  // Reynolds number
196  const ValueType re = abs( diameter * w / (area * mu));
197 
198  if ( re == 0.0 ) {
199  // make sure it is because the mass rate is zero
200  assert(w == 0.);
201  return 0.0;
202  }
203 
204  const ValueType re_value1 = 2000.;
205  const ValueType re_value2 = 4000.;
206 
207  if (re < re_value1) {
208  f = 16. / re;
209  } else if (re > re_value2){
210  f = haalandFormular(re, diameter, roughness);
211  } else { // in between
212  const ValueType f1 = 16. / re_value1;
213  const ValueType f2 = haalandFormular(re_value2, diameter, roughness);
214 
215  f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
216  }
217  return f;
218  }
219 
220 
221 
222 
223 
224 
225  // calculating the friction pressure loss
226  // l is the segment length
227  // area is the segment cross area
228  // diameter is the segment inner diameter
229  // w is mass flow rate through the segment
230  // density is density
231  // roughness is the absolute roughness
232  // mu is the average phase viscosity
233  template <typename ValueType>
234  ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness,
235  const ValueType& density, const ValueType& w, const ValueType& mu)
236  {
237  const ValueType f = calculateFrictionFactor(area, diameter, w, roughness, mu);
238  // \Note: a factor of 2 needs to be here based on the dimensional analysis
239  return 2. * f * l * w * w / (area * area * diameter * density);
240  }
241 
242 
243  template <typename ValueType>
244  ValueType valveContrictionPressureLoss(const ValueType& mass_rate, const ValueType& density,
245  const double area_con, const double cv)
246  {
247  // the formulation is adjusted a little bit for convinience
248  // velocity = mass_rate / (density * area) is applied to the original formulation
249  const double area = (area_con > 1.e-10 ? area_con : 1.e-10);
250  return mass_rate * mass_rate / (2. * density * cv * cv * area * area);
251  }
252 
253 
254  template <typename ValueType>
255  ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density)
256  {
257  // \Note: a factor of 2 is added to the formulation in order to match results from the
258  // reference simulator. This is inline with what is done for the friction loss.
259  return (mass_rate * mass_rate / (area * area * density));
260  }
261 
262 
263 
264  // water in oil emulsion viscosity
265  // TODO: maybe it should be two different ValueTypes. When we calculate the viscosity for transitional zone
266  template <typename ValueType>
267  ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const ValueType& water_liquid_fraction,
268  const double max_visco_ratio)
269  {
270  const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) );
271  const ValueType viscosity_ratio = pow(temp_value, 2.5);
272 
273  if (viscosity_ratio <= max_visco_ratio) {
274  return oil_viscosity * viscosity_ratio;
275  } else {
276  return oil_viscosity * max_visco_ratio;
277  }
278  }
279 
280 
281 
282 
283 
284  // oil in water emulsion viscosity
285  template <typename ValueType>
286  ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const ValueType& water_liquid_fraction,
287  const double max_visco_ratio)
288  {
289  const ValueType temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) );
290  const ValueType viscosity_ratio = pow(temp_value, 2.5);
291 
292  if (viscosity_ratio <= max_visco_ratio) {
293  return water_viscosity * viscosity_ratio;
294  } else {
295  return water_viscosity * max_visco_ratio;
296  }
297  }
298 
299 
300 
301 
302 
303  // calculating the viscosity of oil-water emulsion at local conditons
304  template <typename ValueType>
305  ValueType emulsionViscosity(const ValueType& water_fraction, const ValueType& water_viscosity,
306  const ValueType& oil_fraction, const ValueType& oil_viscosity,
307  const SICD& sicd)
308  {
309  const double width_transition = sicd.widthTransitionRegion();
310 
311  // it is just for now, we should be able to treat it.
312  if (width_transition <= 0.) {
313  OPM_THROW(std::runtime_error, "Not handling non-positive transition width now");
314  }
315 
316  const double critical_value = sicd.criticalValue();
317  const ValueType transition_start_value = critical_value - width_transition / 2.0;
318  const ValueType transition_end_value = critical_value + width_transition / 2.0;
319 
320  const ValueType liquid_fraction = water_fraction + oil_fraction;
321  // if there is no liquid, we just return zero
322  if (liquid_fraction == 0.) {
323  return 0.;
324  }
325 
326  const ValueType water_liquid_fraction = water_fraction / liquid_fraction;
327 
328  const double max_visco_ratio = sicd.maxViscosityRatio();
329  if (water_liquid_fraction <= transition_start_value) {
330  return WIOEmulsionViscosity(oil_viscosity, water_liquid_fraction, max_visco_ratio);
331  } else if(water_liquid_fraction >= transition_end_value) {
332  return OIWEmulsionViscosity(water_viscosity, water_liquid_fraction, max_visco_ratio);
333  } else { // in the transition region
334  const ValueType viscosity_start_transition = WIOEmulsionViscosity(oil_viscosity, transition_start_value, max_visco_ratio);
335  const ValueType viscosity_end_transition = OIWEmulsionViscosity(water_viscosity, transition_end_value, max_visco_ratio);
336  const ValueType emulsion_viscosity = (viscosity_start_transition * (transition_end_value - water_liquid_fraction)
337  + viscosity_end_transition * (water_liquid_fraction - transition_start_value) ) / width_transition;
338  return emulsion_viscosity;
339  }
340  }
341 
342 } // namespace mswellhelpers
343 
344 }
345 
346 #endif
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26