diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 12d58957d..24a5a435d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -83,6 +83,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/GasLiftStage2.cpp opm/simulators/wells/GlobalWellInfo.cpp opm/simulators/wells/GroupState.cpp + opm/simulators/wells/MSWellHelpers.cpp opm/simulators/wells/MultisegmentWellEval.cpp opm/simulators/wells/MultisegmentWellGeneric.cpp opm/simulators/wells/ParallelWellInfo.cpp diff --git a/opm/simulators/wells/MSWellHelpers.cpp b/opm/simulators/wells/MSWellHelpers.cpp new file mode 100644 index 000000000..89a8941c7 --- /dev/null +++ b/opm/simulators/wells/MSWellHelpers.cpp @@ -0,0 +1,377 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + Copyright 2020 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include + +#include +#include + +#include + +#include + +#include +#include + +#include +#include + +#if HAVE_UMFPACK +#include +#endif // HAVE_UMFPACK + +#include + +namespace { + +template +ValueType haalandFormular(const ValueType& re, + const double diameter, + const double roughness) +{ + const ValueType value = -3.6 * log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) ); + + // sqrt(1/f) should be non-positive + assert(value >= 0.0); + + return 1. / (value * value); +} + +// water in oil emulsion viscosity +// TODO: maybe it should be two different ValueTypes. When we calculate the viscosity for transitional zone +template +ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, + const ValueType& water_liquid_fraction, + const double max_visco_ratio) +{ + const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) ); + const ValueType viscosity_ratio = pow(temp_value, 2.5); + + if (viscosity_ratio <= max_visco_ratio) { + return oil_viscosity * viscosity_ratio; + } else { + return oil_viscosity * max_visco_ratio; + } +} + +// oil in water emulsion viscosity +template +ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, + const ValueType& water_liquid_fraction, + const double max_visco_ratio) +{ + const ValueType temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) ); + const ValueType viscosity_ratio = pow(temp_value, 2.5); + + if (viscosity_ratio <= max_visco_ratio) { + return water_viscosity * viscosity_ratio; + } else { + return water_viscosity * max_visco_ratio; + } +} + +} + +namespace Opm { + +namespace mswellhelpers +{ + + /// Applies umfpack and checks for singularity +template +VectorType +applyUMFPack(const MatrixType& D, + std::shared_ptr>& linsolver, + VectorType x) +{ +#if HAVE_UMFPACK + if (!linsolver) + { + linsolver = std::make_shared>(D, 0); + } + + // The copy of x seems mandatory for calling UMFPack! + VectorType y(x.size()); + y = 0.; + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult res; + + // Solve + linsolver->apply(y, x, res); + + // Checking if there is any inf or nan in y + // it will be the solution before we find a way to catch the singularity of the matrix + for (size_t i_block = 0; i_block < y.size(); ++i_block) { + for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) { + if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) { + const std::string msg{"nan or inf value found after UMFPack solve due to singular matrix"}; + OpmLog::debug(msg); + OPM_THROW_NOLOG(NumericalIssue, msg); + } + } + } + return y; +#else + // this is not thread safe + OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. " + "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); +#endif // HAVE_UMFPACK +} + +template +Dune::Matrix +invertWithUMFPack(const MatrixType& D, std::shared_ptr >& linsolver) +{ +#if HAVE_UMFPACK + const int sz = D.M(); + const int bsz = D[0][0].M(); + VectorType e(sz); + e = 0.0; + + // Make a full block matrix. + Dune::Matrix inv(sz, sz); + + // Create inverse by passing basis vectors to the solver. + for (int ii = 0; ii < sz; ++ii) { + for (int jj = 0; jj < bsz; ++jj) { + e[ii][jj] = 1.0; + auto col = applyUMFPack(D, linsolver, e); + for (int cc = 0; cc < sz; ++cc) { + for (int dd = 0; dd < bsz; ++dd) { + inv[cc][ii][dd][jj] = col[cc][dd]; + } + } + e[ii][jj] = 0.0; + } + } + + return inv; +#else + // this is not thread safe + OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. " + "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); +#endif // HAVE_UMFPACK +} + +template +VectorType +invDX(const MatrixType& D, VectorType x, DeferredLogger& deferred_logger) +{ + // the function will change the value of x, so we should not use reference of x here. + + // TODO: store some of the following information to avoid to call it again and again for + // efficiency improvement. + // Bassically, only the solve / apply step is different. + + VectorType y(x.size()); + y = 0.; + + Dune::MatrixAdapter linearOperator(D); + + // Sequential incomplete LU decomposition as the preconditioner +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) + Dune::SeqILU preconditioner(D, 1.0); +#else + Dune::SeqILU0 preconditioner(D, 1.0); +#endif + // Dune::SeqILUn preconditioner(D, 1, 0.92); + // Dune::SeqGS preconditioner(D, 1, 1); + // Dune::SeqJac preconditioner(D, 1, 1); + + // Preconditioned BICGSTAB solver + Dune::BiCGSTABSolver linsolver(linearOperator, + preconditioner, + 1.e-8, // desired residual reduction factor + 250, // maximum number of iterations + 0); // verbosity of the solver */ + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult res; + + // Solve + linsolver.apply(y, x, res); + + if ( !res.converged ) { + OPM_DEFLOG_THROW(NumericalIssue, "the invDX did not converge ", deferred_logger); + } + + return y; +} + +template +ValueType calculateFrictionFactor(const double area, const double diameter, + const ValueType& w, const double roughness, + const ValueType& mu) +{ + + ValueType f = 0.; + // Reynolds number + const ValueType re = abs( diameter * w / (area * mu)); + + if ( re == 0.0 ) { + // make sure it is because the mass rate is zero + assert(w == 0.); + return 0.0; + } + + const ValueType re_value1 = 2000.; + const ValueType re_value2 = 4000.; + + if (re < re_value1) { + f = 16. / re; + } else if (re > re_value2){ + f = haalandFormular(re, diameter, roughness); + } else { // in between + const ValueType f1 = 16. / re_value1; + const ValueType f2 = haalandFormular(re_value2, diameter, roughness); + + f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1; + } + return f; +} + +template +ValueType frictionPressureLoss(const double l, const double diameter, + const double area, const double roughness, + const ValueType& density, + const ValueType& w, const ValueType& mu) +{ + const ValueType f = calculateFrictionFactor(area, diameter, w, roughness, mu); + // \Note: a factor of 2 needs to be here based on the dimensional analysis + return 2. * f * l * w * w / (area * area * diameter * density); +} + +template +ValueType valveContrictionPressureLoss(const ValueType& mass_rate, + const ValueType& density, + const double area_con, const double cv) +{ + // the formulation is adjusted a little bit for convinience + // velocity = mass_rate / (density * area) is applied to the original formulation + const double area = (area_con > 1.e-10 ? area_con : 1.e-10); + return mass_rate * mass_rate / (2. * density * cv * cv * area * area); +} + +template +ValueType velocityHead(const double area, const ValueType& mass_rate, + const ValueType& density) +{ + // \Note: a factor of 2 is added to the formulation in order to match results from the + // reference simulator. This is inline with what is done for the friction loss. + return (mass_rate * mass_rate / (area * area * density)); +} + +template +ValueType emulsionViscosity(const ValueType& water_fraction, + const ValueType& water_viscosity, + const ValueType& oil_fraction, + const ValueType& oil_viscosity, + const SICD& sicd) +{ + const double width_transition = sicd.widthTransitionRegion(); + + // it is just for now, we should be able to treat it. + if (width_transition <= 0.) { + OPM_THROW(std::runtime_error, "Not handling non-positive transition width now"); + } + + const double critical_value = sicd.criticalValue(); + const ValueType transition_start_value = critical_value - width_transition / 2.0; + const ValueType transition_end_value = critical_value + width_transition / 2.0; + + const ValueType liquid_fraction = water_fraction + oil_fraction; + // if there is no liquid, we just return zero + if (liquid_fraction == 0.) { + return 0.; + } + + const ValueType water_liquid_fraction = water_fraction / liquid_fraction; + + const double max_visco_ratio = sicd.maxViscosityRatio(); + if (water_liquid_fraction <= transition_start_value) { + return WIOEmulsionViscosity(oil_viscosity, water_liquid_fraction, max_visco_ratio); + } else if (water_liquid_fraction >= transition_end_value) { + return OIWEmulsionViscosity(water_viscosity, water_liquid_fraction, max_visco_ratio); + } else { // in the transition region + const ValueType viscosity_start_transition = WIOEmulsionViscosity(oil_viscosity, transition_start_value, max_visco_ratio); + const ValueType viscosity_end_transition = OIWEmulsionViscosity(water_viscosity, transition_end_value, max_visco_ratio); + const ValueType emulsion_viscosity = (viscosity_start_transition * (transition_end_value - water_liquid_fraction) + + viscosity_end_transition * (water_liquid_fraction - transition_start_value) ) / width_transition; + return emulsion_viscosity; + } +} + +template +using Vec = Dune::BlockVector>; +template +using Mat = Dune::BCRSMatrix>; + +#define INSTANCE_UMF(Dim) \ + template Vec applyUMFPack,Vec>(const Mat&, \ + std::shared_ptr>>&, \ + Vec); \ + template Dune::Matrix::block_type> \ + invertWithUMFPack,Vec>(const Mat& D, \ + std::shared_ptr>>&); + +INSTANCE_UMF(2) +INSTANCE_UMF(3) +INSTANCE_UMF(4) + +#define INSTANCE_IMPL(...) \ + template __VA_ARGS__ \ + frictionPressureLoss<__VA_ARGS__>(const double, \ + const double, \ + const double, \ + const double, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&); \ + template __VA_ARGS__ \ + valveContrictionPressureLoss<__VA_ARGS__>(const __VA_ARGS__& mass_rate, \ + const __VA_ARGS__& density, \ + const double, const double); \ + template __VA_ARGS__ \ + velocityHead<__VA_ARGS__>(const double, const __VA_ARGS__&, const __VA_ARGS__&); \ + template __VA_ARGS__ \ + emulsionViscosity<__VA_ARGS__>(const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const __VA_ARGS__&, \ + const SICD&); + +#define INSTANCE_EVAL(Dim) \ + INSTANCE_IMPL(DenseAd::Evaluation) + +INSTANCE_EVAL(3) +INSTANCE_EVAL(4) +INSTANCE_EVAL(5) +INSTANCE_EVAL(6) +INSTANCE_EVAL(7) +INSTANCE_EVAL(8) +INSTANCE_EVAL(9) + +} // namespace mswellhelpers + +} // namespace Opm diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index bf1462268..695dcd95d 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -23,203 +23,49 @@ #ifndef OPM_MSWELLHELPERS_HEADER_INCLUDED #define OPM_MSWELLHELPERS_HEADER_INCLUDED -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include -#if HAVE_UMFPACK -#include -#endif // HAVE_UMFPACK -#include +namespace Dune { +template class UMFPack; +} + +#include namespace Opm { +class DeferredLogger; +class SICD; + namespace mswellhelpers { /// Applies umfpack and checks for singularity template VectorType - applyUMFPack(const MatrixType& D, std::shared_ptr >& linsolver, VectorType x) - { -#if HAVE_UMFPACK - if (!linsolver) - { - linsolver.reset(new Dune::UMFPack(D, 0)); - } - - // The copy of x seems mandatory for calling UMFPack! - VectorType y(x.size()); - y = 0.; - - // Object storing some statistics about the solving process - Dune::InverseOperatorResult res; - - // Solve - linsolver->apply(y, x, res); - - // Checking if there is any inf or nan in y - // it will be the solution before we find a way to catch the singularity of the matrix - for (size_t i_block = 0; i_block < y.size(); ++i_block) { - for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) { - if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) { - const std::string msg{"nan or inf value found after UMFPack solve due to singular matrix"}; - OpmLog::debug(msg); - OPM_THROW_NOLOG(NumericalIssue, msg); - } - } - } - return y; -#else - // this is not thread safe - OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. " - "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); -#endif // HAVE_UMFPACK - } + applyUMFPack(const MatrixType& D, + std::shared_ptr>& linsolver, + VectorType x); /// Applies umfpack and checks for singularity template Dune::Matrix - invertWithUMFPack(const MatrixType& D, std::shared_ptr >& linsolver) - { -#if HAVE_UMFPACK - const int sz = D.M(); - const int bsz = D[0][0].M(); - VectorType e(sz); - e = 0.0; - - // Make a full block matrix. - Dune::Matrix inv(sz, sz); - - // Create inverse by passing basis vectors to the solver. - for (int ii = 0; ii < sz; ++ii) { - for (int jj = 0; jj < bsz; ++jj) { - e[ii][jj] = 1.0; - auto col = applyUMFPack(D, linsolver, e); - for (int cc = 0; cc < sz; ++cc) { - for (int dd = 0; dd < bsz; ++dd) { - inv[cc][ii][dd][jj] = col[cc][dd]; - } - } - e[ii][jj] = 0.0; - } - } - - return inv; -#else - // this is not thread safe - OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. " - "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); -#endif // HAVE_UMFPACK - } + invertWithUMFPack(const MatrixType& D, + std::shared_ptr >& linsolver); // obtain y = D^-1 * x with a BICSSTAB iterative solver template VectorType - invDX(const MatrixType& D, VectorType x, DeferredLogger& deferred_logger) - { - // the function will change the value of x, so we should not use reference of x here. - - // TODO: store some of the following information to avoid to call it again and again for - // efficiency improvement. - // Bassically, only the solve / apply step is different. - - VectorType y(x.size()); - y = 0.; - - Dune::MatrixAdapter linearOperator(D); - - // Sequential incomplete LU decomposition as the preconditioner -#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - Dune::SeqILU preconditioner(D, 1.0); -#else - Dune::SeqILU0 preconditioner(D, 1.0); -#endif - // Dune::SeqILUn preconditioner(D, 1, 0.92); - // Dune::SeqGS preconditioner(D, 1, 1); - // Dune::SeqJac preconditioner(D, 1, 1); - - // Preconditioned BICGSTAB solver - Dune::BiCGSTABSolver linsolver(linearOperator, - preconditioner, - 1.e-8, // desired residual reduction factor - 250, // maximum number of iterations - 0); // verbosity of the solver */ - - // Object storing some statistics about the solving process - Dune::InverseOperatorResult res; - - // Solve - linsolver.apply(y, x, res); - - if ( !res.converged ) { - OPM_DEFLOG_THROW(NumericalIssue, "the invDX did not converge ", deferred_logger); - } - - return y; - } - - + invDX(const MatrixType& D, VectorType x, DeferredLogger& deferred_logger); template - inline ValueType haalandFormular(const ValueType& re, const double diameter, const double roughness) - { - const ValueType value = -3.6 * log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) ); - - // sqrt(1/f) should be non-positive - assert(value >= 0.0); - - return 1. / (value * value); - } - - - - - template - inline ValueType calculateFrictionFactor(const double area, const double diameter, - const ValueType& w, const double roughness, const ValueType& mu) - { - - ValueType f = 0.; - // Reynolds number - const ValueType re = abs( diameter * w / (area * mu)); - - if ( re == 0.0 ) { - // make sure it is because the mass rate is zero - assert(w == 0.); - return 0.0; - } - - const ValueType re_value1 = 2000.; - const ValueType re_value2 = 4000.; - - if (re < re_value1) { - f = 16. / re; - } else if (re > re_value2){ - f = haalandFormular(re, diameter, roughness); - } else { // in between - const ValueType f1 = 16. / re_value1; - const ValueType f2 = haalandFormular(re_value2, diameter, roughness); - - f = (f2 - f1) / (re_value2 - re_value1) * (re - re_value1) + f1; - } - return f; - } - - - - + ValueType calculateFrictionFactor(const double area, const double diameter, + const ValueType& w, const double roughness, + const ValueType& mu); // calculating the friction pressure loss @@ -231,116 +77,32 @@ namespace mswellhelpers // roughness is the absolute roughness // mu is the average phase viscosity template - ValueType frictionPressureLoss(const double l, const double diameter, const double area, const double roughness, - const ValueType& density, const ValueType& w, const ValueType& mu) - { - const ValueType f = calculateFrictionFactor(area, diameter, w, roughness, mu); - // \Note: a factor of 2 needs to be here based on the dimensional analysis - return 2. * f * l * w * w / (area * area * diameter * density); - } + ValueType frictionPressureLoss(const double l, const double diameter, + const double area, const double roughness, + const ValueType& density, + const ValueType& w, const ValueType& mu); template - ValueType valveContrictionPressureLoss(const ValueType& mass_rate, const ValueType& density, - const double area_con, const double cv) - { - // the formulation is adjusted a little bit for convinience - // velocity = mass_rate / (density * area) is applied to the original formulation - const double area = (area_con > 1.e-10 ? area_con : 1.e-10); - return mass_rate * mass_rate / (2. * density * cv * cv * area * area); - } + ValueType valveContrictionPressureLoss(const ValueType& mass_rate, + const ValueType& density, + const double area_con, const double cv); template - ValueType velocityHead(const double area, const ValueType& mass_rate, const ValueType& density) - { - // \Note: a factor of 2 is added to the formulation in order to match results from the - // reference simulator. This is inline with what is done for the friction loss. - return (mass_rate * mass_rate / (area * area * density)); - } - - - - // water in oil emulsion viscosity - // TODO: maybe it should be two different ValueTypes. When we calculate the viscosity for transitional zone - template - ValueType WIOEmulsionViscosity(const ValueType& oil_viscosity, const ValueType& water_liquid_fraction, - const double max_visco_ratio) - { - const ValueType temp_value = 1. / (1. - (0.8415 / 0.7480 * water_liquid_fraction) ); - const ValueType viscosity_ratio = pow(temp_value, 2.5); - - if (viscosity_ratio <= max_visco_ratio) { - return oil_viscosity * viscosity_ratio; - } else { - return oil_viscosity * max_visco_ratio; - } - } - - - - - - // oil in water emulsion viscosity - template - ValueType OIWEmulsionViscosity(const ValueType& water_viscosity, const ValueType& water_liquid_fraction, - const double max_visco_ratio) - { - const ValueType temp_value = 1. / (1. - (0.6019 / 0.6410) * (1. - water_liquid_fraction) ); - const ValueType viscosity_ratio = pow(temp_value, 2.5); - - if (viscosity_ratio <= max_visco_ratio) { - return water_viscosity * viscosity_ratio; - } else { - return water_viscosity * max_visco_ratio; - } - } - - - + ValueType velocityHead(const double area, const ValueType& mass_rate, + const ValueType& density); // calculating the viscosity of oil-water emulsion at local conditons template - ValueType emulsionViscosity(const ValueType& water_fraction, const ValueType& water_viscosity, - const ValueType& oil_fraction, const ValueType& oil_viscosity, - const SICD& sicd) - { - const double width_transition = sicd.widthTransitionRegion(); - - // it is just for now, we should be able to treat it. - if (width_transition <= 0.) { - OPM_THROW(std::runtime_error, "Not handling non-positive transition width now"); - } - - const double critical_value = sicd.criticalValue(); - const ValueType transition_start_value = critical_value - width_transition / 2.0; - const ValueType transition_end_value = critical_value + width_transition / 2.0; - - const ValueType liquid_fraction = water_fraction + oil_fraction; - // if there is no liquid, we just return zero - if (liquid_fraction == 0.) { - return 0.; - } - - const ValueType water_liquid_fraction = water_fraction / liquid_fraction; - - const double max_visco_ratio = sicd.maxViscosityRatio(); - if (water_liquid_fraction <= transition_start_value) { - return WIOEmulsionViscosity(oil_viscosity, water_liquid_fraction, max_visco_ratio); - } else if(water_liquid_fraction >= transition_end_value) { - return OIWEmulsionViscosity(water_viscosity, water_liquid_fraction, max_visco_ratio); - } else { // in the transition region - const ValueType viscosity_start_transition = WIOEmulsionViscosity(oil_viscosity, transition_start_value, max_visco_ratio); - const ValueType viscosity_end_transition = OIWEmulsionViscosity(water_viscosity, transition_end_value, max_visco_ratio); - const ValueType emulsion_viscosity = (viscosity_start_transition * (transition_end_value - water_liquid_fraction) - + viscosity_end_transition * (water_liquid_fraction - transition_start_value) ) / width_transition; - return emulsion_viscosity; - } - } + ValueType emulsionViscosity(const ValueType& water_fraction, + const ValueType& water_viscosity, + const ValueType& oil_fraction, + const ValueType& oil_viscosity, + const SICD& sicd); } // namespace mswellhelpers - -} +} // namespace Opm #endif diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index e9a3eafd4..79970efbe 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -21,6 +21,8 @@ #include #include +#include + #include #include diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index 48a444550..6e6be6b49 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -32,11 +32,14 @@ #include #include #include -#include #include #include +namespace Dune { +template class UMFPack; +} + namespace Opm {