22 #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
23 #define OPM_MSWELLHELPERS_HEADER_INCLUDED
25 #include <opm/common/ErrorMacros.hpp>
26 #include <dune/istl/solvers.hh>
28 #include <dune/istl/umfpack.hh>
29 #endif // HAVE_UMFPACK
34 namespace mswellhelpers
38 template <
typename MatrixType,
typename VectorType>
40 invDXDirect(
const MatrixType& D, VectorType x)
42 VectorType y(x.size());
45 Dune::UMFPack<MatrixType> linsolver(D, 0);
48 Dune::InverseOperatorResult res;
51 linsolver.apply(y, x, res);
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");
65 #endif // HAVE_UMFPACK
72 template <
typename MatrixType,
typename VectorType>
74 invDX(
const MatrixType& D, VectorType x)
82 VectorType y(x.size());
85 Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
88 Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
94 Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
101 Dune::InverseOperatorResult res;
104 linsolver.apply(y, x, res);
106 if ( !res.converged ) {
107 OPM_THROW(Opm::NumericalProblem,
"the invDX does not get converged! ");
117 inline double haalandFormular(
const double re,
const double diameter,
const double roughness)
119 const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
122 assert(value >= 0.0);
124 return 1. / (value * value);
131 inline double calculateFrictionFactor(
const double area,
const double diameter,
132 const double w,
const double roughness,
const double mu)
137 const double re = std::abs(diameter * w / (area * mu));
145 const double re_value1 = 200.;
146 const double re_value2 = 4000.;
148 if (re < re_value1) {
150 }
else if (re > re_value2){
151 f = haalandFormular(re, diameter, roughness);
153 const double f1 = 16. / re_value1;
154 const double f2 = haalandFormular(re_value2, diameter, roughness);
156 f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
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)
178 const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value());
180 return 2. * f * l * w * w / (area * area * diameter * density);
187 template <
typename ValueType>
188 ValueType velocityHead(
const double area,
const ValueType& mass_rate,
const ValueType& density)
190 return (0.5 * mass_rate * mass_rate / (area * area * density));