23 #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED
24 #define OPM_MSWELLHELPERS_HEADER_INCLUDED
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>
32 #include<dune/istl/matrix.hh>
33 #include <dune/istl/preconditioners.hh>
34 #include <dune/istl/solvers.hh>
37 #include <dune/istl/umfpack.hh>
43 namespace mswellhelpers
47 template <
typename MatrixType,
typename VectorType>
49 applyUMFPack(
const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver, VectorType x)
54 linsolver.reset(
new Dune::UMFPack<MatrixType>(D, 0));
58 VectorType y(x.size());
62 Dune::InverseOperatorResult res;
65 linsolver->apply(y, x, res);
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"};
74 OPM_THROW_NOLOG(NumericalIssue, msg);
81 OPM_THROW(std::runtime_error,
"Cannot use applyUMFPack() without UMFPACK. "
82 "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
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)
95 const int bsz = D[0][0].M();
100 Dune::Matrix<typename MatrixType::block_type> inv(sz, sz);
103 for (
int ii = 0; ii < sz; ++ii) {
104 for (
int jj = 0; jj < bsz; ++jj) {
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];
119 OPM_THROW(std::runtime_error,
"Cannot use invertWithUMFPack() without UMFPACK. "
120 "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
127 template <
typename MatrixType,
typename VectorType>
129 invDX(
const MatrixType& D, VectorType x, DeferredLogger& deferred_logger)
137 VectorType y(x.size());
140 Dune::MatrixAdapter<MatrixType, VectorType, VectorType> linearOperator(D);
143 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
144 Dune::SeqILU<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
146 Dune::SeqILU0<MatrixType, VectorType, VectorType> preconditioner(D, 1.0);
153 Dune::BiCGSTABSolver<VectorType> linsolver(linearOperator,
160 Dune::InverseOperatorResult res;
163 linsolver.apply(y, x, res);
165 if ( !res.converged ) {
166 OPM_DEFLOG_THROW(NumericalIssue,
"the invDX did not converge ", deferred_logger);
175 template <
typename ValueType>
176 inline ValueType haalandFormular(
const ValueType& re,
const double diameter,
const double roughness)
178 const ValueType value = -3.6 * log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
181 assert(value >= 0.0);
183 return 1. / (value * value);
189 template <
typename ValueType>
190 inline ValueType calculateFrictionFactor(
const double area,
const double diameter,
191 const ValueType& w,
const double roughness,
const ValueType& mu)
196 const ValueType re = abs( diameter * w / (area * mu));
204 const ValueType re_value1 = 2000.;
205 const ValueType re_value2 = 4000.;
207 if (re < re_value1) {
209 }
else if (re > re_value2){
210 f = haalandFormular(re, diameter, roughness);
212 const ValueType f1 = 16. / re_value1;
213 const ValueType f2 = haalandFormular(re_value2, diameter, roughness);
215 f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1;
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)
237 const ValueType f = calculateFrictionFactor(area, diameter, w, roughness, mu);
239 return 2. * f * l * w * w / (area * area * diameter * density);
243 template <
typename ValueType>
244 ValueType valveContrictionPressureLoss(
const ValueType& mass_rate,
const ValueType& density,
245 const double area_con,
const double cv)
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);
254 template <
typename ValueType>
255 ValueType velocityHead(
const double area,
const ValueType& mass_rate,
const ValueType& density)
259 return (mass_rate * mass_rate / (area * area * density));
266 template <
typename ValueType>
267 ValueType WIOEmulsionViscosity(
const ValueType& oil_viscosity,
const ValueType& water_liquid_fraction,
268 const double max_visco_ratio)
270 const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) );
271 const ValueType viscosity_ratio = pow(temp_value, 2.5);
273 if (viscosity_ratio <= max_visco_ratio) {
274 return oil_viscosity * viscosity_ratio;
276 return oil_viscosity * max_visco_ratio;
285 template <
typename ValueType>
286 ValueType OIWEmulsionViscosity(
const ValueType& water_viscosity,
const ValueType& water_liquid_fraction,
287 const double max_visco_ratio)
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);
292 if (viscosity_ratio <= max_visco_ratio) {
293 return water_viscosity * viscosity_ratio;
295 return water_viscosity * max_visco_ratio;
304 template <
typename ValueType>
305 ValueType emulsionViscosity(
const ValueType& water_fraction,
const ValueType& water_viscosity,
306 const ValueType& oil_fraction,
const ValueType& oil_viscosity,
309 const double width_transition = sicd.widthTransitionRegion();
312 if (width_transition <= 0.) {
313 OPM_THROW(std::runtime_error,
"Not handling non-positive transition width now");
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;
320 const ValueType liquid_fraction = water_fraction + oil_fraction;
322 if (liquid_fraction == 0.) {
326 const ValueType water_liquid_fraction = water_fraction / liquid_fraction;
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);
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;
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26