mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
commit
bbd1395d3f
@ -28,11 +28,13 @@
|
|||||||
#include <opm/core/utility/StopWatch.hpp>
|
#include <opm/core/utility/StopWatch.hpp>
|
||||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||||
|
|
||||||
|
#include <opm/core/utility/platform_dependent/disable_warnings.h>
|
||||||
#if HAVE_SUITESPARSE_UMFPACK_H
|
#if HAVE_SUITESPARSE_UMFPACK_H
|
||||||
#include <Eigen/UmfPackSupport>
|
#include <Eigen/UmfPackSupport>
|
||||||
#else
|
#else
|
||||||
#include <Eigen/IterativeLinearSolvers>
|
#include <Eigen/IterativeLinearSolvers>
|
||||||
#endif
|
#endif
|
||||||
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <cstdlib>
|
#include <cstdlib>
|
||||||
|
@ -53,7 +53,7 @@ namespace Opm {
|
|||||||
/// \param[in] fluid fluid properties
|
/// \param[in] fluid fluid properties
|
||||||
/// \param[in] geo rock properties
|
/// \param[in] geo rock properties
|
||||||
/// \param[in] rock_comp_props if non-null, rock compressibility properties
|
/// \param[in] rock_comp_props if non-null, rock compressibility properties
|
||||||
/// \param[in] wells well structure
|
/// \param[in] wells_arg well structure
|
||||||
/// \param[in] linsolver linear solver
|
/// \param[in] linsolver linear solver
|
||||||
/// \param[in] eclState eclipse state
|
/// \param[in] eclState eclipse state
|
||||||
/// \param[in] has_disgas turn on dissolved gas
|
/// \param[in] has_disgas turn on dissolved gas
|
||||||
@ -64,13 +64,13 @@ namespace Opm {
|
|||||||
const BlackoilPropsAdInterface& fluid,
|
const BlackoilPropsAdInterface& fluid,
|
||||||
const DerivedGeology& geo,
|
const DerivedGeology& geo,
|
||||||
const RockCompressibility* rock_comp_props,
|
const RockCompressibility* rock_comp_props,
|
||||||
const Wells* wells,
|
const Wells* wells_arg,
|
||||||
const NewtonIterationBlackoilInterface& linsolver,
|
const NewtonIterationBlackoilInterface& linsolver,
|
||||||
Opm::EclipseStateConstPtr eclState,
|
Opm::EclipseStateConstPtr eclState,
|
||||||
const bool has_disgas,
|
const bool has_disgas,
|
||||||
const bool has_vapoil,
|
const bool has_vapoil,
|
||||||
const bool terminal_output)
|
const bool terminal_output)
|
||||||
: Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
|
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
|
||||||
eclState, has_disgas, has_vapoil, terminal_output)
|
eclState, has_disgas, has_vapoil, terminal_output)
|
||||||
{
|
{
|
||||||
}
|
}
|
||||||
|
@ -92,6 +92,7 @@ typedef Eigen::Array<double,
|
|||||||
namespace detail {
|
namespace detail {
|
||||||
|
|
||||||
|
|
||||||
|
inline
|
||||||
std::vector<int>
|
std::vector<int>
|
||||||
buildAllCells(const int nc)
|
buildAllCells(const int nc)
|
||||||
{
|
{
|
||||||
@ -146,7 +147,7 @@ namespace detail {
|
|||||||
const BlackoilPropsAdInterface& fluid,
|
const BlackoilPropsAdInterface& fluid,
|
||||||
const DerivedGeology& geo ,
|
const DerivedGeology& geo ,
|
||||||
const RockCompressibility* rock_comp_props,
|
const RockCompressibility* rock_comp_props,
|
||||||
const Wells* wells,
|
const Wells* wells_arg,
|
||||||
const NewtonIterationBlackoilInterface& linsolver,
|
const NewtonIterationBlackoilInterface& linsolver,
|
||||||
Opm::EclipseStateConstPtr eclState,
|
Opm::EclipseStateConstPtr eclState,
|
||||||
const bool has_disgas,
|
const bool has_disgas,
|
||||||
@ -156,7 +157,7 @@ namespace detail {
|
|||||||
, fluid_ (fluid)
|
, fluid_ (fluid)
|
||||||
, geo_ (geo)
|
, geo_ (geo)
|
||||||
, rock_comp_props_(rock_comp_props)
|
, rock_comp_props_(rock_comp_props)
|
||||||
, wells_ (wells)
|
, wells_ (wells_arg)
|
||||||
, vfp_properties_(eclState->getTableManager()->getVFPInjTables(), eclState->getTableManager()->getVFPProdTables())
|
, vfp_properties_(eclState->getTableManager()->getVFPInjTables(), eclState->getTableManager()->getVFPProdTables())
|
||||||
, linsolver_ (linsolver)
|
, linsolver_ (linsolver)
|
||||||
, active_(detail::activePhases(fluid.phaseUsage()))
|
, active_(detail::activePhases(fluid.phaseUsage()))
|
||||||
@ -657,6 +658,7 @@ namespace detail {
|
|||||||
}
|
}
|
||||||
|
|
||||||
namespace detail {
|
namespace detail {
|
||||||
|
inline
|
||||||
double getGravity(const double* g, const int dim) {
|
double getGravity(const double* g, const int dim) {
|
||||||
double grav = 0.0;
|
double grav = 0.0;
|
||||||
if (g) {
|
if (g) {
|
||||||
@ -1103,6 +1105,7 @@ namespace detail {
|
|||||||
|
|
||||||
namespace detail
|
namespace detail
|
||||||
{
|
{
|
||||||
|
inline
|
||||||
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
|
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
|
||||||
const int well,
|
const int well,
|
||||||
const int num_phases,
|
const int num_phases,
|
||||||
@ -1117,6 +1120,7 @@ namespace detail {
|
|||||||
return rate;
|
return rate;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
inline
|
||||||
bool constraintBroken(const std::vector<double>& bhp,
|
bool constraintBroken(const std::vector<double>& bhp,
|
||||||
const std::vector<double>& thp,
|
const std::vector<double>& thp,
|
||||||
const std::vector<double>& well_phase_flow_rate,
|
const std::vector<double>& well_phase_flow_rate,
|
||||||
@ -1194,6 +1198,7 @@ namespace detail {
|
|||||||
* @param well_perforation_densities Densities at well perforations
|
* @param well_perforation_densities Densities at well perforations
|
||||||
* @param gravity Gravitational constant (e.g., 9.81...)
|
* @param gravity Gravitational constant (e.g., 9.81...)
|
||||||
*/
|
*/
|
||||||
|
inline
|
||||||
double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
|
double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
|
||||||
const ADB::V& well_perforation_densities, const double gravity) {
|
const ADB::V& well_perforation_densities, const double gravity) {
|
||||||
const double well_ref_depth = wells.depth_ref[w];
|
const double well_ref_depth = wells.depth_ref[w];
|
||||||
@ -1205,6 +1210,7 @@ namespace detail {
|
|||||||
return dp;
|
return dp;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
inline
|
||||||
ADB::V computeHydrostaticCorrection(const Wells& wells, const ADB::V vfp_ref_depth,
|
ADB::V computeHydrostaticCorrection(const Wells& wells, const ADB::V vfp_ref_depth,
|
||||||
const ADB::V& well_perforation_densities, const double gravity) {
|
const ADB::V& well_perforation_densities, const double gravity) {
|
||||||
const int nw = wells.number_of_wells;
|
const int nw = wells.number_of_wells;
|
||||||
@ -1646,6 +1652,7 @@ namespace detail {
|
|||||||
/// \param a The container to compute the infinity norm on.
|
/// \param a The container to compute the infinity norm on.
|
||||||
/// It has to have one entry for each cell.
|
/// It has to have one entry for each cell.
|
||||||
/// \param info In a parallel this holds the information about the data distribution.
|
/// \param info In a parallel this holds the information about the data distribution.
|
||||||
|
inline
|
||||||
double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() )
|
double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() )
|
||||||
{
|
{
|
||||||
static_cast<void>(pinfo); // Suppress warning in non-MPI case.
|
static_cast<void>(pinfo); // Suppress warning in non-MPI case.
|
||||||
@ -1673,6 +1680,7 @@ namespace detail {
|
|||||||
/// \brief Compute the L-infinity norm of a vector representing a well equation.
|
/// \brief Compute the L-infinity norm of a vector representing a well equation.
|
||||||
/// \param a The container to compute the infinity norm on.
|
/// \param a The container to compute the infinity norm on.
|
||||||
/// \param info In a parallel this holds the information about the data distribution.
|
/// \param info In a parallel this holds the information about the data distribution.
|
||||||
|
inline
|
||||||
double infinityNormWell( const ADB& a, const boost::any& pinfo )
|
double infinityNormWell( const ADB& a, const boost::any& pinfo )
|
||||||
{
|
{
|
||||||
static_cast<void>(pinfo); // Suppress warning in non-MPI case.
|
static_cast<void>(pinfo); // Suppress warning in non-MPI case.
|
||||||
|
@ -698,29 +698,29 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
const ADB& sg,
|
const ADB& sg,
|
||||||
const Cells& cells) const
|
const Cells& cells) const
|
||||||
{
|
{
|
||||||
const int numCells = cells.size();
|
const int nCells = cells.size();
|
||||||
const int numActivePhases = numPhases();
|
const int nActivePhases = numPhases();
|
||||||
const int numBlocks = so.numBlocks();
|
const int nBlocks = so.numBlocks();
|
||||||
|
|
||||||
Block activeSat(numCells, numActivePhases);
|
Block activeSat(nCells, nActivePhases);
|
||||||
if (phase_usage_.phase_used[Water]) {
|
if (phase_usage_.phase_used[Water]) {
|
||||||
assert(sw.value().size() == numCells);
|
assert(sw.value().size() == nCells);
|
||||||
activeSat.col(phase_usage_.phase_pos[Water]) = sw.value();
|
activeSat.col(phase_usage_.phase_pos[Water]) = sw.value();
|
||||||
}
|
}
|
||||||
if (phase_usage_.phase_used[Oil]) {
|
if (phase_usage_.phase_used[Oil]) {
|
||||||
assert(so.value().size() == numCells);
|
assert(so.value().size() == nCells);
|
||||||
activeSat.col(phase_usage_.phase_pos[Oil]) = so.value();
|
activeSat.col(phase_usage_.phase_pos[Oil]) = so.value();
|
||||||
} else {
|
} else {
|
||||||
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
||||||
}
|
}
|
||||||
if (phase_usage_.phase_used[Gas]) {
|
if (phase_usage_.phase_used[Gas]) {
|
||||||
assert(sg.value().size() == numCells);
|
assert(sg.value().size() == nCells);
|
||||||
activeSat.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
activeSat.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
||||||
}
|
}
|
||||||
|
|
||||||
Block pc(numCells, numActivePhases);
|
Block pc(nCells, nActivePhases);
|
||||||
Block dpc(numCells, numActivePhases*numActivePhases);
|
Block dpc(nCells, nActivePhases*nActivePhases);
|
||||||
satprops_->capPress(numCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
|
satprops_->capPress(nCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
|
||||||
|
|
||||||
std::vector<ADB> adbCapPressures;
|
std::vector<ADB> adbCapPressures;
|
||||||
adbCapPressures.reserve(3);
|
adbCapPressures.reserve(3);
|
||||||
@ -728,18 +728,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
||||||
if (phase_usage_.phase_used[phase1]) {
|
if (phase_usage_.phase_used[phase1]) {
|
||||||
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
||||||
std::vector<ADB::M> jacs(numBlocks);
|
std::vector<ADB::M> jacs(nBlocks);
|
||||||
for (int block = 0; block < numBlocks; ++block) {
|
for (int block = 0; block < nBlocks; ++block) {
|
||||||
jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
|
jacs[block] = ADB::M(nCells, s[phase1]->derivative()[block].cols());
|
||||||
}
|
}
|
||||||
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
||||||
if (!phase_usage_.phase_used[phase2])
|
if (!phase_usage_.phase_used[phase2])
|
||||||
continue;
|
continue;
|
||||||
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
||||||
// Assemble dpc1/ds2.
|
// Assemble dpc1/ds2.
|
||||||
const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
const int column = phase1_pos + nActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
||||||
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
|
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
|
||||||
for (int block = 0; block < numBlocks; ++block) {
|
for (int block = 0; block < nBlocks; ++block) {
|
||||||
ADB::M temp;
|
ADB::M temp;
|
||||||
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp);
|
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp);
|
||||||
jacs[block] += temp;
|
jacs[block] += temp;
|
||||||
|
@ -76,14 +76,14 @@ namespace Opm {
|
|||||||
const DerivedGeology& geo,
|
const DerivedGeology& geo,
|
||||||
const RockCompressibility* rock_comp_props,
|
const RockCompressibility* rock_comp_props,
|
||||||
const SolventPropsAdFromDeck& solvent_props,
|
const SolventPropsAdFromDeck& solvent_props,
|
||||||
const Wells* wells,
|
const Wells* wells_arg,
|
||||||
const NewtonIterationBlackoilInterface& linsolver,
|
const NewtonIterationBlackoilInterface& linsolver,
|
||||||
const EclipseStateConstPtr eclState,
|
const EclipseStateConstPtr eclState,
|
||||||
const bool has_disgas,
|
const bool has_disgas,
|
||||||
const bool has_vapoil,
|
const bool has_vapoil,
|
||||||
const bool terminal_output,
|
const bool terminal_output,
|
||||||
const bool has_solvent)
|
const bool has_solvent)
|
||||||
: Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
|
: Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver,
|
||||||
eclState, has_disgas, has_vapoil, terminal_output),
|
eclState, has_disgas, has_vapoil, terminal_output),
|
||||||
has_solvent_(has_solvent),
|
has_solvent_(has_solvent),
|
||||||
solvent_pos_(detail::solventPos(fluid.phaseUsage())),
|
solvent_pos_(detail::solventPos(fluid.phaseUsage())),
|
||||||
|
@ -31,11 +31,13 @@
|
|||||||
#include <opm/core/utility/Exceptions.hpp>
|
#include <opm/core/utility/Exceptions.hpp>
|
||||||
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
||||||
|
|
||||||
|
#include <opm/core/utility/platform_dependent/disable_warnings.h>
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
#include <Eigen/UmfPackSupport>
|
#include <Eigen/UmfPackSupport>
|
||||||
#else
|
#else
|
||||||
#include <Eigen/SparseLU>
|
#include <Eigen/SparseLU>
|
||||||
#endif
|
#endif
|
||||||
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
@ -52,10 +54,10 @@ namespace Opm
|
|||||||
|
|
||||||
/// Construct a system solver.
|
/// Construct a system solver.
|
||||||
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
|
NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param,
|
||||||
const boost::any& parallelInformation)
|
const boost::any& parallelInformation_arg)
|
||||||
: cpr_param_( param ),
|
: cpr_param_( param ),
|
||||||
iterations_( 0 ),
|
iterations_( 0 ),
|
||||||
parallelInformation_(parallelInformation),
|
parallelInformation_(parallelInformation_arg),
|
||||||
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
||||||
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
||||||
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
||||||
|
@ -87,18 +87,18 @@ namespace Opm
|
|||||||
template<int category=Dune::SolverCategory::sequential, class O, class P>
|
template<int category=Dune::SolverCategory::sequential, class O, class P>
|
||||||
void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
|
void constructPreconditionerAndSolve(O& opA, DuneMatrix& istlAe,
|
||||||
Vector& x, Vector& istlb,
|
Vector& x, Vector& istlb,
|
||||||
const P& parallelInformation,
|
const P& parallelInformation_arg,
|
||||||
const P& parallelInformationAe,
|
const P& parallelInformationAe,
|
||||||
Dune::InverseOperatorResult& result) const
|
Dune::InverseOperatorResult& result) const
|
||||||
{
|
{
|
||||||
typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
|
typedef Dune::ScalarProductChooser<Vector,P,category> ScalarProductChooser;
|
||||||
std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
|
std::unique_ptr<typename ScalarProductChooser::ScalarProduct>
|
||||||
sp(ScalarProductChooser::construct(parallelInformation));
|
sp(ScalarProductChooser::construct(parallelInformation_arg));
|
||||||
// Construct preconditioner.
|
// Construct preconditioner.
|
||||||
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
|
// typedef Dune::SeqILU0<Mat,Vector,Vector> Preconditioner;
|
||||||
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
|
typedef Opm::CPRPreconditioner<Mat,Vector,Vector,P> Preconditioner;
|
||||||
parallelInformation.copyOwnerToAll(istlb, istlb);
|
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||||
Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation,
|
Preconditioner precond(cpr_param_, opA.getmat(), istlAe, parallelInformation_arg,
|
||||||
parallelInformationAe);
|
parallelInformationAe);
|
||||||
|
|
||||||
// TODO: Revise when linear solvers interface opm-core is done
|
// TODO: Revise when linear solvers interface opm-core is done
|
||||||
|
@ -30,11 +30,13 @@
|
|||||||
#include <opm/core/utility/Exceptions.hpp>
|
#include <opm/core/utility/Exceptions.hpp>
|
||||||
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
#include <opm/core/linalg/ParallelIstlInformation.hpp>
|
||||||
|
|
||||||
|
#include <opm/core/utility/platform_dependent/disable_warnings.h>
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
#include <Eigen/UmfPackSupport>
|
#include <Eigen/UmfPackSupport>
|
||||||
#else
|
#else
|
||||||
#include <Eigen/SparseLU>
|
#include <Eigen/SparseLU>
|
||||||
#endif
|
#endif
|
||||||
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
@ -51,9 +53,9 @@ namespace Opm
|
|||||||
|
|
||||||
/// Construct a system solver.
|
/// Construct a system solver.
|
||||||
NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
|
NewtonIterationBlackoilInterleaved::NewtonIterationBlackoilInterleaved(const parameter::ParameterGroup& param,
|
||||||
const boost::any& parallelInformation)
|
const boost::any& parallelInformation_arg)
|
||||||
: iterations_( 0 ),
|
: iterations_( 0 ),
|
||||||
parallelInformation_(parallelInformation),
|
parallelInformation_(parallelInformation_arg),
|
||||||
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
newton_use_gmres_( param.getDefault("newton_use_gmres", false ) ),
|
||||||
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
linear_solver_reduction_( param.getDefault("linear_solver_reduction", 1e-2 ) ),
|
||||||
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
linear_solver_maxiter_( param.getDefault("linear_solver_maxiter", 50 ) ),
|
||||||
|
@ -92,19 +92,19 @@ namespace Opm
|
|||||||
template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
|
template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
|
||||||
void constructPreconditionerAndSolve(O& opA,
|
void constructPreconditionerAndSolve(O& opA,
|
||||||
Vector& x, Vector& istlb,
|
Vector& x, Vector& istlb,
|
||||||
const POrComm& parallelInformation,
|
const POrComm& parallelInformation_arg,
|
||||||
Dune::InverseOperatorResult& result) const
|
Dune::InverseOperatorResult& result) const
|
||||||
{
|
{
|
||||||
// Construct scalar product.
|
// Construct scalar product.
|
||||||
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
|
||||||
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
|
||||||
SPPointer sp(ScalarProductChooser::construct(parallelInformation));
|
SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
|
||||||
|
|
||||||
// Construct preconditioner.
|
// Construct preconditioner.
|
||||||
auto precond = constructPrecond(opA, parallelInformation);
|
auto precond = constructPrecond(opA, parallelInformation_arg);
|
||||||
|
|
||||||
// Communicate if parallel.
|
// Communicate if parallel.
|
||||||
parallelInformation.copyOwnerToAll(istlb, istlb);
|
parallelInformation_arg.copyOwnerToAll(istlb, istlb);
|
||||||
|
|
||||||
// Solve.
|
// Solve.
|
||||||
solve(opA, x, istlb, *sp, precond, result);
|
solve(opA, x, istlb, *sp, precond, result);
|
||||||
|
@ -33,8 +33,8 @@ namespace Opm
|
|||||||
/// \param[in] parallelInformation In the case of a parallel run
|
/// \param[in] parallelInformation In the case of a parallel run
|
||||||
/// with dune-istl the information about the parallelization.
|
/// with dune-istl the information about the parallelization.
|
||||||
NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
|
NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param,
|
||||||
const boost::any& parallelInformation)
|
const boost::any& parallelInformation_arg)
|
||||||
: iterations_( 0 ), parallelInformation_(parallelInformation)
|
: iterations_( 0 ), parallelInformation_(parallelInformation_arg)
|
||||||
{
|
{
|
||||||
linsolver_.reset(new LinearSolverFactory(param));
|
linsolver_.reset(new LinearSolverFactory(param));
|
||||||
}
|
}
|
||||||
|
@ -24,11 +24,13 @@
|
|||||||
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
||||||
#include <opm/core/utility/ErrorMacros.hpp>
|
#include <opm/core/utility/ErrorMacros.hpp>
|
||||||
|
|
||||||
|
#include <opm/core/utility/platform_dependent/disable_warnings.h>
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
#include <Eigen/UmfPackSupport>
|
#include <Eigen/UmfPackSupport>
|
||||||
#else
|
#else
|
||||||
#include <Eigen/SparseLU>
|
#include <Eigen/SparseLU>
|
||||||
#endif
|
#endif
|
||||||
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
|
@ -77,7 +77,7 @@ namespace Opm
|
|||||||
bool isOscillate = false;
|
bool isOscillate = false;
|
||||||
bool isStagnate = false;
|
bool isStagnate = false;
|
||||||
const enum RelaxType relaxtype = relaxType();
|
const enum RelaxType relaxtype = relaxType();
|
||||||
int linearIterations = 0;
|
int linIters = 0;
|
||||||
|
|
||||||
// ---------- Main Newton loop ----------
|
// ---------- Main Newton loop ----------
|
||||||
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
|
while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) {
|
||||||
@ -85,7 +85,7 @@ namespace Opm
|
|||||||
V dx = model_->solveJacobianSystem();
|
V dx = model_->solveJacobianSystem();
|
||||||
|
|
||||||
// Store number of linear iterations used.
|
// Store number of linear iterations used.
|
||||||
linearIterations += model_->linearIterationsLastSolve();
|
linIters += model_->linearIterationsLastSolve();
|
||||||
|
|
||||||
// Stabilize the Newton update.
|
// Stabilize the Newton update.
|
||||||
detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
|
detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate);
|
||||||
@ -119,15 +119,15 @@ namespace Opm
|
|||||||
return -1; // -1 indicates that the solver has to be restarted
|
return -1; // -1 indicates that the solver has to be restarted
|
||||||
}
|
}
|
||||||
|
|
||||||
linearIterations_ += linearIterations;
|
linearIterations_ += linIters;
|
||||||
newtonIterations_ += iteration;
|
newtonIterations_ += iteration;
|
||||||
linearIterationsLast_ = linearIterations;
|
linearIterationsLast_ = linIters;
|
||||||
newtonIterationsLast_ = iteration;
|
newtonIterationsLast_ = iteration;
|
||||||
|
|
||||||
// Do model-specific post-step actions.
|
// Do model-specific post-step actions.
|
||||||
model_->afterStep(dt, reservoir_state, well_state);
|
model_->afterStep(dt, reservoir_state, well_state);
|
||||||
|
|
||||||
return linearIterations;
|
return linIters;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -178,7 +178,7 @@ namespace Opm
|
|||||||
template <class PhysicalModel>
|
template <class PhysicalModel>
|
||||||
void
|
void
|
||||||
NewtonSolver<PhysicalModel>::detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
|
NewtonSolver<PhysicalModel>::detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
|
||||||
const int it, const double relaxRelTol,
|
const int it, const double relaxRelTol_arg,
|
||||||
bool& oscillate, bool& stagnate) const
|
bool& oscillate, bool& stagnate) const
|
||||||
{
|
{
|
||||||
// The detection of oscillation in two primary variable results in the report of the detection
|
// The detection of oscillation in two primary variable results in the report of the detection
|
||||||
@ -201,7 +201,7 @@ namespace Opm
|
|||||||
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
|
const double d1 = std::abs((F0[p] - F2[p]) / F0[p]);
|
||||||
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
|
const double d2 = std::abs((F0[p] - F1[p]) / F0[p]);
|
||||||
|
|
||||||
oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2);
|
oscillatePhase += (d1 < relaxRelTol_arg) && (relaxRelTol_arg < d2);
|
||||||
|
|
||||||
// Process is 'stagnate' unless at least one phase
|
// Process is 'stagnate' unless at least one phase
|
||||||
// exhibits significant residual change.
|
// exhibits significant residual change.
|
||||||
|
@ -180,11 +180,11 @@ public:
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
template<class B, class T>
|
template<class B, class T>
|
||||||
void scatter(B& buffer, const T& e, std::size_t size)
|
void scatter(B& buffer, const T& e, std::size_t size_arg)
|
||||||
{
|
{
|
||||||
assert( T::codimension == 0);
|
assert( T::codimension == 0);
|
||||||
assert( int(size) == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
|
assert( int(size_arg) == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
|
||||||
static_cast<void>(size);
|
static_cast<void>(size_arg);
|
||||||
|
|
||||||
double val;
|
double val;
|
||||||
for ( int i=0; i<recvState_.numPhases(); ++i )
|
for ( int i=0; i<recvState_.numPhases(); ++i )
|
||||||
@ -283,10 +283,10 @@ public:
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
template<class B, class T>
|
template<class B, class T>
|
||||||
void scatter(B& buffer, const T& e, std::size_t size)
|
void scatter(B& buffer, const T& e, std::size_t size_arg)
|
||||||
{
|
{
|
||||||
assert( T::codimension == 0);
|
assert( T::codimension == 0);
|
||||||
assert( size==size_ ); (void) size;
|
assert( size_arg==size_ ); (void) size_arg;
|
||||||
double val;
|
double val;
|
||||||
buffer.read(val);
|
buffer.read(val);
|
||||||
recvProps_.cellPvtRegionIdx_[e.index()]=val;
|
recvProps_.cellPvtRegionIdx_[e.index()]=val;
|
||||||
|
@ -910,6 +910,7 @@ inline double findTHP(
|
|||||||
}
|
}
|
||||||
//Canary in a coal mine: shouldn't really be required
|
//Canary in a coal mine: shouldn't really be required
|
||||||
assert(found == true);
|
assert(found == true);
|
||||||
|
static_cast<void>(found); //Silence compiler warning
|
||||||
|
|
||||||
const double& x0 = thp_array[i ];
|
const double& x0 = thp_array[i ];
|
||||||
const double& x1 = thp_array[i+1];
|
const double& x1 = thp_array[i+1];
|
||||||
|
@ -74,7 +74,7 @@ VFPInjProperties::VFPInjProperties(const std::map<int, VFPInjTable>& tables) {
|
|||||||
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
||||||
const Wells& wells,
|
const Wells& wells,
|
||||||
const ADB& qs,
|
const ADB& qs,
|
||||||
const ADB& thp) const {
|
const ADB& thp_val) const {
|
||||||
const int nw = wells.number_of_wells;
|
const int nw = wells.number_of_wells;
|
||||||
|
|
||||||
//Short-hands for water / oil / gas phases
|
//Short-hands for water / oil / gas phases
|
||||||
@ -84,7 +84,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
|||||||
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
|
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
|
||||||
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
|
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
|
||||||
|
|
||||||
return bhp(table_id, w, o, g, thp);
|
return bhp(table_id, w, o, g, thp_val);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -97,16 +97,16 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
|||||||
const ADB& aqua,
|
const ADB& aqua,
|
||||||
const ADB& liquid,
|
const ADB& liquid,
|
||||||
const ADB& vapour,
|
const ADB& vapour,
|
||||||
const ADB& thp) const {
|
const ADB& thp_arg) const {
|
||||||
const int nw = thp.size();
|
const int nw = thp_arg.size();
|
||||||
|
|
||||||
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp);
|
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg);
|
||||||
|
|
||||||
assert(static_cast<int>(table_id.size()) == nw);
|
assert(static_cast<int>(table_id.size()) == nw);
|
||||||
assert(aqua.size() == nw);
|
assert(aqua.size() == nw);
|
||||||
assert(liquid.size() == nw);
|
assert(liquid.size() == nw);
|
||||||
assert(vapour.size() == nw);
|
assert(vapour.size() == nw);
|
||||||
assert(thp.size() == nw);
|
assert(thp_arg.size() == nw);
|
||||||
|
|
||||||
//Allocate data for bhp's and partial derivatives
|
//Allocate data for bhp's and partial derivatives
|
||||||
ADB::V value = ADB::V::Zero(nw);
|
ADB::V value = ADB::V::Zero(nw);
|
||||||
@ -130,7 +130,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
|||||||
if (table != nullptr) {
|
if (table != nullptr) {
|
||||||
//First, find the values to interpolate between
|
//First, find the values to interpolate between
|
||||||
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
|
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
|
||||||
auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis());
|
auto thp_i = detail::findInterpData(thp_arg.value()[i], table->getTHPAxis());
|
||||||
|
|
||||||
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
|
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
|
||||||
|
|
||||||
@ -155,8 +155,8 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
|
|||||||
//but may not save too much on that.
|
//but may not save too much on that.
|
||||||
jacs[block] = ADB::M(nw, block_pattern[block]);
|
jacs[block] = ADB::M(nw, block_pattern[block]);
|
||||||
|
|
||||||
if (!thp.derivative().empty()) {
|
if (!thp_arg.derivative().empty()) {
|
||||||
jacs[block] += dthp_diag * thp.derivative()[block];
|
jacs[block] += dthp_diag * thp_arg.derivative()[block];
|
||||||
}
|
}
|
||||||
if (!flo.derivative().empty()) {
|
if (!flo.derivative().empty()) {
|
||||||
jacs[block] += dflo_diag * flo.derivative()[block];
|
jacs[block] += dflo_diag * flo.derivative()[block];
|
||||||
@ -175,10 +175,10 @@ double VFPInjProperties::bhp(int table_id,
|
|||||||
const double& aqua,
|
const double& aqua,
|
||||||
const double& liquid,
|
const double& liquid,
|
||||||
const double& vapour,
|
const double& vapour,
|
||||||
const double& thp) const {
|
const double& thp_arg) const {
|
||||||
const VFPInjTable* table = detail::getTable(m_tables, table_id);
|
const VFPInjTable* table = detail::getTable(m_tables, table_id);
|
||||||
|
|
||||||
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp);
|
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg);
|
||||||
return retval.value;
|
return retval.value;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -193,7 +193,7 @@ double VFPInjProperties::thp(int table_id,
|
|||||||
const double& aqua,
|
const double& aqua,
|
||||||
const double& liquid,
|
const double& liquid,
|
||||||
const double& vapour,
|
const double& vapour,
|
||||||
const double& bhp) const {
|
const double& bhp_arg) const {
|
||||||
const VFPInjTable* table = detail::getTable(m_tables, table_id);
|
const VFPInjTable* table = detail::getTable(m_tables, table_id);
|
||||||
const VFPInjTable::array_type& data = table->getTable();
|
const VFPInjTable::array_type& data = table->getTable();
|
||||||
|
|
||||||
@ -215,8 +215,8 @@ double VFPInjProperties::thp(int table_id,
|
|||||||
bhp_array[i] = detail::interpolate(data, flo_i, thp_i).value;
|
bhp_array[i] = detail::interpolate(data, flo_i, thp_i).value;
|
||||||
}
|
}
|
||||||
|
|
||||||
double thp = detail::findTHP(bhp_array, thp_array, bhp);
|
double retval = detail::findTHP(bhp_array, thp_array, bhp_arg);
|
||||||
return thp;
|
return retval;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
@ -59,7 +59,7 @@ VFPProdProperties::VFPProdProperties(const std::map<int, VFPProdTable>& tables)
|
|||||||
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
||||||
const Wells& wells,
|
const Wells& wells,
|
||||||
const ADB& qs,
|
const ADB& qs,
|
||||||
const ADB& thp,
|
const ADB& thp_arg,
|
||||||
const ADB& alq) const {
|
const ADB& alq) const {
|
||||||
const int nw = wells.number_of_wells;
|
const int nw = wells.number_of_wells;
|
||||||
|
|
||||||
@ -70,7 +70,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
|||||||
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
|
const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
|
||||||
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
|
const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
|
||||||
|
|
||||||
return bhp(table_id, w, o, g, thp, alq);
|
return bhp(table_id, w, o, g, thp_arg, alq);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
@ -82,17 +82,17 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
|||||||
const ADB& aqua,
|
const ADB& aqua,
|
||||||
const ADB& liquid,
|
const ADB& liquid,
|
||||||
const ADB& vapour,
|
const ADB& vapour,
|
||||||
const ADB& thp,
|
const ADB& thp_arg,
|
||||||
const ADB& alq) const {
|
const ADB& alq) const {
|
||||||
const int nw = thp.size();
|
const int nw = thp_arg.size();
|
||||||
|
|
||||||
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp, alq);
|
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg, alq);
|
||||||
|
|
||||||
assert(static_cast<int>(table_id.size()) == nw);
|
assert(static_cast<int>(table_id.size()) == nw);
|
||||||
assert(aqua.size() == nw);
|
assert(aqua.size() == nw);
|
||||||
assert(liquid.size() == nw);
|
assert(liquid.size() == nw);
|
||||||
assert(vapour.size() == nw);
|
assert(vapour.size() == nw);
|
||||||
assert(thp.size() == nw);
|
assert(thp_arg.size() == nw);
|
||||||
assert(alq.size() == nw);
|
assert(alq.size() == nw);
|
||||||
|
|
||||||
//Allocate data for bhp's and partial derivatives
|
//Allocate data for bhp's and partial derivatives
|
||||||
@ -123,7 +123,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
|||||||
//First, find the values to interpolate between
|
//First, find the values to interpolate between
|
||||||
//Value of FLO is negative in OPM for producers, but positive in VFP table
|
//Value of FLO is negative in OPM for producers, but positive in VFP table
|
||||||
auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
|
auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis());
|
||||||
auto thp_i = detail::findInterpData( thp.value()[i], table->getTHPAxis());
|
auto thp_i = detail::findInterpData( thp_arg.value()[i], table->getTHPAxis());
|
||||||
auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
|
auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis());
|
||||||
auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
|
auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis());
|
||||||
auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis());
|
auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis());
|
||||||
@ -157,8 +157,8 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
|
|||||||
//but may not save too much on that.
|
//but may not save too much on that.
|
||||||
jacs[block] = ADB::M(nw, block_pattern[block]);
|
jacs[block] = ADB::M(nw, block_pattern[block]);
|
||||||
|
|
||||||
if (!thp.derivative().empty()) {
|
if (!thp_arg.derivative().empty()) {
|
||||||
jacs[block] += dthp_diag * thp.derivative()[block];
|
jacs[block] += dthp_diag * thp_arg.derivative()[block];
|
||||||
}
|
}
|
||||||
if (!wfr.derivative().empty()) {
|
if (!wfr.derivative().empty()) {
|
||||||
jacs[block] += dwfr_diag * wfr.derivative()[block];
|
jacs[block] += dwfr_diag * wfr.derivative()[block];
|
||||||
@ -184,11 +184,11 @@ double VFPProdProperties::bhp(int table_id,
|
|||||||
const double& aqua,
|
const double& aqua,
|
||||||
const double& liquid,
|
const double& liquid,
|
||||||
const double& vapour,
|
const double& vapour,
|
||||||
const double& thp,
|
const double& thp_arg,
|
||||||
const double& alq) const {
|
const double& alq) const {
|
||||||
const VFPProdTable* table = detail::getTable(m_tables, table_id);
|
const VFPProdTable* table = detail::getTable(m_tables, table_id);
|
||||||
|
|
||||||
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp, alq);
|
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq);
|
||||||
return retval.value;
|
return retval.value;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -198,7 +198,7 @@ double VFPProdProperties::thp(int table_id,
|
|||||||
const double& aqua,
|
const double& aqua,
|
||||||
const double& liquid,
|
const double& liquid,
|
||||||
const double& vapour,
|
const double& vapour,
|
||||||
const double& bhp,
|
const double& bhp_arg,
|
||||||
const double& alq) const {
|
const double& alq) const {
|
||||||
const VFPProdTable* table = detail::getTable(m_tables, table_id);
|
const VFPProdTable* table = detail::getTable(m_tables, table_id);
|
||||||
const VFPProdTable::array_type& data = table->getTable();
|
const VFPProdTable::array_type& data = table->getTable();
|
||||||
@ -227,8 +227,8 @@ double VFPProdProperties::thp(int table_id,
|
|||||||
bhp_array[i] = detail::interpolate(data, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
|
bhp_array[i] = detail::interpolate(data, flo_i, thp_i, wfr_i, gfr_i, alq_i).value;
|
||||||
}
|
}
|
||||||
|
|
||||||
double thp = detail::findTHP(bhp_array, thp_array, bhp);
|
double retval = detail::findTHP(bhp_array, thp_array, bhp_arg);
|
||||||
return thp;
|
return retval;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user