[cleanup] Remove Eigen from BlackoilModelEbos.

This commit is contained in:
Robert Kloefkorn
2016-10-20 17:47:45 +02:00
parent 2a9a803135
commit c8374a4b95
5 changed files with 127 additions and 45 deletions

View File

@@ -291,6 +291,91 @@ namespace detail {
}
}
template <class Scalar>
inline
double
convergenceReduction(const std::vector< std::vector< Scalar > >& B,
const std::vector< std::vector< Scalar > >& tempV,
const std::vector< std::vector< Scalar > >& R,
std::vector< Scalar >& R_sum,
std::vector< Scalar >& maxCoeff,
std::vector< Scalar >& B_avg,
std::vector< Scalar >& maxNormWell,
const int nc,
const int np,
const std::vector< Scalar >& pv,
const std::vector< Scalar >& residual_well)
{
const int nw = residual_well.size() / np;
assert(nw * np == int(residual_well.size()));
// Do the global reductions
#if 0 // HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
// Compute the global number of cells and porevolume
std::vector<int> v(nc, 1);
auto nc_and_pv = std::tuple<int, double>(0, 0.0);
auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
auto nc_and_pv_containers = std::make_tuple(v, pv);
info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
for ( int idx = 0; idx < np; ++idx )
{
auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
auto containers = std::make_tuple(B.col(idx),
tempV.col(idx),
R.col(idx));
auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
Opm::Reduction::makeGlobalMaxFunctor<double>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
info.computeReduction(containers, operators, values);
B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
maxCoeff[idx] = std::get<1>(values);
R_sum[idx] = std::get<2>(values);
assert(np >= np);
if (idx < np) {
maxNormWell[idx] = 0.0;
for ( int w = 0; w < nw; ++w ) {
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
}
}
}
info.communicator().max(maxNormWell.data(), np);
// Compute pore volume
return std::get<1>(nc_and_pv);
}
else
#endif
{
B_avg.resize(np);
maxCoeff.resize(np);
R_sum.resize(np);
maxNormWell.resize(np);
for ( int idx = 0; idx < np; ++idx )
{
B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(), 0.0 ) / nc;
R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(), 0.0 );
maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() ));
assert(np >= np);
if (idx < np) {
maxNormWell[idx] = 0.0;
for ( int w = 0; w < nw; ++w ) {
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
}
}
}
// Compute total pore volume
return std::accumulate(pv.begin(), pv.end(), 0.0);
}
}
/// \brief Compute the L-infinity norm of a vector representing a well equation.
/// \param a The container to compute the infinity norm on.
/// \param info In a parallel this holds the information about the data distribution.

View File

@@ -51,7 +51,6 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
@@ -104,10 +103,6 @@ namespace Opm {
typedef BlackoilModelEbos ThisType;
public:
// --------- Types and enums ---------
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoil WellState;
typedef BlackoilModelParameters ModelParameters;
@@ -124,8 +119,8 @@ namespace Opm {
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef double Scalar;
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
typedef Dune::FieldVector<Scalar, 3 > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, 3, 3 > MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
//typedef typename SolutionVector :: value_type PrimaryVariables ;
@@ -635,6 +630,8 @@ namespace Opm {
/// \param[in] iteration current iteration number
bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector<double>& residual_norms)
{
typedef std::vector< double > Vector;
const double dt = timer.currentStepLength();
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
@@ -643,50 +640,56 @@ namespace Opm {
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = numPhases();
const V& pv = geo_.poreVolume();
const auto& pv = geo_.poreVolume();
std::vector<double> R_sum(np);
std::vector<double> B_avg(np);
std::vector<double> maxCoeff(np);
std::vector<double> maxNormWell(np);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, np);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, np);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R2(nc, np);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, np);
Vector R_sum(np);
Vector B_avg(np);
Vector maxCoeff(np);
Vector maxNormWell(np);
std::vector< Vector > B( np, Vector( nc ) );
std::vector< Vector > R( np, Vector( nc ) );
std::vector< Vector > R2( np, Vector( nc ) );
std::vector< Vector > tempV( np, Vector( nc ) );
const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
auto ebosResid = ebosSimulator_.model().linearizer().residual();
for ( int idx = 0; idx < np; ++idx )
{
V b(nc);
V r(nc);
Vector& R2_idx = R2[ idx ];
Vector& B_idx = B[ idx ];
const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx);
const int ebosCompIdx = flowPhaseToEbosCompIdx(idx);
for (int cell_idx = 0; cell_idx < nc; ++cell_idx) {
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(idx);
int ebosCompIdx = flowPhaseToEbosCompIdx(idx);
b[cell_idx] = 1 / fs.invB(ebosPhaseIdx).value;
r[cell_idx] = ebosResid[cell_idx][ebosCompIdx];
B_idx [cell_idx] = 1 / fs.invB(ebosPhaseIdx).value;
R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx];
}
R2.col(idx) = r;
B.col(idx) = b;
}
for ( int idx = 0; idx < np; ++idx )
{
tempV.col(idx) = R2.col(idx).abs()/pv;
//tempV.col(idx) = R2.col(idx).abs()/pv;
Vector& tempV_idx = tempV[ idx ];
Vector& R2_idx = R2[ idx ];
for( int cell_idx = 0; cell_idx < nc; ++cell_idx )
{
tempV_idx[ cell_idx ] = std::abs( R2_idx[ cell_idx ] ) / pv[ cell_idx ];
}
}
std::vector<double> pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
Vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
Vector wellResidual = wellModel().residual();
const double pvSum = detail::convergenceReduction(B, tempV, R2,
R_sum, maxCoeff, B_avg, maxNormWell,
nc, np, pv_vector, wellModel().residual());
nc, np, pv_vector, wellResidual );
std::vector<double> CNV(np);
std::vector<double> mass_balance_residual(np);
std::vector<double> well_flux_residual(np);
Vector CNV(np);
Vector mass_balance_residual(np);
Vector well_flux_residual(np);
bool converged_MB = true;
bool converged_CNV = true;
@@ -787,13 +790,6 @@ namespace Opm {
protected:
// --------- Types and enums ---------
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
// --------- Data members ---------
Simulator& ebosSimulator_;

View File

@@ -215,7 +215,8 @@ public:
props_.permeability(),
dynamic_list_econ_limited,
is_parallel_run_,
well_potentials);
well_potentials );
const Wells* wells = wells_manager.c_wells();
WellState well_state;
well_state.init(wells, state, prev_well_state);

View File

@@ -30,7 +30,6 @@
#include <opm/output/vtk/writeVtkData.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BackupRestore.hpp>

View File

@@ -452,7 +452,7 @@ namespace Opm {
/// Diff to bhp for each well perforation.
std::vector<double> wellPerforationPressureDiffs() const
{
{
return well_perforation_pressure_diffs_;
}
@@ -683,9 +683,10 @@ namespace Opm {
const int nw = wells().number_of_wells;
std::vector<double> res(np*nw);
for( int p=0; p<np; ++p) {
const int ebosCompIdx = flowPhaseToEbosCompIdx(p);
for (int i = 0; i < nw; ++i) {
int idx = i + nw*p;
res[idx] = resWell_[i][flowPhaseToEbosCompIdx(p)];
res[idx] = resWell_[ i ][ ebosCompIdx ];
}
}
return res;
@@ -1054,7 +1055,7 @@ namespace Opm {
}
break;
}
}
}
}
}