diff --git a/opm/autodiff/BlackoilDetails.hpp b/opm/autodiff/BlackoilDetails.hpp index 7f72af998..9763e0817 100644 --- a/opm/autodiff/BlackoilDetails.hpp +++ b/opm/autodiff/BlackoilDetails.hpp @@ -220,10 +220,10 @@ namespace detail { int nc, int np, const std::vector pv, - LinearisedBlackoilResidual residual) + std::vector residual_well) { - const int nw = residual.well_flux_eq.size() / np; - assert(nw * np == int(residual.well_flux_eq.size())); + const int nw = residual_well.size() / np; + assert(nw * np == int(residual_well.size())); // Do the global reductions #if HAVE_MPI @@ -257,7 +257,7 @@ namespace detail { if (idx < np) { maxNormWell[idx] = 0.0; for ( int w = 0; w < nw; ++w ) { - maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w])); + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); } } } @@ -282,7 +282,7 @@ namespace detail { if (idx < np) { maxNormWell[idx] = 0.0; for ( int w = 0; w < nw; ++w ) { - maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual.well_flux_eq.value()[nw*idx + w])); + maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w])); } } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c4a4b2ee5..106724468 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -116,6 +116,12 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; + + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; //typedef typename SolutionVector :: value_type PrimaryVariables ; // --------- Public methods --------- @@ -175,9 +181,7 @@ namespace Opm { { const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); - well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth); - wellModel().setWellsActive( localWellsActive() ); global_nc_ = Opm::AutoDiffGrid::numCells(grid_); } @@ -229,7 +233,7 @@ namespace Opm { const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged); if (must_solve) { // enable single precision for solvers when dt is smaller then 20 days - residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ; + //residual_.singlePrecision = (unit::convert::to(dt, unit::day) < 20.) ; // Compute the nonlinear update. V dx = solveJacobianSystem(); @@ -283,10 +287,22 @@ namespace Opm { // -------- Mass balance equations -------- assembleMassBalanceEq(timer, iterationIdx, reservoir_state); - // -------- Well equations ---------- double dt = timer.currentStepLength(); + IterationReport iter_report = wellModel().assemble(ebosSimulator_, iterationIdx, dt, well_state, residual_); + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + + typedef Dune::MatrixAdapter Operator; + auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + convertResults2(ebosResid, ebosJac); + wellModel().addRhs(ebosResid, ebosJac); return iter_report; } @@ -392,7 +408,82 @@ namespace Opm { /// r is the residual. V solveJacobianSystem() const { - return linsolver_.computeNewtonIncrement(residual_); + + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + + typedef Dune::MatrixAdapter Operator; + auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator_.model().linearizer().residual(); + + Operator opA(ebosJac); + const double relax = 0.9; + typedef Dune::SeqILU0 SeqPreconditioner; + SeqPreconditioner precond(opA.getmat(), relax); + std::cout << "hei" << std::endl; + + Dune::SeqScalarProduct sp; + + + wellModel().apply(ebosJac, ebosResid); + Dune::BiCGSTABSolver linsolve(opA, sp, precond, + 0.001, + 100, + false); + + + const int np = numPhases(); + const int nc = AutoDiffGrid::numCells(grid_); +// std::cout << "ebosResid" << std::endl; +// for( int p=0; p maxNormWell(np); Eigen::Array B(nc, np); Eigen::Array R(nc, np); + Eigen::Array R2(nc, np); Eigen::Array tempV(nc, np); + auto ebosResid = ebosSimulator_.model().linearizer().residual(); for ( int idx = 0; idx < np; ++idx ) { V b(nc); + V r(nc); 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]; + } + R2.col(idx) = r; B.col(idx) = b; } - for ( int idx = 0; idx < np; ++idx ) { //const ADB& tempB = rq_[idx].b; //B.col(idx) = 1./tempB.value(); R.col(idx) = residual_.material_balance_eq[idx].value(); tempV.col(idx) = R.col(idx).abs()/pv; + std::cout << "------R------- " << idx << std::endl; + std::cout << R.col(idx)[0] << std::endl; + std::cout << "------R2------- " << idx << std::endl; + std::cout << R2.col(idx)[0] << std::endl; } std::vector pv_vector (geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); - const double pvSum = detail::convergenceReduction(B, tempV, R, + const double pvSum = detail::convergenceReduction(B, tempV, R2, R_sum, maxCoeff, B_avg, maxNormWell, - nc, np, pv_vector, residual_); + nc, np, pv_vector, wellModel().residual()); std::vector CNV(np); std::vector mass_balance_residual(np); @@ -1006,13 +1107,11 @@ namespace Opm { } public: - /* int ebosCompToFlowPhaseIdx( const int compIdx ) const { const int compToPhase[ 3 ] = { Oil, Water, Gas }; return compToPhase[ compIdx ]; } - */ int flowToEbosPvIdx( const int flowPv ) const { @@ -1032,6 +1131,7 @@ namespace Opm { + private: void convertResults(const Simulator& simulator) { @@ -1136,6 +1236,52 @@ namespace Opm { std::move(adbJacs[eqIdx])); } } + void convertResults2(BVector& ebosResid, Mat& ebosJac) const + { + const int numPhases = wells().number_of_phases; + const int numCells = ebosJac.N(); + const int cols = ebosJac.M(); + assert( numCells == cols ); + + // write the right-hand-side values from the ebosJac into the objects + // allocated above. + const auto endrow = ebosJac.end(); + for( int cellIdx = 0; cellIdx < numCells; ++cellIdx ) + { + const double cellVolume = ebosSimulator_.model().dofTotalVolume(cellIdx); + auto& cellRes = ebosResid[ cellIdx ]; + + for( int flowPhaseIdx = 0; flowPhaseIdx < numPhases; ++flowPhaseIdx ) + { + const double refDens = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( flowPhaseIdx ), 0 ); + cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] /= refDens; + cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] *= cellVolume; + } + } + + for( auto row = ebosJac.begin(); row != endrow; ++row ) + { + const int rowIdx = row.index(); + const double cellVolume = ebosSimulator_.model().dofTotalVolume(rowIdx); + + + // translate the Jacobian of the residual from the format used by ebos to + // the one expected by flow + const auto endcol = row->end(); + for( auto col = row->begin(); col != endcol; ++col ) + { + for( int flowPhaseIdx = 0; flowPhaseIdx < numPhases; ++flowPhaseIdx ) + { + const double refDens = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( flowPhaseIdx ), 0 ); + for( int pvIdx=0; pvIdx& Je = eqs[eq].derivative(); + Sp Bb; const M& B = Je[n]; + B.toSparse(Bb); + //std::cout << "B eigen" << std::endl; + //std::cout << Bb << std::endl; + // Update right hand side. vals[eq] = eqs[eq].value().matrix() - B * Dibn; + //std::cout << "vals " << eq << std::endl; + //std::cout << vals[eq][0] << std::endl; } for (int var = 0; var < num_eq; ++var) { if (var == n) { @@ -149,6 +156,11 @@ namespace Opm V equation_value = equation.value(); ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs))); const M& C = eq_coll.derivative()[0]; + typedef Eigen::SparseMatrix Sp; + Sp Cc; + C.toSparse(Cc); + //std::cout << "C eigen" << std::endl; + //std::cout << Cc < Sp; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index d1b2a1f41..1ebb7eb08 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -47,9 +47,12 @@ #include #include #include +#include - +#include +#include +#include #include #include @@ -110,6 +113,11 @@ namespace Opm { using Vector = ADB::V; using V = ADB::V; typedef ADB::M M; + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; // copied from BlackoilModelBase @@ -138,6 +146,10 @@ namespace Opm { , param_(param) , terminal_output_(terminal_output) , pv_(pv) + //, resWell(wells_->number_of_wells) + //, duneD( Mat(2, 2, 3, 0.4, Mat::implicit)) + //, duneB( Mat(300, 2, 6, 0.4, Mat::implicit)) + //, duneC( Mat(2, 300, 6, 0.4, Mat::implicit)) { } @@ -163,21 +175,140 @@ namespace Opm { return iter_report; } //wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells); + if (param_.solve_welleq_initially_ && iterationIdx == 0) { // solve the well equations as a pre-processing step iter_report = solveWellEq(ebosSimulator, dt, well_state); } + std::vector cq_s; - computeWellFluxDense(ebosSimulator, cq_s, 4); + int nw = wells_->number_of_wells; + int nc = numCells(); + Mat duneD (nw, nw, 9, 0.4, Mat::implicit); + Mat duneB (nc, nw, 9, 0.4, Mat::implicit); + Mat duneC (nw, nc, 9, 0.4, Mat::implicit); + Mat duneA (nc, nc, 9, 0.4, Mat::implicit); + BVector rhs(nc); + + computeWellFluxDense(ebosSimulator, cq_s, 4, duneA, duneB, duneC, duneD, rhs); + //std::cout << cq_s[0] << std::endl; + //std::cout << cq_s[1] << std::endl; + //std::cout << cq_s[2] << std::endl; + + updatePerfPhaseRatesAndPressures(cq_s, well_state); - addWellFluxEq(cq_s, dt, 4, residual); - addWellContributionToMassBalanceEq(cq_s, residual); + addWellContributionToMassBalanceEq(cq_s,residual); + + addWellFluxEq(cq_s, dt, 4, residual, duneD); + //std::cout << residual.well_flux_eq << std::endl; + duneA.compress(); + duneB.compress(); + duneC.compress(); + duneD.compress(); + + //std::cout << "B" << std::endl; + //print(duneB); + //std::cout << "C" << std::endl; + //print(duneC); + //print(duneD); + + + V resWellEigen = residual.well_flux_eq.value(); + const int np = numPhases(); + BVector resWell(nw); + for (int i = 0; i < nw; ++i){ + for( int p = 0; p < np; ++p ) { + int idx = i + p * nw; + resWell[i][flowPhaseToEbosCompIdx(p)] = resWellEigen(idx); + } + } + + resWell_ = resWell; + rhs_ = rhs; + duneB_ = duneB; + duneC_ = duneC; + localInvert(duneD); + invDuneD_ = duneD; + duneA_ = duneA; + if (param_.compute_well_potentials_) { //wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); } return iter_report; } + void localInvert(Mat& istlA) const { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + //std::cout << (*col) << std::endl; + (*col).invert(); + } + } + } + + void print(Mat& istlA) const { + for (auto row = istlA.begin(), rowend = istlA.end(); row != rowend; ++row ) { + for (auto col = row->begin(), colend = row->end(); col != colend; ++col ) { + std::cout << (*col) << std::endl; + } + } + } + void addRhs(BVector& x, Mat& jac) const { + assert(x.size() == rhs.size()); + x += rhs_; + jac += duneA_; + } + + void apply( Mat& A, + BVector& res) const { + + Mat BmultinvD; + Mat duneA; + + Dune::matMultMat(BmultinvD, duneB_ , invDuneD_); + Dune::matMultMat(duneA, BmultinvD, duneC_); + A -= duneA; + BmultinvD.mmv(resWell_, res); + } + + void recoverVariable(const BVector& x, BVector& xw) const { + BVector resWell = resWell_; + duneC_.mmv(x, resWell); + invDuneD_.mv(resWell, xw); + } + + + + void addDuneMatrix(Eigen::SparseMatrix eigen, Mat& dune, int p1, int p2) { + //Mat dune( eigen.rows(), eigen.cols(), eigen.nonZeros(), 0.4,Mat::implicit) ; + const int* ia = eigen.outerIndexPtr(); + const int* ja = eigen.innerIndexPtr(); + for (int row = 0; row < eigen.rows(); ++row) { + dune.entry(ia[row],ja[row])[p1][p2] = eigen.coeff(ia[row] ,ja[row]); + } + } + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + return phaseToComp[ phaseIdx ]; + } + + int flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + return flowToEbos[ flowPv ]; + } + + int ebosCompToFlowPhaseIdx( const int compIdx ) const + { + const int compToPhase[ 3 ] = { Oil, Water, Gas }; + return compToPhase[ compIdx ]; + } @@ -193,7 +324,7 @@ namespace Opm { phase_condition_ = pc_arg; vfp_properties_ = vfp_properties_arg; gravity_ = gravity_arg; - perf_cell_depth_ = subset(depth_arg, wellOps().well_cells); + perf_cell_depth_ = subset(depth_arg, wellOps().well_cells); } const WellOps& wellOps() const @@ -429,14 +560,9 @@ namespace Opm { addWellFluxEq(std::vector cq_s, const double dt, const int numBlocks, - LinearisedBlackoilResidual& residual) + LinearisedBlackoilResidual& residual, + Mat& duneD) { - if( !localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } const int np = wells().number_of_phases; const int nw = wells().number_of_wells; @@ -458,13 +584,14 @@ namespace Opm { // res -= cq_s[perf*np + p]; //} res_vec[w] = res; - + for (int i = 0; i < np; ++i) { + duneD.entry(w,w)[flowPhaseToEbosCompIdx(p)][flowToEbosPvIdx(i)] += res.derivatives[i+3]; + } } + ADB tmp = convertToADBWell(res_vec, numBlocks); qs += superset(tmp,Span(nw,1,p*nw), nw*np); } - //std::cout << residual_.well_flux_eq << std::endl; - //wellModel().convertToADB(res_vec, well_cells, nc, well_id, nw, numBlocks); //ADB qs = state.qs; @@ -474,7 +601,6 @@ namespace Opm { } residual.well_flux_eq = qs; - //std::cout << "etter Ebos" << residual_.well_flux_eq << std::endl; } @@ -521,7 +647,12 @@ namespace Opm { void computeWellFluxDense(const Simulator& ebosSimulator, std::vector& cq_s, - const int numBlocks) const + const int numBlocks, + Mat& duneA, + Mat& duneB, + Mat& duneC, + Mat& duneD, + BVector& rhs) const { if( ! localWellsActive() ) return ; const int np = wells().number_of_phases; @@ -668,10 +799,22 @@ namespace Opm { EvalWell cqt_is = cqt_i/volumeRatio; //std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl; for (int phase = 0; phase < np; ++phase) { - cq_s_dense[phase][perf] = cmix_s[phase] * cqt_is; // * b_perfcells_dense[phase]; + cq_s_dense[phase][perf] = cmix_s[phase] * cqt_is; // * b_perfcells_dense[phase]; } } + + for (int p1 = 0; p1 < np; ++p1) { + EvalWell tmp = cq_s_dense[p1][perf]; + for (int p2 = 0; p2 < np; ++p2) { + duneC.entry(w, cell_idx)[flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= tmp.derivatives[p2]; + duneB.entry(cell_idx, w)[flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= tmp.derivatives[p2+3]; + duneA.entry(cell_idx, cell_idx)[flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= tmp.derivatives[p2]; + duneD.entry(w, w)[flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= tmp.derivatives[p2+3]; + } + rhs[cell_idx][flowPhaseToEbosCompIdx(p1)] -= tmp.value; + } } + } cq_s.resize(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { @@ -700,10 +843,31 @@ namespace Opm { bool converged; do { // bhp and Q for the wells - computeWellFluxDense(ebosSimulator, cq_s, 1); + int nw = wells().number_of_wells; + int nc = numCells(); + Mat duneDlo( nw, nw, 9, 0.4,Mat::implicit); + Mat duneBlo( nc, nw, 9, 0.4, Mat::implicit); + Mat duneClo( nw, nc, 9, 0.4, Mat::implicit); + Mat duneAlo( nc, nc, 9, 0.4, Mat::implicit); + BVector rhslo(nc); + computeWellFluxDense(ebosSimulator, cq_s, 1, duneAlo, duneBlo, duneClo, duneDlo, rhslo); updatePerfPhaseRatesAndPressures(cq_s, well_state); - addWellFluxEq(cq_s, dt, 1, residual); - converged = getWellConvergence(ebosSimulator, it, residual); + addWellFluxEq(cq_s, dt, 1, residual, duneDlo); + V resWellEigen = residual.well_flux_eq.value(); + //std::cout << "resWellEigen " << resWellEigen << std::endl; + BVector resWell(nw); + for (int i = 0; i < nw; ++i){ + for( int p = 0; p < np; ++p ) { + int idx = i + p * nw; + resWell[i][flowPhaseToEbosCompIdx(p)] = resWellEigen(idx); + } + } + duneDlo.compress(); + //print(duneDlo); + localInvert(duneDlo); + + resWell_ = resWell; + converged = getWellConvergence(ebosSimulator, it); if (converged) { break; @@ -712,20 +876,46 @@ namespace Opm { ++it; if( localWellsActive() ) { - std::vector eqs; - eqs.reserve(1); - eqs.push_back(residual.well_flux_eq); - //eqs.push_back(residual_.well_eq); - ADB total_residual = vertcatCollapseJacs(eqs); - const std::vector& Jn = total_residual.derivative(); - typedef Eigen::SparseMatrix Sp; - Sp Jn0; - Jn[0].toSparse(Jn0); - const Eigen::SparseLU< Sp > solver(Jn0); - ADB::V total_residual_v = total_residual.value(); - const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); +// std::vector eqs; +// eqs.reserve(1); +// eqs.push_back(residual.well_flux_eq); +// //eqs.push_back(residual_.well_eq); +// ADB total_residual = vertcatCollapseJacs(eqs); +// const std::vector& Jn = total_residual.derivative(); +// typedef Eigen::SparseMatrix Sp; +// Sp Jn0; +// Jn[0].toSparse(Jn0); +// std::cout << Jn0 << std::endl; +// const Eigen::SparseLU< Sp > solver(Jn0); +// ADB::V total_residual_v = total_residual.value(); +// std::cout << "tot res " < residual() { + + const int np = numPhases(); + const int nw = wells().number_of_wells; + std::vector res(np*nw); + for( int p=0; p bool getWellConvergence(Simulator& ebosSimulator, - const int iteration, - const LinearisedBlackoilResidual& residual) + const int iteration) { const int np = numPhases(); const int nc = numCells(); @@ -773,7 +976,7 @@ namespace Opm { B.col(idx) = b; } - detail::convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc, np, pv_, residual); + detail::convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc, np, pv_, residual()); std::vector well_flux_residual(np); bool converged_Well = true; @@ -1533,16 +1736,6 @@ namespace Opm { return flowToEbos[ phaseIdx ]; } - int flowToEbosPvIdx( const int flowPv ) const - { - const int flowToEbos[ 3 ] = { - BlackoilIndices::pressureSwitchIdx, - BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx - }; - return flowToEbos[ flowPv ]; - } - /// upate the dynamic lists related to economic limits template void @@ -1698,6 +1891,12 @@ namespace Opm { std::vector wellVariables_; std::vector F0_; const std::vector& pv_; + Mat duneA_; + Mat duneB_; + Mat duneC_; + Mat invDuneD_; + BVector resWell_; + BVector rhs_; // protected methods EvalWell getBhp(const int wellIdx) const {