00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
00023 #define OPM_MSWELLHELPERS_HEADER_INCLUDED
00024
00025 #include <opm/common/ErrorMacros.hpp>
00026 #include <dune/istl/solvers.hh>
00027 #if HAVE_UMFPACK
00028 #include <dune/istl/umfpack.hh>
00029 #endif // HAVE_UMFPACK
00030 #include <cmath>
00031
00032 namespace Opm {
00033
00034 namespace mswellhelpers
00035 {
00036 #if HAVE_UMFPACK
00037
00038 template <typename MatrixType, typename VectorType>
00039 VectorType
00040 invDXDirect(const MatrixType& D, VectorType x)
00041 {
00042 VectorType y(x.size());
00043 y = 0.;
00044
00045 Dune::UMFPack<MatrixType> linsolver(D, 0);
00046
00047
00048 Dune::InverseOperatorResult res;
00049
00050
00051 linsolver.apply(y, x, res);
00052
00053
00054
00055 for (size_t i_block = 0; i_block < y.size(); ++i_block) {
00056 for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) {
00057 if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) {
00058 OPM_THROW(Opm::NumericalProblem, "nan or inf value found in invDXDirect due to singular matrix");
00059 }
00060 }
00061 }
00062
00063 return y;
00064 }
00065 #endif // HAVE_UMFPACK
00066
00067
00068
00069
00070
00071
00072 template <typename MatrixType, typename VectorType>
00073 VectorType
00074 invDX(const MatrixType& D, VectorType x)
00075 {
00076
00077
00078
00079
00080
00081
00082 VectorType y(x.size());
00083 y = 0.;
00084
00085 Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
00086
00087
00088 Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
00089
00090
00091
00092
00093
00094 Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
00095 preconditioner,
00096 1.e-8,
00097 250,
00098 0);
00099
00100
00101 Dune::InverseOperatorResult res;
00102
00103
00104 linsolver.apply(y, x, res);
00105
00106 if ( !res.converged ) {
00107 OPM_THROW(Opm::NumericalProblem, "the invDX does not get converged! ");
00108 }
00109
00110 return y;
00111 }
00112
00113
00114
00115
00116
00117 inline double haalandFormular(const double re, const double diameter, const double roughness)
00118 {
00119 const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
00120
00121
00122 assert(value >= 0.0);
00123
00124 return 1. / (value * value);
00125 }
00126
00127
00128
00129
00130
00131 inline double calculateFrictionFactor(const double area, const double diameter,
00132 const double w, const double roughness, const double mu)
00133 {
00134
00135 double f = 0.;
00136
00137 const double re = std::abs(diameter * w / (area * mu));
00138
00139 if ( re == 0.0 ) {
00140
00141 assert(w == 0.);
00142 return 0.0;
00143 }
00144
00145 const double re_value1 = 200.;
00146 const double re_value2 = 4000.;
00147
00148 if (re < re_value1) {
00149 f = 16. / re;
00150 } else if (re > re_value2){
00151 f = haalandFormular(re, diameter, roughness);
00152 } else {
00153 const double f1 = 16. / re_value1;
00154 const double f2 = haalandFormular(re_value2, diameter, roughness);
00155
00156 f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
00157 }
00158 return f;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 template <typename ValueType>
00175 ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness,
00176 const ValueType& density, const ValueType& w, const ValueType& mu)
00177 {
00178 const double f = calculateFrictionFactor(area, diameter, w.value(), roughness, mu.value());
00179
00180 return 2. * f * l * w * w / (area * area * diameter * density);
00181 }
00182
00183
00184
00185
00186
00187 template <typename ValueType>
00188 ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density)
00189 {
00190 return (0.5 * mass_rate * mass_rate / (area * area * density));
00191 }
00192
00193
00194 }
00195
00196 }
00197
00198 #endif