diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 9be9915cd..b23f370a1 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -101,7 +101,7 @@ main(int argc, char* argv[]) Opm::LinearSolverFactory linsolver(param); - Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, *wells, linsolver); + Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver); Opm::BlackoilState state; initStateBasic(*g, props0, param, 0.0, state); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 914d88586..969a4e96a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -187,11 +188,13 @@ namespace Opm { FullyImplicitBlackoilSolver(const UnstructuredGrid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, const Wells& wells, const LinearSolverInterface& linsolver) : grid_ (grid) , fluid_ (fluid) , geo_ (geo) + , rock_comp_props_(rock_comp_props) , wells_ (wells) , linsolver_ (linsolver) , active_(activePhases(fluid.phaseUsage())) @@ -523,12 +526,14 @@ namespace Opm { const std::vector& sat = state.saturation; const ADB& rs = state.rs; + const ADB pv_mult = poroMult(press); + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; for (int phase = 0; phase < maxnp; ++phase) { if (active_[ phase ]) { const int pos = pu.phase_pos[ phase ]; rq_[pos].b = fluidReciprocFVF(phase, press, rs, cells_); - rq_[pos].accum[aix] = rq_[pos].b * sat[pos]; + rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos]; // DUMP(rq_[pos].b); // DUMP(rq_[pos].accum[aix]); } @@ -539,7 +544,7 @@ namespace Opm { const int po = pu.phase_pos[ Oil ]; const int pg = pu.phase_pos[ Gas ]; - rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; + rq_[pg].accum[aix] += pv_mult * state.rs * rq_[po].accum[aix]; // DUMP(rq_[pg].accum[aix]); } } @@ -707,7 +712,7 @@ namespace Opm { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; const ADB rs_perf = subset(state.rs, well_cells); - well_rates_all += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np); + well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np); // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs); residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs; } @@ -992,8 +997,9 @@ namespace Opm { { const int phase = canph_[ actph ]; const ADB mu = fluidViscosity(phase, state.pressure, state.rs, cells_); + const ADB tr_mult = transMult(state.pressure); - rq_[ actph ].mob = kr[ phase ] / mu; + rq_[ actph ].mob = tr_mult * kr[ phase ] / mu; const ADB rho = fluidDensity(phase, state.pressure, state.rs, cells_); const ADB gflux = grav_ * rho; @@ -1123,4 +1129,58 @@ namespace Opm { } + + + + ADB + FullyImplicitBlackoilSolver::poroMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V pm(n); + V dpm(n); + for (int i = 0; i < n; ++i) { + pm[i] = rock_comp_props_->poroMult(p.value()[i]); + dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]); + } + ADB::M dpm_diag = spdiag(dpm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dpm_diag * p.derivative()[block]; + } + return ADB::function(pm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + + + + + ADB + FullyImplicitBlackoilSolver::transMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V tm(n); + V dtm(n); + for (int i = 0; i < n; ++i) { + tm[i] = rock_comp_props_->transMult(p.value()[i]); + dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]); + } + ADB::M dtm_diag = spdiag(dtm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dtm_diag * p.derivative()[block]; + } + return ADB::function(tm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + } // namespace Opm diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 1c6e4b874..a3859f472 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -30,6 +30,7 @@ struct Wells; namespace Opm { class DerivedGeology; + class RockCompressibility; class LinearSolverInterface; class BlackoilState; class WellState; @@ -42,6 +43,7 @@ namespace Opm { FullyImplicitBlackoilSolver(const UnstructuredGrid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, const Wells& wells, const LinearSolverInterface& linsolver); @@ -98,6 +100,7 @@ namespace Opm { const UnstructuredGrid& grid_; const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; + const RockCompressibility* rock_comp_props_; const Wells& wells_; const LinearSolverInterface& linsolver_; // For each canonical phase -> true if active @@ -187,6 +190,12 @@ namespace Opm { ADB fluidRsMax(const ADB& p, const std::vector& cells) const; + + ADB + poroMult(const ADB& p) const; + + ADB + transMult(const ADB& p) const; }; } // namespace Opm diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index e689e6192..cb6d8be5d 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -249,7 +249,7 @@ namespace Opm bcs_(bcs), gravity_(gravity), geo_(grid_, props_, gravity_), - solver_(grid_, props_, geo_, *wells_manager.c_wells(), linsolver) + solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) /* param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10),