20#ifndef OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED
21#define OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED
26#include <opm/core/linalg/LinearSolverIstl.hpp>
29#include <opm/common/ErrorMacros.hpp>
30#include <opm/core/utility/SparseTable.hpp>
40 template <
class GridInterface,
53 : pgrid_(0), prock_(0), pfluid_(0)
60 void init(
const Opm::parameter::ParameterGroup& param)
63 typename FluidInterface::CompVec mix(0.0);
64 const int nc = FluidInterface::numComponents;
65 double inflow_mixture_gas = param.getDefault(
"inflow_mixture_gas", 1.0);
66 double inflow_mixture_oil = param.getDefault(
"inflow_mixture_oil", 0.0);
69 mix[FluidInterface::Gas] = inflow_mixture_gas;
70 mix[FluidInterface::Oil] = inflow_mixture_oil;
73 double inflow_mixture_water = param.getDefault(
"inflow_mixture_water", 0.0);
74 mix[FluidInterface::Water] = inflow_mixture_water;
75 mix[FluidInterface::Gas] = inflow_mixture_gas;
76 mix[FluidInterface::Oil] = inflow_mixture_oil;
80 OPM_THROW(std::runtime_error,
"Unhandled number of components: " << nc);
82 inflow_mixture_ = mix;
83 linsolver_.reset(
new LinearSolverIstl(param));
84 flux_rel_tol_ = param.getDefault(
"flux_rel_tol", 1e-5);
85 press_rel_tol_ = param.getDefault(
"press_rel_tol", 1e-5);
86 max_num_iter_ = param.getDefault(
"max_num_iter", 15);
87 max_relative_voldiscr_ = param.getDefault(
"max_relative_voldiscr", 0.15);
88 relax_time_voldiscr_ = param.getDefault(
"relax_time_voldiscr", 0.0);
89 relax_weight_pressure_iteration_ = param.getDefault(
"relax_weight_pressure_iteration", 1.0);
90 experimental_jacobian_ = param.getDefault(
"experimental_jacobian",
false);
91 nonlinear_residual_tolerance_ = param.getDefault(
"nonlinear_residual_tolerance", 0.0);
92 output_residual_ = param.getDefault(
"output_residual",
false);
93 voldisc_factor_ = param.getDefault(
"voldisc_factor", 1.0);
103 return inflow_mixture_;
137 void setup(
const GridInterface& grid,
138 const RockInterface& rock,
139 const FluidInterface& fluid,
140 const WellsInterface& wells,
141 const typename GridInterface::Vector& grav,
142 const BCInterface& bc,
143 const std::vector<typename FluidInterface::PhaseVec>* face_pressure = 0)
152 const double* perm = &(rock.permeability(0)(0,0));
154 poro_.resize(grid.numCells(), 1.0);
155 for (
int i = 0; i < grid.numCells(); ++i) {
156 poro_[i] = rock.porosity(i);
159 psolver_.
init(grid, wells, perm, &poro_[0], grav);
162 int num_faces = grid.numFaces();
166 bcvalues_.resize(num_faces, 0.0);
167 for (
int face = 0; face < num_faces; ++face) {
168 int bid = pgrid_->boundaryId(face);
173 FlowBC face_bc = bc.flowCond(bid);
176 bcvalues_[face] = face_bc.
pressure();
178 bcvalues_[face] = (*face_pressure)[face][FluidInterface::Liquid];
182 bcvalues_[face] = face_bc.
outflux();
183 if (bcvalues_[face] != 0.0) {
184 OPM_THROW(std::runtime_error,
"Nonzero Neumann conditions not yet properly implemented "
185 "(signs must be fixed, also face pressures are not correctly computed for this case)");
188 OPM_THROW(std::runtime_error,
"Unhandled boundary condition type.");
198 int num_wells = pwells_->numWells();
199 for (
int well = 0; well < num_wells; ++well) {
200 int num_perf = pwells_->numPerforations(well);
201 for (
int perf = 0; perf < num_perf; ++perf) {
202 int cell = pwells_->wellCell(well, perf);
203 perf_wells_.push_back(well);
204 perf_cells_.push_back(cell);
207 int num_perf = perf_wells_.size();
208 perf_A_.resize(num_perf*numPhases*numComponents);
209 perf_mob_.resize(num_perf*numPhases);
210 perf_sat_.resize(num_perf);
218 return max_relative_voldiscr_;
233 const std::vector<typename FluidInterface::PhaseVec>& face_pressure,
234 const std::vector<double>& well_perf_pressure,
235 const std::vector<typename FluidInterface::CompVec>& cell_z,
238 computeFluidProps(cell_pressure, face_pressure, well_perf_pressure, cell_z, dt);
240 if (rel_voldiscr > max_relative_voldiscr_) {
241 std::cout <<
" Relative volume discrepancy too large: " << rel_voldiscr << std::endl;
244 std::cout <<
" Relative volume discrepancy ok: " << rel_voldiscr << std::endl;
291 std::vector<typename FluidInterface::PhaseVec>& face_pressure,
292 const std::vector<typename FluidInterface::CompVec>& cell_z,
293 std::vector<double>& face_flux,
294 std::vector<double>& well_bhp_pressures,
295 std::vector<double>& well_perf_pressures,
296 std::vector<double>& well_perf_fluxes,
297 const std::vector<double>& src,
304 int num_cells = cell_pressure.size();
305 int num_faces = face_pressure.size();
306 state.cell_pressure.resize(num_cells);
307 for (
int cell = 0; cell < num_cells; ++cell) {
308 state.cell_pressure[cell] = cell_pressure[cell][FluidInterface::Liquid];
310 state.face_pressure.resize(num_faces);
311 for (
int face = 0; face < num_faces; ++face) {
312 state.face_pressure[face] = face_pressure[face][FluidInterface::Liquid];
314 state.face_flux.clear();
315 state.face_flux.resize(num_faces);
316 state.well_bhp_pressure = well_bhp_pressures;
317 state.well_perf_pressure = well_perf_pressures;
318 state.well_perf_flux = well_perf_fluxes;
322 if (nonlinear_residual_tolerance_ == 0.0) {
324 retcode = solveImpl(cell_z, src, dt, state);
328 retcode = solveImplNew(cell_z, src, dt, state);
335 for (
int cell = 0; cell < num_cells; ++cell) {
336 cell_pressure[cell] = state.cell_pressure[cell];
338 for (
int face = 0; face < num_faces; ++face) {
339 face_pressure[face] = state.face_pressure[face];
341 face_flux = state.face_flux;
342 well_bhp_pressures = state.well_bhp_pressure;
343 well_perf_pressures = state.well_perf_pressure;
344 well_perf_fluxes = state.well_perf_flux;
357 &(pfluid_->surfaceDensities()[0]));
363 void doStepIMPES(std::vector<typename FluidInterface::CompVec>& cell_z,
372 const GridInterface* pgrid_;
373 const RockInterface* prock_;
374 const FluidInterface* pfluid_;
375 const WellsInterface* pwells_;
376 typename GridInterface::Vector gravity_;
379 std::vector<double> poro_;
381 std::unique_ptr<LinearSolverIstl> linsolver_;
382 std::vector<PressureAssembler::FlowBCTypes> bctypes_;
383 std::vector<double> bcvalues_;
385 typename FluidInterface::CompVec inflow_mixture_;
386 double flux_rel_tol_;
387 double press_rel_tol_;
389 double max_relative_voldiscr_;
390 double relax_time_voldiscr_;
391 double relax_weight_pressure_iteration_;
392 bool experimental_jacobian_;
393 bool output_residual_;
394 double nonlinear_residual_tolerance_;
395 double voldisc_factor_;
397 typedef typename FluidInterface::PhaseVec PhaseVec;
398 typedef typename FluidInterface::CompVec CompVec;
399 enum { numPhases = FluidInterface::numPhases,
400 numComponents = FluidInterface::numComponents };
402 std::vector<int> perf_wells_;
403 std::vector<int> perf_cells_;
404 std::vector<double> perf_A_;
405 std::vector<double> perf_mob_;
406 std::vector<PhaseVec> perf_sat_;
407 std::vector<double> perf_gpot_;
410 mutable std::vector<PhaseVec> helper_cell_pressure_;
411 mutable std::vector<PhaseVec> helper_face_pressure_;
415 std::vector<double> cell_pressure;
416 std::vector<double> face_pressure;
417 std::vector<double> face_flux;
418 std::vector<double> well_bhp_pressure;
419 std::vector<double> well_perf_pressure;
420 std::vector<double> well_perf_flux;
424 void computeFluidProps(
const std::vector<typename FluidInterface::PhaseVec>& phase_pressure,
425 const std::vector<typename FluidInterface::PhaseVec>& phase_pressure_face,
426 const std::vector<double>& well_perf_pressure,
427 const std::vector<typename FluidInterface::CompVec>& cell_z,
430 fp_.
computeNew(*pgrid_, *prock_, *pfluid_, gravity_, phase_pressure, phase_pressure_face, cell_z, inflow_mixture_, dt);
435 unsigned int perfcount = 0;
436 int num_wells = pwells_->numWells();
437 for (
int well = 0; well < num_wells; ++well) {
438 bool inj = pwells_->type(well) == WellsInterface::Injector;
439 int num_perf = pwells_->numPerforations(well);
440 for (
int perf = 0; perf < num_perf; ++perf) {
441 int cell = pwells_->wellCell(well, perf);
443 PhaseVec well_pressure = inj ? PhaseVec(well_perf_pressure[perfcount]) : phase_pressure[cell];
444 CompVec well_mixture = inj ? pwells_->injectionMixture(cell) : cell_z[cell];
445 typename FluidInterface::FluidState state = pfluid_->computeState(well_pressure, well_mixture);
446 std::copy(&state.phase_to_comp_[0][0], &state.phase_to_comp_[0][0] + numComponents*numPhases,
447 &perf_A_[perfcount*numPhases*numComponents]);
448 std::copy(state.mobility_.begin(), state.mobility_.end(),
449 &perf_mob_[perfcount*numPhases]);
450 perf_sat_[perfcount] = state.saturation_;
454 assert(perfcount == perf_wells_.size());
460 void computeFluidPropsScalarPress(
const std::vector<double>& cell_pressure_scalar,
461 const std::vector<double>& face_pressure_scalar,
462 const std::vector<double>& well_perf_pressure,
463 const std::vector<typename FluidInterface::CompVec>& cell_z,
467 int num_cells = pgrid_->numCells();
468 int num_faces = pgrid_->numFaces();
469 helper_cell_pressure_.resize(num_cells);
470 helper_face_pressure_.resize(num_faces);
471 for (
int cell = 0; cell < num_cells; ++cell) {
472 helper_cell_pressure_[cell] = cell_pressure_scalar[cell];
474 for (
int face = 0; face < num_faces; ++face) {
475 helper_face_pressure_[face] = face_pressure_scalar[face];
477 computeFluidProps(helper_cell_pressure_, helper_face_pressure_, well_perf_pressure, cell_z, dt);
483 void computeLinearResidual(
const PressureAssembler::LinearSystem& s, std::vector<double>& res)
486 for (
int row = 0; row < s.n; ++row) {
487 res[row] = -s.b[row];
488 for (
int i = s.ia[row]; i < s.ia[row + 1]; ++i) {
489 res[row] += s.sa[i]*s.x[s.ja[i]];
497 void computeResidualJacobian(
const std::vector<double>& initial_voldiscr,
498 const std::vector<double>& cell_pressure_scalar,
499 const std::vector<double>& cell_pressure_scalar_initial,
500 const std::vector<double>& well_bhp,
501 const std::vector<double>& src,
503 PressureAssembler::LinearSystem& linsys,
504 std::vector<double>& res)
506 std::vector<double>
zero(initial_voldiscr.size(), 0.0);
508 psolver_.assemble(&src[0],
519 &cell_pressure_scalar_initial[0],
522 &(pfluid_->surfaceDensities()[0]));
523 psolver_.linearSystem(linsys);
526 int num_cells = pgrid_->numCells();
527 std::copy(cell_pressure_scalar.begin(), cell_pressure_scalar.end(), linsys.x);
528 std::copy(well_bhp.begin(), well_bhp.end(), linsys.x + num_cells);
529 computeLinearResidual(linsys, res);
532 for (
int cell = 0; cell < num_cells; ++cell) {
534 double dres = tc*(cell_pressure_scalar[cell] - cell_pressure_scalar_initial[cell]);
536 dres *= pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
540 for (
int cell = 0; cell < num_cells; ++cell) {
541 for (
int i = linsys.ia[cell]; i < linsys.ia[cell + 1]; ++i) {
542 if (linsys.ja[i] == cell) {
544 linsys.sa[i] -= tc*pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
546 linsys.sa[i] += ex*pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
556 void computeWellPotentials(std::vector<double>& wellperf_gpot)
const
558 int num_perf = perf_cells_.size();
559 wellperf_gpot.resize(num_perf*numPhases);
560 for (
int perf = 0; perf < num_perf; ++perf) {
561 int well = perf_wells_[perf];
562 int cell = perf_cells_[perf];
563 typename GridInterface::Vector pos = pgrid_->cellCentroid(cell);
565 assert(gravity_[0] == 0.0 && gravity_[1] == 0.0);
566 double depth_delta = pos[2] - pwells_->referenceDepth(well);
567 double gh = gravity_[2]*depth_delta;
569 const double* At = &perf_A_[perf*numPhases*numComponents];
570 PhaseVec rho = pfluid_->phaseDensities(At);
571 for (
int phase = 0; phase < numPhases; ++phase) {
573 wellperf_gpot[numPhases*perf + phase] = rho[phase]*gh;
581 static std::pair<double, double>
582 computeFluxPressChanges(
const std::vector<double>& face_flux,
583 const std::vector<double>& well_perf_fluxes,
584 const std::vector<double>& cell_pressure_scalar,
585 const std::vector<double>& start_face_flux,
586 const std::vector<double>& start_perf_flux,
587 const std::vector<double>& start_cell_press)
589 int num_faces = face_flux.size();
590 int num_perf = well_perf_fluxes.size();
591 int num_cells = cell_pressure_scalar.size();
592 double max_flux_face = std::max(std::fabs(*std::min_element(face_flux.begin(), face_flux.end())),
593 std::fabs(*std::max_element(face_flux.begin(), face_flux.end())));
594 double max_flux_perf = num_perf == 0 ? 0.0
595 : std::max(std::fabs(*std::min_element(well_perf_fluxes.begin(), well_perf_fluxes.end())),
596 std::fabs(*std::max_element(well_perf_fluxes.begin(), well_perf_fluxes.end())));
597 double max_flux = std::max(max_flux_face, max_flux_perf);
598 double max_press = std::max(std::fabs(*std::min_element(cell_pressure_scalar.begin(),
599 cell_pressure_scalar.end())),
600 std::fabs(*std::max_element(cell_pressure_scalar.begin(),
601 cell_pressure_scalar.end())));
602 double flux_change_infnorm = 0.0;
603 double press_change_infnorm = 0.0;
604 for (
int face = 0; face < num_faces; ++face) {
605 flux_change_infnorm = std::max(flux_change_infnorm,
606 std::fabs(face_flux[face] - start_face_flux[face]));
608 for (
int perf = 0; perf < num_perf; ++perf) {
609 flux_change_infnorm = std::max(flux_change_infnorm,
610 std::fabs(well_perf_fluxes[perf] - start_perf_flux[perf]));
612 for (
int cell = 0; cell < num_cells; ++cell) {
613 press_change_infnorm = std::max(press_change_infnorm,
614 std::fabs(cell_pressure_scalar[cell] - start_cell_press[cell]));
616 double flux_rel_difference = flux_change_infnorm/max_flux;
617 double press_rel_difference = press_change_infnorm/max_press;
618 return std::make_pair(flux_rel_difference, press_rel_difference);
624 void computeWellPerfPressures(
const std::vector<double>& well_perf_fluxes,
625 const std::vector<double>& well_bhp,
626 const std::vector<double>& wellperf_gpot,
627 std::vector<double>& well_perf_pressures)
const
632 int num_perf = well_perf_fluxes.size();
633 std::vector<PhaseVec> well_sat(pwells_->numWells(), PhaseVec(0.0));
634 std::vector<double> well_flux(pwells_->numWells(), 0.0);
635 for (
int perf = 0; perf < num_perf; ++perf) {
636 int well = perf_wells_[perf];
637 double flux = well_perf_fluxes[perf];
638 well_flux[well] += flux;
639 PhaseVec tmp = perf_sat_[perf];
641 well_sat[well] += tmp;
643 for (
int well = 0; well < pwells_->numWells(); ++well) {
644 well_sat[well] *= 1.0/well_flux[well];
648 for (
int perf = 0; perf < num_perf; ++perf) {
649 well_perf_pressures[perf] = well_bhp[perf_wells_[perf]];
650 PhaseVec sat = well_sat[perf_wells_[perf]];
651 for (
int phase = 0; phase < numPhases; ++phase) {
652 well_perf_pressures[perf]
653 += sat[phase]*wellperf_gpot[numPhases*perf + phase];
661 bool getInitialVolumeDiscrepancy(
const double dt, std::vector<double>& voldiscr_initial)
665 if (rel_voldiscr > max_relative_voldiscr_) {
666 std::cout <<
"Initial relative volume discrepancy too large: " << rel_voldiscr << std::endl;
669 if (relax_time_voldiscr_ > 0.0) {
670 double relax = std::min(1.0,dt/relax_time_voldiscr_);
671 std::transform(voldiscr_initial.begin(), voldiscr_initial.end(), voldiscr_initial.begin(),
672 std::binder1st<std::multiplies<double> >(std::multiplies<double>() , relax));
674 if (voldisc_factor_ != 1.0) {
675 std::transform(voldiscr_initial.begin(), voldiscr_initial.end(), voldiscr_initial.begin(),
676 std::binder1st<std::multiplies<double> >(std::multiplies<double>(), voldisc_factor_));
687 void relaxState(
const double weight,
688 const SolverState& start_state,
691 int num_cells = pgrid_->numCells();
692 int num_faces = pgrid_->numFaces();
693 for (
int cell = 0; cell < num_cells; ++cell) {
694 state.cell_pressure[cell] = weight*state.cell_pressure[cell] + (1.0-weight)*start_state.cell_pressure[cell];
696 for (
int face = 0; face < num_faces; ++face) {
697 state.face_pressure[face] = weight*state.face_pressure[face] + (1.0-weight)*start_state.face_pressure[face];
698 state.face_flux[face] = weight*state.face_flux[face] + (1.0-weight)*start_state.face_flux[face];
705 ReturnCode solveImpl(
const std::vector<typename FluidInterface::CompVec>& cell_z,
706 const std::vector<double>& src,
710 int num_cells = pgrid_->numCells();
711 std::vector<double> cell_pressure_initial = state.cell_pressure;
712 std::vector<double> voldiscr_initial;
713 SolverState start_state;
716 for (
int iter = 0; iter < max_num_iter_; ++iter) {
719 computeFluidPropsScalarPress(state.cell_pressure, state.face_pressure, state.well_perf_pressure, cell_z, dt);
723 bool voldisc_ok = getInitialVolumeDiscrepancy(dt, voldiscr_initial);
731 computeWellPotentials(perf_gpot_);
734 if (experimental_jacobian_) {
736 PressureAssembler::LinearSystem s;
737 std::vector<double> residual;
738 computeResidualJacobian(voldiscr_initial, state.cell_pressure, cell_pressure_initial,
739 state.well_bhp_pressure, src, dt, s, residual);
741 if (output_residual_) {
743 static int psolve_iter = -1;
747 std::ostringstream oss;
748 oss <<
"residual-" << psolve_iter <<
'-' << iter <<
".dat";
749 std::ofstream outres(oss.str().c_str());
750 std::copy(residual.begin(), residual.end(), std::ostream_iterator<double>(outres,
"\n"));
754 LinearSolverIstl::LinearSolverReport result
755 = linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, &residual[0], s.x);
756 if (!result.converged) {
757 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
758 <<
"Residual reduction achieved is " << result.residual_reduction <<
'\n');
762 for (
int cell = 0; cell < num_cells; ++cell) {
763 s.x[cell] = state.cell_pressure[cell] - s.x[cell];
765 for (
int well = 0; well < pwells_->numWells(); ++well) {
766 s.x[num_cells + well] = state.well_bhp_pressure[well] - s.x[num_cells + well];
770 psolver_.assemble(&src[0],
775 &voldiscr_initial[0],
781 &cell_pressure_initial[0],
784 &(pfluid_->surfaceDensities()[0]));
785 PressureAssembler::LinearSystem s;
786 psolver_.linearSystem(s);
788 LinearSolverIstl::LinearSolverReport res = linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x);
789 if (!res.converged) {
790 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << res.iterations <<
" iterations.\n"
791 <<
"Residual reduction achieved is " << res.residual_reduction <<
'\n');
796 psolver_.computePressuresAndFluxes(state.cell_pressure, state.face_pressure, state.face_flux,
797 state.well_bhp_pressure, state.well_perf_flux);
800 if (relax_weight_pressure_iteration_ != 1.0) {
801 relaxState(relax_weight_pressure_iteration_, start_state, state);
805 computeWellPerfPressures(state.well_perf_flux, state.well_bhp_pressure,
806 perf_gpot_, state.well_perf_pressure);
809 std::pair<double, double> rel_changes
810 = computeFluxPressChanges(state.face_flux, state.well_perf_flux, state.cell_pressure,
811 start_state.face_flux, start_state.well_perf_flux, start_state.cell_pressure);
812 double flux_rel_difference = rel_changes.first;
813 double press_rel_difference = rel_changes.second;
817 std::cout <<
"Iteration Rel. flux change Rel. pressure change\n";
819 std::cout.precision(5);
820 std::cout << std::setw(6) << iter
821 << std::setw(24) << flux_rel_difference
822 << std::setw(24) << press_rel_difference << std::endl;
823 std::cout.precision(16);
825 if (flux_rel_difference < flux_rel_tol_ || press_rel_difference < press_rel_tol_) {
826 std::cout <<
"Pressure solver converged. Number of iterations: " << iter + 1 <<
'\n' << std::endl;
838 ReturnCode solveImplNew(
const std::vector<typename FluidInterface::CompVec>& cell_z,
839 const std::vector<double>& src,
843 int num_cells = pgrid_->numCells();
844 int num_faces = pgrid_->numFaces();
845 std::vector<double> cell_pressure_initial = state.cell_pressure;
846 std::vector<double> voldiscr_initial;
847 SolverState start_state;
850 for (
int iter = 0; iter < max_num_iter_; ++iter) {
853 computeFluidPropsScalarPress(state.cell_pressure, state.face_pressure, state.well_perf_pressure, cell_z, dt);
857 bool voldisc_ok = getInitialVolumeDiscrepancy(dt, voldiscr_initial);
865 computeWellPotentials(perf_gpot_);
869 PressureAssembler::LinearSystem s;
870 std::vector<double> residual;
871 computeResidualJacobian(voldiscr_initial, state.cell_pressure, cell_pressure_initial,
872 state.well_bhp_pressure, src, dt, s, residual);
873 if (output_residual_) {
875 static int psolve_iter = -1;
879 std::ostringstream oss;
880 oss <<
"residual-" << psolve_iter <<
'-' << iter <<
".dat";
881 std::ofstream outres(oss.str().c_str());
882 std::copy(residual.begin(), residual.end(), std::ostream_iterator<double>(outres,
"\n"));
886 double maxres = std::max(std::fabs(*std::max_element(residual.begin(), residual.end())),
887 std::fabs(*std::min_element(residual.begin(), residual.end())));
890 LinearSolverIstl::LinearSolverReport result
891 = linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, &residual[0], s.x);
892 if (!result.converged) {
893 OPM_THROW(std::runtime_error,
"Linear solver failed to converge in " << result.iterations <<
" iterations.\n"
894 <<
"Residual reduction achieved is " << result.residual_reduction <<
'\n');
899 for (
int cell = 0; cell < num_cells; ++cell) {
900 s.x[cell] = state.cell_pressure[cell] - s.x[cell];
902 for (
int well = 0; well < pwells_->numWells(); ++well) {
903 s.x[num_cells + well] = state.well_bhp_pressure[well] - s.x[num_cells + well];
907 psolver_.computePressuresAndFluxes(state.cell_pressure, state.face_pressure, state.face_flux,
908 state.well_bhp_pressure, state.well_perf_flux);
911 if (relax_weight_pressure_iteration_ != 1.0) {
912 double ww = relax_weight_pressure_iteration_;
913 for (
int cell = 0; cell < num_cells; ++cell) {
914 state.cell_pressure[cell] = ww*state.cell_pressure[cell] + (1.0-ww)*start_state.cell_pressure[cell];
917 for (
int face = 0; face < num_faces; ++face) {
918 state.face_pressure[face] = ww*state.face_pressure[face] + (1.0-ww)*start_state.face_pressure[face];
919 state.face_flux[face] = ww*state.face_flux[face] + (1.0-ww)*start_state.face_flux[face];
925 computeWellPerfPressures(state.well_perf_flux, state.well_bhp_pressure,
926 perf_gpot_, state.well_perf_pressure);
929 std::pair<double, double> rel_changes
930 = computeFluxPressChanges(state.face_flux, state.well_perf_flux, state.cell_pressure,
931 start_state.face_flux, start_state.well_perf_flux, start_state.cell_pressure);
932 double flux_rel_difference = rel_changes.first;
933 double press_rel_difference = rel_changes.second;
937 std::cout <<
"Iteration Rel. flux change Rel. pressure change Residual\n";
939 std::cout.precision(5);
940 std::cout << std::setw(6) << iter
941 << std::setw(24) << flux_rel_difference
942 << std::setw(24) << press_rel_difference
943 << std::setw(24) << maxres << std::endl;
944 std::cout.precision(16);
946 if (maxres < nonlinear_residual_tolerance_) {
947 std::cout <<
"Pressure solver converged. Number of iterations: " << iter + 1 <<
'\n' << std::endl;
bool isDirichlet() const
Type query.
Definition: BoundaryConditions.hpp:85
bool isNeumann() const
Type query.
Definition: BoundaryConditions.hpp:91
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:122
double pressure() const
Query a Dirichlet condition.
Definition: BoundaryConditions.hpp:145
double outflux() const
Query a Neumann condition.
Definition: BoundaryConditions.hpp:154
Definition: TpfaCompressible.hpp:46
void setup(const GridInterface &grid, const RockInterface &rock, const FluidInterface &fluid, const WellsInterface &wells, const typename GridInterface::Vector &grav, const BCInterface &bc, const std::vector< typename FluidInterface::PhaseVec > *face_pressure=0)
Setup routine, does grid/rock-dependent initialization.
Definition: TpfaCompressible.hpp:137
double volumeDiscrepancyLimit() const
Definition: TpfaCompressible.hpp:216
ReturnCode solve(std::vector< typename FluidInterface::PhaseVec > &cell_pressure, std::vector< typename FluidInterface::PhaseVec > &face_pressure, const std::vector< typename FluidInterface::CompVec > &cell_z, std::vector< double > &face_flux, std::vector< double > &well_bhp_pressures, std::vector< double > &well_perf_pressures, std::vector< double > &well_perf_fluxes, const std::vector< double > &src, const double dt)
Construct and solve system of linear equations for the phase pressure values on cells and faces,...
Definition: TpfaCompressible.hpp:290
TpfaCompressible()
Default constructor. Does nothing.
Definition: TpfaCompressible.hpp:52
void init(const Opm::parameter::ParameterGroup ¶m)
Initializes run-time parameters of the solver.
Definition: TpfaCompressible.hpp:60
TpfaCompressibleAssembler PressureAssembler
Definition: TpfaCompressible.hpp:48
ReturnCode
Definition: TpfaCompressible.hpp:249
@ SolveOk
Definition: TpfaCompressible.hpp:249
@ VolumeDiscrepancyTooLarge
Definition: TpfaCompressible.hpp:249
@ FailedToConverge
Definition: TpfaCompressible.hpp:249
const std::vector< double > & faceTransmissibilities()
Definition: TpfaCompressible.hpp:224
void doStepIMPES(std::vector< typename FluidInterface::CompVec > &cell_z, const double dt)
Do an IMPES step using the facilities of the pressure solver.
Definition: TpfaCompressible.hpp:363
FluidInterface::CompVec inflowMixture() const
Accessor for the inflow mixture.
Definition: TpfaCompressible.hpp:101
bool volumeDiscrepancyAcceptable(const std::vector< typename FluidInterface::PhaseVec > &cell_pressure, const std::vector< typename FluidInterface::PhaseVec > &face_pressure, const std::vector< double > &well_perf_pressure, const std::vector< typename FluidInterface::CompVec > &cell_z, const double dt)
Definition: TpfaCompressible.hpp:232
double stableStepIMPES()
Call this function after solve().
Definition: TpfaCompressible.hpp:352
Encapsulates the cfs_tpfa (= compressible flow solver two-point flux approximation) solver modules.
Definition: TpfaCompressibleAssembler.hpp:39
void init(const Grid &grid, const Wells &wells, const double *perm, const double *porosity, const typename Grid::Vector &gravity)
Initialize the solver's structures for a given grid, for well setup also call initWells().
Definition: TpfaCompressibleAssembler.hpp:73
@ FBC_FLUX
Definition: TpfaCompressibleAssembler.hpp:138
@ FBC_PRESSURE
Definition: TpfaCompressibleAssembler.hpp:137
@ FBC_UNSET
Definition: TpfaCompressibleAssembler.hpp:136
double explicitTimestepLimit(const double *faceA, const double *phasemobf, const double *phasemobf_deriv, const double *surf_dens)
Explicit IMPES time step limit.
Definition: TpfaCompressibleAssembler.hpp:310
void explicitTransport(const double dt, double *cell_surfvols)
Explicit IMPES transport.
Definition: TpfaCompressibleAssembler.hpp:329
const std::vector< double > & faceTransmissibilities() const
Definition: TpfaCompressibleAssembler.hpp:389
Definition: BlackoilFluid.hpp:32
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
Definition: BlackoilFluid.hpp:406
FaceFluidData face_data
Definition: BlackoilFluid.hpp:413
void computeNew(const Grid &grid, const Rock &rock, const BlackoilFluid &fluid, const typename Grid::Vector gravity, const std::vector< PhaseVec > &cell_pressure, const std::vector< PhaseVec > &face_pressure, const std::vector< CompVec > &cell_z, const CompVec &bdy_z, const double dt)
Definition: BlackoilFluid.hpp:417
std::vector< double > voldiscr
Definition: BlackoilFluid.hpp:409
AllFluidStates cell_data
Definition: BlackoilFluid.hpp:408
std::vector< double > relvoldiscr
Definition: BlackoilFluid.hpp:410
std::vector< Scalar > total_compressibility
Definition: FluidStateBlackoil.hpp:78
std::vector< Scalar > total_phase_volume_density
Definition: FluidStateBlackoil.hpp:75
std::vector< Scalar > experimental_term
Definition: FluidStateBlackoil.hpp:79
std::vector< PhaseToCompMatrix > state_matrix
Definition: FluidStateBlackoil.hpp:73
std::vector< PhaseVec > gravity_potential
Definition: BlackoilFluid.hpp:401
std::vector< PhaseVec > mobility
Definition: BlackoilFluid.hpp:397
std::vector< PhaseJacobian > mobility_deriv
Definition: BlackoilFluid.hpp:398
std::vector< PhaseToCompMatrix > state_matrix
Definition: BlackoilFluid.hpp:394