OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-07 04:53:03 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5376 from akva2/various_template_scalar

Browse Source 
Various classes: use Scalar type
This commit is contained in:
Bård Skaflestad 2024-05-23 09:48:33 +02:00 committed by GitHub
commit 8edcc6a060
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
25 changed files with 448 additions and 384 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
opm/simulators/flow/BlackoilModel.hpp
View File

@ -356,7 +356,7 @@ namespace Opm {
}
// ----------- Check if converged -----------
std::vector<double> residual_norms;
std::vector<Scalar> residual_norms;
perfTimer.reset();
perfTimer.start();
// the step is not considered converged until at least minIter iterations is done
@ -527,7 +527,7 @@ namespace Opm {
}
// compute the "relative" change of the solution between time steps
double relativeChange() const
Scalar relativeChange() const
{
Scalar resultDelta = 0.0;
Scalar resultDenom = 0.0;
@ -714,17 +714,17 @@ namespace Opm {
return terminal_output_;
}
std::tuple<double,double> convergenceReduction(Parallel::Communication comm,
const double pvSumLocal,
const double numAquiferPvSumLocal,
std::tuple<Scalar,Scalar> convergenceReduction(Parallel::Communication comm,
const Scalar pvSumLocal,
const Scalar numAquiferPvSumLocal,
std::vector< Scalar >& R_sum,
std::vector< Scalar >& maxCoeff,
std::vector< Scalar >& B_avg)
{
OPM_TIMEBLOCK(convergenceReduction);
// Compute total pore volume (use only owned entries)
double pvSum = pvSumLocal;
double numAquiferPvSum = numAquiferPvSumLocal;
Scalar pvSum = pvSumLocal;
Scalar numAquiferPvSum = numAquiferPvSumLocal;
if( comm.size() > 1 )
{
@ -777,14 +777,14 @@ namespace Opm {
/// \brief Get reservoir quantities on this process needed for convergence calculations.
/// \return A pair of the local pore volume of interior cells and the pore volumes
/// of the cells associated with a numerical aquifer.
std::pair<double,double> localConvergenceData(std::vector<Scalar>& R_sum,
std::pair<Scalar,Scalar> localConvergenceData(std::vector<Scalar>& R_sum,
std::vector<Scalar>& maxCoeff,
std::vector<Scalar>& B_avg,
std::vector<int>& maxCoeffCell)
{
OPM_TIMEBLOCK(localConvergenceData);
double pvSumLocal = 0.0;
double numAquiferPvSumLocal = 0.0;
Scalar pvSumLocal = 0.0;
Scalar numAquiferPvSumLocal = 0.0;
const auto& model = simulator_.model();
const auto& problem = simulator_.problem();
@ -830,10 +830,10 @@ namespace Opm {
/// \brief Compute the total pore volume of cells violating CNV that are not part
/// of a numerical aquifer.
double computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
Scalar computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
{
OPM_TIMEBLOCK(computeCnvErrorPv);
double errorPV{};
Scalar errorPV{};
const auto& model = simulator_.model();
const auto& problem = simulator_.problem();
const auto& residual = simulator_.model().linearizer().residual();
@ -852,7 +852,7 @@ namespace Opm {
elemCtx.updatePrimaryStencil(elem);
// elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const double pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) * model.dofTotalVolume(cell_idx);
const Scalar pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) * model.dofTotalVolume(cell_idx);
const auto& cellResidual = residual[cell_idx];
bool cnvViolated = false;
@ -939,8 +939,8 @@ namespace Opm {
OpmLog::debug(message);
}
}
const double tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
const double tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
const Scalar tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
const Scalar tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
// Finish computation
std::vector<Scalar> CNV(numComp);
@ -958,10 +958,10 @@ namespace Opm {
using CR = ConvergenceReport;
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
const double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
const Scalar res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
const CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
CR::ReservoirFailure::Type::Cnv };
const double tol[2] = { tol_mb, tol_cnv };
const Scalar tol[2] = { tol_mb, tol_cnv };
for (int ii : {0, 1}) {
if (std::isnan(res[ii])) {
@ -1031,7 +1031,7 @@ namespace Opm {
ConvergenceReport getConvergence(const SimulatorTimerInterface& timer,
const int iteration,
const int maxIter,
std::vector<double>& residual_norms)
std::vector<Scalar>& residual_norms)
{
OPM_TIMEBLOCK(getConvergence);
// Get convergence reports for reservoir and wells.
@ -1055,20 +1055,20 @@ namespace Opm {
/// Wrapper required due to not following generic API
template<class T>
std::vector<std::vector<double> >
std::vector<std::vector<Scalar> >
computeFluidInPlace(const T&, const std::vector<int>& fipnum) const
{
return computeFluidInPlace(fipnum);
}
/// Should not be called
std::vector<std::vector<double> >
std::vector<std::vector<Scalar> >
computeFluidInPlace(const std::vector<int>& /*fipnum*/) const
{
OPM_TIMEBLOCK(computeFluidInPlace);
//assert(true)
//return an empty vector
std::vector<std::vector<double> > regionValues(0, std::vector<double>(0,0.0));
std::vector<std::vector<Scalar> > regionValues(0, std::vector<Scalar>(0,0.0));
return regionValues;
}
@ -1148,8 +1148,8 @@ namespace Opm {
/// \brief The number of cells of the global grid.
long int global_nc_;
std::vector<std::vector<double>> residual_norms_history_;
double current_relaxation_;
std::vector<std::vector<Scalar>> residual_norms_history_;
Scalar current_relaxation_;
BVector dx_old_;
std::vector<StepReport> convergence_reports_;
@ -1198,7 +1198,7 @@ namespace Opm {
const auto R2 = modelResid[cell_idx][compIdx];
R_sum[compIdx] += R2;
const double Rval = std::abs(R2) / pvValue;
const Scalar Rval = std::abs(R2) / pvValue;
if (Rval > maxCoeff[compIdx]) {
maxCoeff[compIdx] = Rval;
maxCoeffCell[compIdx] = cell_idx;
@ -1304,10 +1304,10 @@ namespace Opm {
}
private:
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; }
double maxResidualAllowed() const { return param_.max_residual_allowed_; }
Scalar dpMaxRel() const { return param_.dp_max_rel_; }
Scalar dsMax() const { return param_.ds_max_; }
Scalar drMaxRel() const { return param_.dr_max_rel_; }
Scalar maxResidualAllowed() const { return param_.max_residual_allowed_; }
double linear_solve_setup_time_;
public:

42
opm/simulators/flow/BlackoilModelNldd.hpp
View File

@ -389,7 +389,7 @@ private:
}
detailTimer.reset();
detailTimer.start();
std::vector<double> resnorms;
std::vector<Scalar> resnorms;
auto convreport = this->getDomainConvergence(domain, timer, 0, logger, resnorms);
if (convreport.converged()) {
// TODO: set more info, timing etc.
@ -532,7 +532,7 @@ private:
}
//! \brief Get reservoir quantities on this process needed for convergence calculations.
std::pair<double, double> localDomainConvergenceData(const Domain& domain,
std::pair<Scalar, Scalar> localDomainConvergenceData(const Domain& domain,
std::vector<Scalar>& R_sum,
std::vector<Scalar>& maxCoeff,
std::vector<Scalar>& B_avg,
@ -540,8 +540,8 @@ private:
{
const auto& modelSimulator = model_.simulator();
double pvSumLocal = 0.0;
double numAquiferPvSumLocal = 0.0;
Scalar pvSumLocal = 0.0;
Scalar numAquiferPvSumLocal = 0.0;
const auto& model = modelSimulator.model();
const auto& problem = modelSimulator.problem();
@ -617,12 +617,12 @@ private:
iteration >= model_.param().min_strict_cnv_iter_;
// Tighter bound for local convergence should increase the
// likelyhood of: local convergence => global convergence
const double tol_cnv = model_.param().local_tolerance_scaling_cnv_
const Scalar tol_cnv = model_.param().local_tolerance_scaling_cnv_
* (use_relaxed_cnv ? model_.param().tolerance_cnv_relaxed_
: model_.param().tolerance_cnv_);
const bool use_relaxed_mb = iteration >= model_.param().min_strict_mb_iter_;
const double tol_mb = model_.param().local_tolerance_scaling_mb_
const Scalar tol_mb = model_.param().local_tolerance_scaling_mb_
* (use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
// Finish computation
@ -639,10 +639,10 @@ private:
ConvergenceReport report{reportTime};
using CR = ConvergenceReport;
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
Scalar res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
CR::ReservoirFailure::Type::Cnv };
double tol[2] = { tol_mb, tol_cnv };
Scalar tol[2] = { tol_mb, tol_cnv };
for (int ii : {0, 1}) {
if (std::isnan(res[ii])) {
report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
@ -704,7 +704,7 @@ private:
const SimulatorTimerInterface& timer,
const int iteration,
DeferredLogger& logger,
std::vector<double>& residual_norms)
std::vector<Scalar>& residual_norms)
{
std::vector<Scalar> B_avg(numEq, 0.0);
auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(),
@ -732,16 +732,16 @@ private:
return domain_order;
} else if (model_.param().local_solve_approach_ == DomainSolveApproach::GaussSeidel) {
// Calculate the measure used to order the domains.
std::vector<double> measure_per_domain(domains_.size());
std::vector<Scalar> measure_per_domain(domains_.size());
switch (model_.param().local_domain_ordering_) {
case DomainOrderingMeasure::AveragePressure: {
// Use average pressures to order domains.
for (const auto& domain : domains_) {
double press_sum = 0.0;
Scalar press_sum = 0.0;
for (const int c : domain.cells) {
press_sum += solution[c][Indices::pressureSwitchIdx];
}
const double avgpress = press_sum / domain.cells.size();
const Scalar avgpress = press_sum / domain.cells.size();
measure_per_domain[domain.index] = avgpress;
}
break;
@ -749,7 +749,7 @@ private:
case DomainOrderingMeasure::MaxPressure: {
// Use max pressures to order domains.
for (const auto& domain : domains_) {
double maxpress = 0.0;
Scalar maxpress = 0.0;
for (const int c : domain.cells) {
maxpress = std::max(maxpress, solution[c][Indices::pressureSwitchIdx]);
}
@ -762,7 +762,7 @@ private:
const auto& residual = modelSimulator.model().linearizer().residual();
const int num_vars = residual[0].size();
for (const auto& domain : domains_) {
double maxres = 0.0;
Scalar maxres = 0.0;
for (const int c : domain.cells) {
for (int ii = 0; ii < num_vars; ++ii) {
maxres = std::max(maxres, std::fabs(residual[c][ii]));
@ -829,8 +829,8 @@ private:
// We do not accept a solution if the wells are unconverged.
if (!convrep.wellFailed()) {
// Calculare the sums of the mb and cnv failures.
double mb_sum = 0.0;
double cnv_sum = 0.0;
Scalar mb_sum = 0.0;
Scalar cnv_sum = 0.0;
for (const auto& rc : convrep.reservoirConvergence()) {
if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) {
mb_sum += rc.value();
@ -839,8 +839,8 @@ private:
}
}
// If not too high, we overrule the convergence failure.
const double acceptable_local_mb_sum = 1e-3;
const double acceptable_local_cnv_sum = 1.0;
const Scalar acceptable_local_mb_sum = 1e-3;
const Scalar acceptable_local_cnv_sum = 1.0;
if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) {
local_report.converged = true;
logger.debug(fmt::format("Accepting solution in unconverged domain {} on rank {}.", domain.index, rank_));
@ -865,17 +865,17 @@ private:
}
}
double computeCnvErrorPvLocal(const Domain& domain,
Scalar computeCnvErrorPvLocal(const Domain& domain,
const std::vector<Scalar>& B_avg, double dt) const
{
double errorPV{};
Scalar errorPV{};
const auto& simulator = model_.simulator();
const auto& model = simulator.model();
const auto& problem = simulator.problem();
const auto& residual = simulator.model().linearizer().residual();
for (const int cell_idx : domain.cells) {
const double pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
const Scalar pvValue = problem.referencePorosity(cell_idx, /*timeIdx=*/0) *
model.dofTotalVolume(cell_idx);
const auto& cellResidual = residual[cell_idx];
bool cnvViolated = false;

40
opm/simulators/flow/BlackoilModelParameters.hpp
View File

@ -488,37 +488,37 @@ namespace Opm
public:
/// Max relative change in bhp in single iteration.
double dbhp_max_rel_;
Scalar dbhp_max_rel_;
/// Max absolute change in well volume fraction in single iteration.
double dwell_fraction_max_;
Scalar dwell_fraction_max_;
/// Absolute max limit for residuals.
double max_residual_allowed_;
Scalar max_residual_allowed_;
//// Max allowed pore volume faction where CNV is violated. Below the
//// relaxed tolerance tolerance_cnv_relaxed_ is used.
double relaxed_max_pv_fraction_;
Scalar relaxed_max_pv_fraction_;
/// Relative mass balance tolerance (total mass balance error).
double tolerance_mb_;
Scalar tolerance_mb_;
/// Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).
double tolerance_mb_relaxed_;
Scalar tolerance_mb_relaxed_;
/// Local convergence tolerance (max of local saturation errors).
double tolerance_cnv_;
Scalar tolerance_cnv_;
/// Relaxed local convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV < relaxed_max_pv_fraction_).
double tolerance_cnv_relaxed_;
Scalar tolerance_cnv_relaxed_;
/// Well convergence tolerance.
double tolerance_wells_;
Scalar tolerance_wells_;
/// Tolerance for the well control equations
// TODO: it might need to distinguish between rate control and pressure control later
double tolerance_well_control_;
Scalar tolerance_well_control_;
/// Tolerance for the pressure equations for multisegment wells
double tolerance_pressure_ms_wells_;
Scalar tolerance_pressure_ms_wells_;
/// Relaxed tolerance for for the well flow residual
double relaxed_tolerance_flow_well_;
Scalar relaxed_tolerance_flow_well_;
/// Relaxed tolerance for the MSW pressure solution
double relaxed_tolerance_pressure_ms_well_;
Scalar relaxed_tolerance_pressure_ms_well_;
/// Maximum pressure change over an iteratio for ms wells
double max_pressure_change_ms_wells_;
Scalar max_pressure_change_ms_wells_;
/// Maximum inner iteration number for ms wells
int max_inner_iter_ms_wells_;
@ -530,7 +530,7 @@ namespace Opm
int strict_outer_iter_wells_;
/// Regularization factor for wells
double regularization_factor_wells_;
Scalar regularization_factor_wells_;
/// Maximum newton iterations with inner well iterations
int max_niter_inner_well_iter_;
@ -546,7 +546,7 @@ namespace Opm
/// Tolerance for time step in seconds where single precision can be used
/// for solving for the Jacobian
double maxSinglePrecisionTimeStep_;
Scalar maxSinglePrecisionTimeStep_;
/// Minimum number of Newton iterations before we can use relaxed CNV convergence criterion
int min_strict_cnv_iter_;
@ -606,12 +606,12 @@ namespace Opm
int max_local_solve_iterations_;
double local_tolerance_scaling_mb_;
double local_tolerance_scaling_cnv_;
Scalar local_tolerance_scaling_mb_;
Scalar local_tolerance_scaling_cnv_;
int nldd_num_initial_newton_iter_{1};
int num_local_domains_{0};
double local_domain_partition_imbalance_{1.03};
Scalar local_domain_partition_imbalance_{1.03};
std::string local_domain_partition_method_;
DomainOrderingMeasure local_domain_ordering_{DomainOrderingMeasure::MaxPressure};
@ -671,7 +671,7 @@ namespace Opm
local_tolerance_scaling_cnv_ = Parameters::get<TypeTag, Properties::LocalToleranceScalingCnv>();
nldd_num_initial_newton_iter_ = Parameters::get<TypeTag, Properties::NlddNumInitialNewtonIter>();
num_local_domains_ = Parameters::get<TypeTag, Properties::NumLocalDomains>();
local_domain_partition_imbalance_ = std::max(1.0, Parameters::get<TypeTag, Properties::LocalDomainsPartitioningImbalance>());
local_domain_partition_imbalance_ = std::max(Scalar{1.0}, Parameters::get<TypeTag, Properties::LocalDomainsPartitioningImbalance>());
local_domain_partition_method_ = Parameters::get<TypeTag, Properties::LocalDomainsPartitioningMethod>();
deck_file_name_ = Parameters::get<TypeTag, Properties::EclDeckFileName>();
network_max_strict_iterations_ = Parameters::get<TypeTag, Properties::NetworkMaxStrictIterations>();

10
opm/simulators/flow/GenericTracerModel.hpp
View File

@ -77,7 +77,7 @@ public:
/*!
* \brief Return well tracer rates
*/
const std::map<std::pair<std::string, std::string>, double>&
const std::map<std::pair<std::string, std::string>, Scalar>&
getWellTracerRates() const {return wellTracerRate_;}
template<class Serializer>
@ -105,9 +105,11 @@ protected:
bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b);
bool linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b);
bool linearSolveBatchwise_(const TracerMatrix& M,
std::vector<TracerVector>& x,
std::vector<TracerVector>& b);
double currentConcentration_(const Well& eclWell, const std::string& name) const;
Scalar currentConcentration_(const Well& eclWell, const std::string& name) const;
const GridView& gridView_;
const EclipseState& eclState_;
@ -120,7 +122,7 @@ protected:
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
// <wellName, tracerIdx> -> wellRate
std::map<std::pair<std::string, std::string>, double> wellTracerRate_;
std::map<std::pair<std::string, std::string>, Scalar> wellTracerRate_;
/// \brief Function returning the cell centers
std::function<std::array<double,dimWorld>(int)> centroids_;
};

2
opm/simulators/flow/GenericTracerModel_impl.hpp
View File

@ -144,7 +144,7 @@ fname(int tracerIdx) const
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
double GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
currentConcentration_(const Well& eclWell, const std::string& name) const
{
return eclWell.getTracerProperties().getConcentration(name);

42
opm/simulators/flow/NonlinearSolver.cpp
View File

@ -31,8 +31,9 @@
namespace Opm {
namespace detail {
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, const int numPhases, const double relaxRelTol,
template<class Scalar>
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, const int numPhases, const Scalar relaxRelTol,
bool& oscillate, bool& stagnate)
{
// The detection of oscillation in two primary variable results in the report of the detection
@ -52,8 +53,8 @@ void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const auto& F1 = residualHistory[it - 1];
const auto& F2 = residualHistory[it - 2];
for (int p = 0; p < numPhases; ++p) {
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
const Scalar d1 = std::abs((F0[p] - F2[p]) / F0[p]);
const Scalar d2 = std::abs((F0[p] - F1[p]) / F0[p]);
oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2);
@ -65,9 +66,9 @@ void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
oscillate = (oscillatePhase > 1);
}
template <class BVector>
template <class BVector, class Scalar>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const double omega,
const Scalar omega,
NonlinearRelaxType relaxType)
{
// The dxOld is updated with dx.
@ -104,16 +105,25 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
return;
}
template<int Size> using BV = Dune::BlockVector<Dune::FieldVector<double,Size>>;
#define INSTANCE(Size) \
template void stabilizeNonlinearUpdate<BV<Size>>(BV<Size>&, BV<Size>&, \
const double, NonlinearRelaxType);
INSTANCE(1)
INSTANCE(2)
INSTANCE(3)
INSTANCE(4)
INSTANCE(5)
INSTANCE(6)
template<class Scalar, int Size>
using BV = Dune::BlockVector<Dune::FieldVector<Scalar,Size>>;
#define INSTANCE(T,Size) \
template void stabilizeNonlinearUpdate<BV<T,Size>,T>(BV<T,Size>&, BV<T,Size>&, \
const T, NonlinearRelaxType);
#define INSTANCE_TYPE(T) \
template void detectOscillations(const std::vector<std::vector<T>>&, \
const int, const int, const T, \
bool&, bool&); \
INSTANCE(T,1) \
INSTANCE(T,2) \
INSTANCE(T,3) \
INSTANCE(T,4) \
INSTANCE(T,5) \
INSTANCE(T,6)
INSTANCE_TYPE(double)
} // namespace detail
} // namespace Opm

27
opm/simulators/flow/NonlinearSolver.hpp
View File

@ -89,15 +89,16 @@ enum class NonlinearRelaxType {
namespace detail {
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
const int it, const int numPhases, const double relaxRelTol,
template<class Scalar>
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, const int numPhases, const Scalar relaxRelTol,
bool& oscillate, bool& stagnate);
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>
template <class BVector, class Scalar>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
const double omega, NonlinearRelaxType relaxType);
const Scalar omega, NonlinearRelaxType relaxType);
}
@ -113,9 +114,9 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
struct SolverParameters
{
NonlinearRelaxType relaxType_;
double relaxMax_;
double relaxIncrement_;
double relaxRelTol_;
Scalar relaxMax_;
Scalar relaxIncrement_;
Scalar relaxRelTol_;
int maxIter_; // max nonlinear iterations
int minIter_; // min nonlinear iterations
@ -277,7 +278,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
int wellIterationsLastStep() const
{ return wellIterationsLast_; }
std::vector<std::vector<double> >
std::vector<std::vector<Scalar> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{ return model_->computeFluidInPlace(fipnum); }
@ -290,7 +291,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
{ return *model_; }
/// Detect oscillation or stagnation in a given residual history.
void detectOscillations(const std::vector<std::vector<double>>& residualHistory,
void detectOscillations(const std::vector<std::vector<Scalar>>& residualHistory,
const int it, bool& oscillate, bool& stagnate) const
{
detail::detectOscillations(residualHistory, it, model_->numPhases(),
@ -301,17 +302,17 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const Scalar omega) const
{
detail::stabilizeNonlinearUpdate(dx, dxOld, omega, this->relaxType());
}
/// The greatest relaxation factor (i.e. smallest factor) allowed.
double relaxMax() const
Scalar relaxMax() const
{ return param_.relaxMax_; }
/// The step-change size for the relaxation factor.
double relaxIncrement() const
Scalar relaxIncrement() const
{ return param_.relaxIncrement_; }
/// The relaxation type (Dampen or SOR).
@ -319,7 +320,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
{ return param_.relaxType_; }
/// The relaxation relative tolerance.
double relaxRelTol() const
Scalar relaxRelTol() const
{ return param_.relaxRelTol_; }
/// The maximum number of nonlinear iterations allowed.

10
opm/simulators/flow/OutputBlackoilModule.hpp
View File

@ -1607,9 +1607,9 @@ private:
sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/false);
}
const double sg = getValue(fs.saturation(gasPhaseIdx));
const double rhog = getValue(fs.density(gasPhaseIdx));
const double xgW = FluidSystem::phaseIsActive(waterPhaseIdx)
const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
const Scalar rhog = getValue(fs.density(gasPhaseIdx));
const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx)
? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
: FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex());
@ -1629,7 +1629,7 @@ private:
}
if (!this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty()) {
const Scalar mobileGas = massGas / mM * std::max(0.0, sg - sgcr);
const Scalar mobileGas = massGas / mM * std::max(Scalar{0.0}, sg - sgcr);
this->fip_[Inplace::Phase::CO2InGasPhaseMob][globalDofIdx] = mobileGas;
}
@ -1657,7 +1657,7 @@ private:
}
if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty()) {
const Scalar mobileMassGas = massGas * std::max(0.0, sg - sgcr);
const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - sgcr);
this->fip_[Inplace::Phase::CO2MassInGasPhaseMob][globalDofIdx] = mobileMassGas;
}

2
opm/simulators/flow/Transmissibility.hpp
View File

@ -227,7 +227,7 @@ protected:
void applyEditNncToGridTransHelper_(const std::unordered_map<std::size_t,int>& globalToLocal,
const std::string& keyword, const std::vector<NNCdata>& nncs,
const std::function<KeywordLocation(const NNCdata&)>& getLocation,
const std::function<void(double&, const double&)>& apply);
const std::function<void(Scalar&, const Scalar&)>& apply);
void extractPermeability_();

24
opm/simulators/flow/Transmissibility_impl.hpp
View File

@ -679,7 +679,13 @@ extractPorosity_()
// over several processes.)
const auto& fp = eclState_.fieldProps();
if (fp.has_double("PORO")) {
porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO");
if constexpr (std::is_same_v<Scalar,double>) {
porosity_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO");
} else {
const auto por = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"PORO");
porosity_.resize(por.size());
std::copy(por.begin(), por.end(), porosity_.begin());
}
}
else
throw std::logic_error("Can't read the porosityfrom the ecl state. "
@ -695,7 +701,13 @@ extractDispersion_()
"contains the DISPERC keyword.");
}
const auto& fp = eclState_.fieldProps();
dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC");
if constexpr (std::is_same_v<Scalar,double>) {
dispersion_ = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC");
} else {
const auto disp = this-> lookUpData_.assignFieldPropsDoubleOnLeaf(fp,"DISPERC");
dispersion_.resize(disp.size());
std::copy(disp.begin(), disp.end(), dispersion_.begin());
}
}
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@ -760,7 +772,7 @@ applyAllZMultipliers_(Scalar& trans,
auto multiplier = transMult.getMultiplier(cartElemIdx, FaceDir::ZPlus);
cartElemIdx += cartDims[0]*cartDims[1];
multiplier *= transMult.getMultiplier(cartElemIdx, FaceDir::ZMinus);
mult = std::min(mult, multiplier);
mult = std::min(mult, static_cast<Scalar>(multiplier));
}
}
@ -1105,7 +1117,7 @@ applyEditNncToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLoca
[&input](const NNCdata& nnc){
return input.edit_location(nnc);},
// Multiply transmissibility with EDITNNC value
[](double& trans, const double& rhs){ trans *= rhs;});
[](Scalar& trans, const Scalar& rhs){ trans *= rhs;});
}
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@ -1118,7 +1130,7 @@ applyEditNncrToGridTrans_(const std::unordered_map<std::size_t,int>& globalToLoc
[&input](const NNCdata& nnc){
return input.editr_location(nnc);},
// Replace Transmissibility with EDITNNCR value
[](double& trans, const double& rhs){ trans = rhs;});
[](Scalar& trans, const Scalar& rhs){ trans = rhs;});
}
template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@ -1127,7 +1139,7 @@ applyEditNncToGridTransHelper_(const std::unordered_map<std::size_t,int>& global
const std::string& keyword,
const std::vector<NNCdata>& nncs,
const std::function<KeywordLocation(const NNCdata&)>& getLocation,
const std::function<void(double&, const double&)>& apply)
const std::function<void(Scalar&, const Scalar&)>& apply)
{
if (nncs.empty())
return;

164
opm/simulators/linalg/WellOperators.hpp
View File

@ -32,9 +32,7 @@
#include <cstddef>
namespace Opm
{
namespace Opm {
//=====================================================================
// Implementation for ISTL-matrix based operators
@ -56,23 +54,27 @@ template <class X, class Y>
class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
{
public:
using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const = 0;
using field_type = typename X::field_type;
using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
virtual void addWellPressureEquations(PressureMatrix& jacobian,
const X& weights,
const bool use_well_weights) const = 0;
virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
virtual int getNumberOfExtraEquations() const = 0;
};
template <class WellModel, class X, class Y>
class WellModelAsLinearOperator : public Opm::LinearOperatorExtra<X, Y>
class WellModelAsLinearOperator : public LinearOperatorExtra<X, Y>
{
public:
using Base = Opm::LinearOperatorExtra<X, Y>;
using Base = LinearOperatorExtra<X, Y>;
using field_type = typename Base::field_type;
using PressureMatrix = typename Base::PressureMatrix;
explicit WellModelAsLinearOperator(const WellModel& wm)
: wellMod_(wm)
{
}
/*! \brief apply operator to x: \f$ y = A(x) \f$
The input vector is consistent and the output must also be
consistent on the interior+border partition.
@ -84,7 +86,7 @@ public:
}
//! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
virtual void applyscaleadd(field_type alpha, const X& x, Y& y) const override
void applyscaleadd(field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
wellMod_.applyScaleAdd(alpha, x, y);
@ -99,16 +101,21 @@ public:
{
return Dune::SolverCategory::sequential;
}
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const override
void addWellPressureEquations(PressureMatrix& jacobian,
const X& weights,
const bool use_well_weights) const override
{
OPM_TIMEBLOCK(addWellPressureEquations);
wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
{
OPM_TIMEBLOCK(addWellPressureEquationsStruct);
wellMod_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const override
{
return wellMod_.numLocalWellsEnd();
@ -133,84 +140,88 @@ template<class M, class X, class Y, bool overlapping >
class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
{
public:
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
using matrix_type = M;
using domain_type = X;
using range_type = Y;
using field_type = typename X::field_type;
using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
#else
typedef Dune::CollectiveCommunication< int > communication_type;
using communication_type = Dune::CollectiveCommunication<int>;
#endif
Dune::SolverCategory::Category category() const override
{
return overlapping ?
Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
}
Dune::SolverCategory::Category category() const override
{
return overlapping ?
Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
}
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A,
const Opm::LinearOperatorExtra<X, Y>& wellOper,
const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
: A_( A ), wellOper_( wellOper ), comm_(comm)
{}
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A,
const LinearOperatorExtra<X, Y>& wellOper,
const std::shared_ptr<communication_type>& comm = {})
: A_( A ), wellOper_( wellOper ), comm_(comm)
{}
void apply( const X& x, Y& y ) const override
{
OPM_TIMEBLOCK(apply);
A_.mv(x, y);
virtual void apply( const X& x, Y& y ) const override
{
OPM_TIMEBLOCK(apply);
A_.mv( x, y );
// add well model modification to y
wellOper_.apply(x, y);
// add well model modification to y
wellOper_.apply(x, y );
#if HAVE_MPI
if (comm_) {
comm_->project(y);
}
#endif
}
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
// y += \alpha * A * x
void applyscaleadd (field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
A_.usmv(alpha, x, y);
// y += \alpha * A * x
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
A_.usmv(alpha,x,y);
// add scaled well model modification to y
wellOper_.applyscaleadd(alpha, x, y);
// add scaled well model modification to y
wellOper_.applyscaleadd( alpha, x, y );
#if HAVE_MPI
if (comm_) {
comm_->project( y );
}
#endif
}
#if HAVE_MPI
if( comm_ )
comm_->project( y );
#endif
}
const matrix_type& getmat() const override { return A_; }
virtual const matrix_type& getmat() const override { return A_; }
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
void addWellPressureEquations(PressureMatrix& jacobian,
const X& weights,
const bool use_well_weights) const
{
OPM_TIMEBLOCK(addWellPressureEquations);
wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
OPM_TIMEBLOCK(addWellPressureEquations);
wellOper_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const
{
return wellOper_.getNumberOfExtraEquations();
}
protected:
const matrix_type& A_ ;
const Opm::LinearOperatorExtra<X, Y>& wellOper_;
std::shared_ptr< communication_type > comm_;
const matrix_type& A_ ;
const LinearOperatorExtra<X, Y>& wellOper_;
std::shared_ptr<communication_type> comm_;
};
/*!
\brief Adapter to combine a matrix and another linear operator into
a combined linear operator.
@ -223,18 +234,17 @@ template<class M, class X, class Y, bool overlapping >
class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
{
public:
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
using matrix_type = M;
using domain_type = X;
using range_type = Y;
using field_type = typename X::field_type;
using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
#if HAVE_MPI
typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
#else
typedef Dune::CollectiveCommunication< int > communication_type;
using communication_type = Dune::CollectiveCommunication<int>;
#endif
Dune::SolverCategory::Category category() const override
{
return overlapping ?
@ -243,12 +253,12 @@ public:
//! constructor: just store a reference to a matrix
WellModelGhostLastMatrixAdapter (const M& A,
const Opm::LinearOperatorExtra<X, Y>& wellOper,
const LinearOperatorExtra<X, Y>& wellOper,
const std::size_t interiorSize )
: A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
{}
virtual void apply( const X& x, Y& y ) const override
void apply(const X& x, Y& y) const override
{
OPM_TIMEBLOCK(apply);
for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
@ -260,13 +270,13 @@ public:
}
// add well model modification to y
wellOper_.apply(x, y );
wellOper_.apply(x, y);
ghostLastProject( y );
ghostLastProject(y);
}
// y += \alpha * A * x
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
void applyscaleadd (field_type alpha, const X& x, Y& y) const override
{
OPM_TIMEBLOCK(applyscaleadd);
for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
@ -276,23 +286,27 @@ public:
(*col).usmv(alpha, x[col.index()], y[row.index()]);
}
// add scaled well model modification to y
wellOper_.applyscaleadd( alpha, x, y );
wellOper_.applyscaleadd(alpha, x, y);
ghostLastProject( y );
ghostLastProject(y);
}
virtual const matrix_type& getmat() const override { return A_; }
const matrix_type& getmat() const override { return A_; }
void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
void addWellPressureEquations(PressureMatrix& jacobian,
const X& weights,
const bool use_well_weights) const
{
OPM_TIMEBLOCK(addWellPressureEquations);
wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
}
void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
{
OPM_TIMEBLOCK(addWellPressureEquationsStruct);
wellOper_.addWellPressureEquationsStruct(jacobian);
}
int getNumberOfExtraEquations() const
{
return wellOper_.getNumberOfExtraEquations();
@ -307,7 +321,7 @@ protected:
}
const matrix_type& A_ ;
const Opm::LinearOperatorExtra< X, Y>& wellOper_;
const LinearOperatorExtra<X, Y>& wellOper_;
std::size_t interiorSize_;
};

2
opm/simulators/wells/BlackoilWellModelGeneric.hpp
View File

@ -549,7 +549,7 @@ protected:
std::vector<ParallelWellInfo<Scalar>> parallel_well_info_;
std::vector<std::reference_wrapper<ParallelWellInfo<Scalar>>> local_parallel_well_info_;
std::vector<WellProdIndexCalculator> prod_index_calc_;
std::vector<WellProdIndexCalculator<Scalar>> prod_index_calc_;
mutable ParallelWBPCalculation<Scalar> wbpCalculationService_;
std::vector<int> pvt_region_idx_;

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -2274,7 +2274,7 @@ namespace Opm {
auto& ws = this->wellState().well(well->indexOfWell());
for (int p = 0; p < np; ++p) {
// make sure the potentials are positive
ws.well_potentials[p] = std::max(0.0, potentials[p]);
ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
}
}

2
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -135,7 +135,7 @@ namespace Opm {
DeferredLogger& deferred_logger) override;
void updateProductivityIndex(const Simulator& simulator,
const WellProdIndexCalculator& wellPICalc,
const WellProdIndexCalculator<Scalar>& wellPICalc,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const override;

2
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -735,7 +735,7 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
updateProductivityIndex(const Simulator& simulator,
const WellProdIndexCalculator& wellPICalc,
const WellProdIndexCalculator<Scalar>& wellPICalc,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
{

2
opm/simulators/wells/StandardWell.hpp
View File

@ -177,7 +177,7 @@ namespace Opm
DeferredLogger& deferred_logger) override; // should be const?
void updateProductivityIndex(const Simulator& simulator,
const WellProdIndexCalculator& wellPICalc,
const WellProdIndexCalculator<Scalar>& wellPICalc,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const override;

2
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -1220,7 +1220,7 @@ namespace Opm
void
StandardWell<TypeTag>::
updateProductivityIndex(const Simulator& simulator,
const WellProdIndexCalculator& wellPICalc,
const WellProdIndexCalculator<Scalar>& wellPICalc,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const
{

2
opm/simulators/wells/VFPHelpers.cpp
View File

@ -557,7 +557,7 @@ intersectWithIPR(const VFPProdTable& table,
if (y0 < 0 && y1 >= 0){
// crossing with positive slope
Scalar w = -y0/(y1-y0);
w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1~0)
w = std::clamp(w, Scalar{0.0}, Scalar{1.0}); // just to be safe (if y0~y1~0)
flo_x = flo0 + w*(flo1 - flo0);
}
flo0 = flo1;

2
opm/simulators/wells/VFPProdProperties.cpp
View File

@ -202,7 +202,7 @@ bhp(const int table_id,
detail::VFPEvaluation bhp_val = VFPHelpers<Scalar>::interpolate(table, flo_i, thp_i, wfr_i,
gfr_i, alq_i);
bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo);
bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(Scalar{0.0}, bhp_val.dflo) * flo);
bhp.setValue(bhp_val.value);
return bhp;
}

2
opm/simulators/wells/WellFilterCake.cpp
View File

@ -69,7 +69,7 @@ updateFiltrationParticleVolume(const WellInterfaceGeneric<Scalar>& well,
const std::size_t np = well_state.numPhases();
for (int perf = 0; perf < well.numPerfs(); ++perf) {
// not considering the production water
const Scalar water_rates = std::max(0., connection_rates[perf * np + water_index]);
const Scalar water_rates = std::max(Scalar{0.}, connection_rates[perf * np + water_index]);
const Scalar filtrate_rate = water_rates * conc;
filtration_particle_volume_[perf] += filtrate_rate * dt;
ws.perf_data.filtrate_data.rates[perf] = filtrate_rate;

2
opm/simulators/wells/WellGroupHelpers.cpp
View File

@ -802,7 +802,7 @@ updateGpMaintTargetForGroups(const Group& group,
} else {
gpm->resetState(gpmaint_state);
}
group_state.update_gpmaint_target(group.name(), std::max(0.0, sign * rate));
group_state.update_gpmaint_target(group.name(), std::max(Scalar{0.0}, sign * rate));
}
template<class Scalar>

2
opm/simulators/wells/WellInterface.hpp
View File

@ -255,7 +255,7 @@ public:
DeferredLogger& deferred_logger) = 0; // should be const?
virtual void updateProductivityIndex(const Simulator& simulator,
const WellProdIndexCalculator& wellPICalc,
const WellProdIndexCalculator<Scalar>& wellPICalc,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const = 0;

171
opm/simulators/wells/WellProdIndexCalculator.cpp
View File

@ -34,90 +34,102 @@
#include <vector>
namespace {
void checkSizeCompatibility(const Opm::WellProdIndexCalculator& wellPICalc,
const std::vector<double>& connMobility)
{
if (connMobility.size() != wellPICalc.numConnections()) {
throw std::logic_error {
"Mobility vector size does not match expected number of connections"
};
}
}
double logRescale(const double r0, const double rw, const double rd, const double S)
{
const auto numerator = std::log(r0 / rw) + S;
const auto denom = std::log(rd / rw) + S;
return numerator / denom;
}
void standardConnFactorsExplicitDrainRadius(const Opm::Well& well,
std::vector<double>& stdConnFact)
{
const auto& connections = well.getConnections();
const auto rdrain = well.getDrainageRadius();
std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
[rdrain](const Opm::Connection& conn)
{
return conn.CF() * logRescale(conn.r0(), conn.rw(), rdrain, conn.skinFactor());
});
}
void standardConnFactorsDrainIsEquivalent(const Opm::Well& well,
std::vector<double>& stdConnFact)
{
const auto& connections = well.getConnections();
std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
[](const Opm::Connection& conn)
{
return conn.CF();
});
}
std::vector<double> calculateStandardConnFactors(const Opm::Well& well)
{
std::vector<double> stdConnFact(well.getConnections().size());
if (well.getDrainageRadius() > 0.0) {
// Well has an explicit drainage radius. Apply logarithmic
// scaling to the CTFs.
standardConnFactorsExplicitDrainRadius(well, stdConnFact);
}
else {
// Unspecified drainage radius. Standard mobility connection
// scaling factors are just the regular CTFs.
standardConnFactorsDrainIsEquivalent(well, stdConnFact);
}
return stdConnFact;
}
} // namespace Anonymous
Opm::WellProdIndexCalculator::WellProdIndexCalculator(const Well& well)
: standardConnFactors_{ calculateStandardConnFactors(well) }
{}
void Opm::WellProdIndexCalculator::reInit(const Well& well)
template<class Scalar>
void checkSizeCompatibility(const Opm::WellProdIndexCalculator<Scalar>& wellPICalc,
const std::vector<Scalar>& connMobility)
{
this->standardConnFactors_ = calculateStandardConnFactors(well);
if (connMobility.size() != wellPICalc.numConnections()) {
throw std::logic_error {
"Mobility vector size does not match expected number of connections"
};
}
}
double
Opm::WellProdIndexCalculator::
template<class Scalar>
Scalar logRescale(const Scalar r0, const Scalar rw, const Scalar rd, const Scalar S)
{
const auto numerator = std::log(r0 / rw) + S;
const auto denom = std::log(rd / rw) + S;
return numerator / denom;
}
template<class Scalar>
void standardConnFactorsExplicitDrainRadius(const Opm::Well& well,
std::vector<Scalar>& stdConnFact)
{
const auto& connections = well.getConnections();
const auto rdrain = well.getDrainageRadius();
std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
[rdrain](const Opm::Connection& conn)
{
return conn.CF() * logRescale(conn.r0(), conn.rw(), rdrain, conn.skinFactor());
});
}
template<class Scalar>
void standardConnFactorsDrainIsEquivalent(const Opm::Well& well,
std::vector<Scalar>& stdConnFact)
{
const auto& connections = well.getConnections();
std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
[](const Opm::Connection& conn)
{
return conn.CF();
});
}
template<class Scalar>
std::vector<Scalar> calculateStandardConnFactors(const Opm::Well& well)
{
std::vector<Scalar> stdConnFact(well.getConnections().size());
if (well.getDrainageRadius() > 0.0) {
// Well has an explicit drainage radius. Apply logarithmic
// scaling to the CTFs.
standardConnFactorsExplicitDrainRadius(well, stdConnFact);
}
else {
// Unspecified drainage radius. Standard mobility connection
// scaling factors are just the regular CTFs.
standardConnFactorsDrainIsEquivalent(well, stdConnFact);
}
return stdConnFact;
}
} // namespace Anonymous
template<class Scalar>
Opm::WellProdIndexCalculator<Scalar>::
WellProdIndexCalculator(const Well& well)
: standardConnFactors_{ calculateStandardConnFactors<Scalar>(well) }
{}
template<class Scalar>
void Opm::WellProdIndexCalculator<Scalar>::
reInit(const Well& well)
{
this->standardConnFactors_ = calculateStandardConnFactors<Scalar>(well);
}
template<class Scalar>
Scalar Opm::WellProdIndexCalculator<Scalar>::
connectionProdIndStandard(const std::size_t connIdx,
const double connMobility) const
const Scalar connMobility) const
{
return this->standardConnFactors_[connIdx] * connMobility;
}
// ===========================================================================
std::vector<double>
Opm::connectionProdIndStandard(const WellProdIndexCalculator& wellPICalc,
const std::vector<double>& connMobility)
template<class Scalar>
std::vector<Scalar>
Opm::connectionProdIndStandard(const WellProdIndexCalculator<Scalar>& wellPICalc,
const std::vector<Scalar>& connMobility)
{
checkSizeCompatibility(wellPICalc, connMobility);
@ -131,10 +143,19 @@ Opm::connectionProdIndStandard(const WellProdIndexCalculator& wellPICalc,
return connPI;
}
double Opm::wellProdIndStandard(const WellProdIndexCalculator& wellPICalc,
const std::vector<double>& connMobility)
template<class Scalar>
Scalar Opm::wellProdIndStandard(const WellProdIndexCalculator<Scalar>& wellPICalc,
const std::vector<Scalar>& connMobility)
{
const auto connPI = connectionProdIndStandard(wellPICalc, connMobility);
return std::accumulate(connPI.begin(), connPI.end(), 0.0);
}
template class Opm::WellProdIndexCalculator<double>;
template std::vector<double>
Opm::connectionProdIndStandard(const WellProdIndexCalculator<double>&,
const std::vector<double>&);
template double
Opm::wellProdIndStandard(const WellProdIndexCalculator<double>&,
const std::vector<double>&);

156
opm/simulators/wells/WellProdIndexCalculator.hpp
View File

@ -29,87 +29,91 @@ namespace Opm {
namespace Opm {
/// Collect per-connection static information to enable calculating
/// connection-level or well-level productivity index values when
/// incorporating dynamic phase mobilities.
class WellProdIndexCalculator
/// Collect per-connection static information to enable calculating
/// connection-level or well-level productivity index values when
/// incorporating dynamic phase mobilities.
template<class Scalar>
class WellProdIndexCalculator
{
public:
/// Constructor
///
/// \param[in] well Individual well for which to collect
/// per-connection static data.
explicit WellProdIndexCalculator(const Well& well);
/// Reinitialization operation
///
/// Needed to repopulate the internal data members in case of
/// changes to the Well's properties, e.g., as a result of the
/// Well's CTFs being rescaled due to WELPI.
///
/// \param[in] well Individual well for which to collect
/// per-connection static data.
void reInit(const Well& well);
/// Compute connection-level steady-state productivity index value
/// using dynamic phase mobility.
///
/// \param[in] connIdx Linear connection index. Must be in the
/// range 0..numConnections() - 1.
///
/// \param[in] connMobility Phase mobility at connection \p connIdx.
/// Typically derived from dynamic flow state conditions in cell
/// intersected by well's connection \p connIdx.
///
/// \return Connection-level steady-state productivity index.
Scalar connectionProdIndStandard(const std::size_t connIdx,
const Scalar connMobility) const;
/// Number of connections in this well.
///
/// Used primarily for consistency checks.
std::size_t numConnections() const
{
public:
/// Constructor
///
/// \param[in] well Individual well for which to collect
/// per-connection static data.
explicit WellProdIndexCalculator(const Well& well);
return this->standardConnFactors_.size();
}
/// Reinitialization operation
///
/// Needed to repopulate the internal data members in case of
/// changes to the Well's properties, e.g., as a result of the
/// Well's CTFs being rescaled due to WELPI.
///
/// \param[in] well Individual well for which to collect
/// per-connection static data.
void reInit(const Well& well);
private:
/// Static, per-connection multiplicative PI factors.
///
/// Corresponds to the well's connection transmissibility factors,
/// multiplied by a ratio of logarithms if the well has an explicit,
/// positive drainage radius.
std::vector<Scalar> standardConnFactors_{};
};
/// Compute connection-level steady-state productivity index value
/// using dynamic phase mobility.
///
/// \param[in] connIdx Linear connection index. Must be in the
/// range 0..numConnections() - 1.
///
/// \param[in] connMobility Phase mobility at connection \p connIdx.
/// Typically derived from dynamic flow state conditions in cell
/// intersected by well's connection \p connIdx.
///
/// \return Connection-level steady-state productivity index.
double connectionProdIndStandard(const std::size_t connIdx,
const double connMobility) const;
/// Compute connection-level productivity index values for all
/// connections in a well.
///
/// \param[in] wellPICalc Productivity index calculator.
///
/// \param[in] connMobility Phase mobility for each connection.
/// Typically derived from dynamic flow state conditions in cells
/// intersected by well's connections. Must have one value for each
/// \code wellPICalc.numConnections() \endcode well connection.
///
/// \return Connection-level steady-state productivity index values for
/// all connections.
template<class Scalar>
std::vector<Scalar>
connectionProdIndStandard(const WellProdIndexCalculator<Scalar>& wellPICalc,
const std::vector<Scalar>& connMobility);
/// Number of connections in this well.
///
/// Used primarily for consistency checks.
std::size_t numConnections() const
{
return this->standardConnFactors_.size();
}
/// Compute well-level productivity index value.
///
/// \param[in] wellPICalc Productivity index calculator.
///
/// \param[in] connMobility Phase mobility for each connection.
/// Typically derived from dynamic flow state conditions in cells
/// intersected by well's connections. Must have one value for each
/// \code wellPICalc.numConnections() \endcode well connection.
///
/// \return Well-level steady-state productivity index value.
template<class Scalar>
Scalar wellProdIndStandard(const WellProdIndexCalculator<Scalar>& wellPICalc,
const std::vector<Scalar>& connMobility);
private:
/// Static, per-connection multiplicative PI factors.
///
/// Corresponds to the well's connection transmissibility factors,
/// multiplied by a ratio of logarithms if the well has an explicit,
/// positive drainage radius.
std::vector<double> standardConnFactors_{};
};
/// Compute connection-level productivity index values for all
/// connections in a well.
///
/// \param[in] wellPICalc Productivity index calculator.
///
/// \param[in] connMobility Phase mobility for each connection.
/// Typically derived from dynamic flow state conditions in cells
/// intersected by well's connections. Must have one value for each
/// \code wellPICalc.numConnections() \endcode well connection.
///
/// \return Connection-level steady-state productivity index values for
/// all connections.
std::vector<double>
connectionProdIndStandard(const WellProdIndexCalculator& wellPICalc,
const std::vector<double>& connMobility);
/// Compute well-level productivity index value.
///
/// \param[in] wellPICalc Productivity index calculator.
///
/// \param[in] connMobility Phase mobility for each connection.
/// Typically derived from dynamic flow state conditions in cells
/// intersected by well's connections. Must have one value for each
/// \code wellPICalc.numConnections() \endcode well connection.
///
/// \return Well-level steady-state productivity index value.
double wellProdIndStandard(const WellProdIndexCalculator& wellPICalc,
const std::vector<double>& connMobility);
} // namespace Opm
#endif // OPM_WELLPRODINDEXCALCULATOR_HEADER_INCLUDED

64
tests/test_wellprodindexcalculator.cpp
View File

@ -151,7 +151,7 @@ BOOST_AUTO_TEST_SUITE(ConnectionLevel)
BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_SameCF())
};
@ -165,7 +165,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
};
@ -179,7 +179,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin2_SameCF())
};
@ -193,7 +193,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin421_DifferentCF())
};
@ -207,7 +207,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_SameCF())
};
@ -221,7 +221,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
};
@ -235,7 +235,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin2_SameCF())
};
@ -249,7 +249,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin421_DifferentCF())
};
@ -269,7 +269,7 @@ BOOST_AUTO_TEST_SUITE(AllConnections)
BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_SameCF())
};
@ -290,7 +290,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
};
@ -311,7 +311,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin2_SameCF())
};
@ -332,7 +332,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin421_DifferentCF())
};
@ -353,7 +353,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_SameCF())
};
@ -374,7 +374,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
};
@ -395,7 +395,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin2_SameCF())
};
@ -416,7 +416,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin421_DifferentCF())
};
@ -445,7 +445,7 @@ BOOST_AUTO_TEST_SUITE(WellLevel)
BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_SameCF())
};
@ -461,7 +461,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
};
@ -477,7 +477,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin2_SameCF())
};
@ -493,7 +493,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(drainRadDefaulted(), skin421_DifferentCF())
};
@ -509,7 +509,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_SameCF())
};
@ -525,7 +525,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
};
@ -541,7 +541,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin2_SameCF())
};
@ -557,7 +557,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
{
const auto wpiCalc = Opm::WellProdIndexCalculator {
const auto wpiCalc = Opm::WellProdIndexCalculator<double> {
createWell(explicitDrainRad(), skin421_DifferentCF())
};
@ -582,7 +582,7 @@ BOOST_AUTO_TEST_SUITE(Re_Init_Connection_Level)
BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
{
auto well = createWell(drainRadDefaulted(), noSkinFactor_SameCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -605,7 +605,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
{
auto well = createWell(drainRadDefaulted(), noSkinFactor_DifferentCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -628,7 +628,7 @@ BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
{
auto well = createWell(drainRadDefaulted(), skin2_SameCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -651,7 +651,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
{
auto well = createWell(drainRadDefaulted(), skin421_DifferentCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -674,7 +674,7 @@ BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
{
auto well = createWell(explicitDrainRad(), noSkinFactor_SameCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -697,7 +697,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
{
auto well = createWell(explicitDrainRad(), noSkinFactor_DifferentCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -720,7 +720,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
{
auto well = createWell(explicitDrainRad(), skin2_SameCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());
@ -743,7 +743,7 @@ BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
{
auto well = createWell(explicitDrainRad(), skin421_DifferentCF());
auto wpiCalc = Opm::WellProdIndexCalculator { well };
auto wpiCalc = Opm::WellProdIndexCalculator<double> { well };
well.updateWellProductivityIndex( );
const auto scalingFactor = well.convertDeckPI(2.0) / (1.0*liquid_PI_unit());