diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index e1de61460..28ad657e2 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -28,7 +28,7 @@ #include #include #include -#include +#include #include #include #include @@ -53,9 +53,10 @@ #include #include -#include +#include #include +#include #include #include #include @@ -570,20 +571,19 @@ main(int argc, char** argv) // Check that rock compressibility is not used with solvers that do not handle it. int nl_pressure_maxiter = 0; - double nl_pressure_tolerance = 0.0; + double nl_pressure_residual_tolerance = 0.0; + double nl_pressure_change_tolerance = 0.0; if (rock_comp->isActive()) { if (!use_reorder) { THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet."); } + nl_pressure_residual_tolerance = param.getDefault("nl_pressure_residual_tolerance", 0.0); + nl_pressure_change_tolerance = param.getDefault("nl_pressure_change_tolerance", 1.0); // In Pascal. nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10); - nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal } // Source-related variables init. int num_cells = grid->c_grid()->number_of_cells; - std::vector totmob; - std::vector omega; // Will remain empty if no gravity. - std::vector rc; // Will remain empty if no rock compressibility. // Extra rock init. std::vector porevol; @@ -594,12 +594,8 @@ main(int argc, char** argv) } double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - // We need a separate reorder_sat, because the reorder - // code expects a scalar sw, not both sw and so. - std::vector reorder_sat(num_cells); - std::vector src(num_cells, 0.0); - // Initialising src + std::vector src(num_cells, 0.0); if (wells->c_wells()) { // Do nothing, wells will be the driving force, not source terms. // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); @@ -625,7 +621,10 @@ main(int argc, char** argv) Opm::LinearSolverFactory linsolver(param); // Pressure solver. const double *grav = use_gravity ? &gravity[0] : 0; - Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver, wells->c_wells()); + Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver, + nl_pressure_residual_tolerance, nl_pressure_change_tolerance, + nl_pressure_maxiter, + grav, wells->c_wells(), src, bcs.c_bcs()); // Reordering solver. const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const int nl_maxiter = param.getDefault("nl_maxiter", 30); @@ -639,7 +638,7 @@ main(int argc, char** argv) THROW("Unknown method: " << method_string); } - Opm::TransportModelPolymer reorder_model(*grid->c_grid(), props->porosity(), &porevol[0], *props, polyprop, + Opm::TransportModelPolymer reorder_model(*grid->c_grid(), *props, polyprop, method, nl_tolerance, nl_maxiter); if (use_gauss_seidel_gravity) { @@ -717,17 +716,15 @@ main(int argc, char** argv) Opm::Watercut watercut; watercut.push(0.0, 0.0, 0.0); Opm::WellReport wellreport; - std::vector well_bhp; - std::vector well_perfrates; + Opm::WellState well_state; + well_state.init(wells->c_wells(), state); std::vector fractional_flows; std::vector well_resflows_phase; int num_wells = 0; if (wells->c_wells()) { num_wells = wells->c_wells()->number_of_wells; - well_bhp.resize(num_wells, 0.0); - well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0); well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); - wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates); + wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); } for (; !simtimer.done(); ++simtimer) { // Report timestep and (optionally) write state to disk. @@ -737,17 +734,6 @@ main(int argc, char** argv) } // Solve pressure. - if (use_gravity) { - computeTotalMobilityOmega(*props, polyprop, allcells, state.saturation(), state.concentration(), state.maxconcentration(), - totmob, omega); - } else { - computeTotalMobility(*props, polyprop, allcells, state.saturation(), state.concentration(), state.maxconcentration(), - totmob); - } - std::vector wdp; - if (wells->c_wells()) { - Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp); - } if (check_well_controls) { computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); } @@ -758,48 +744,27 @@ main(int argc, char** argv) int well_control_iteration = 0; do { pressure_timer.start(); - if (rock_comp->isActive()) { - rc.resize(num_cells); - std::vector initial_pressure = state.pressure(); - std::vector initial_porevolume(num_cells); - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, initial_pressure, initial_porevolume); - std::vector pressure_increment(num_cells + num_wells); - std::vector prev_pressure(num_cells + num_wells); - for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { - - for (int cell = 0; cell < num_cells; ++cell) { - rc[cell] = rock_comp->rockComp(state.pressure()[cell]); - } - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); - std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); - std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells); - // prev_pressure = state.pressure(); - - // compute pressure increment - psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, - prev_pressure, initial_porevolume, simtimer.currentStepLength(), - pressure_increment); - - double max_change = 0.0; - for (int cell = 0; cell < num_cells; ++cell) { - state.pressure()[cell] += pressure_increment[cell]; - max_change = std::max(max_change, std::fabs(pressure_increment[cell])); - } - for (int well = 0; well < num_wells; ++well) { - well_bhp[well] += pressure_increment[num_cells + well]; - max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); - } - - std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; - if (max_change < nl_pressure_tolerance) { - break; - } + std::vector initial_pressure = state.pressure(); + psolver.solve(simtimer.currentStepLength(), state, well_state); + if (!rock_comp->isActive()) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.; + double av_press = 0.; + double tot_vol = 0.; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid->c_grid()->cell_volumes[cell]; + av_press += state.pressure()[cell]*grid->c_grid()->cell_volumes[cell]; + tot_vol += grid->c_grid()->cell_volumes[cell]; + } + // Renormalization constant + const double ren_const = (av_prev_press - av_press)/tot_vol; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += ren_const; + } + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; } - psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); - } else { - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); } pressure_timer.stop(); double pt = pressure_timer.secsSinceStart(); @@ -810,11 +775,11 @@ main(int argc, char** argv) if (check_well_controls) { Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), fractional_flows, - well_perfrates, + well_state.perfRates(), well_resflows_phase); std::cout << "Checking well conditions." << std::endl; // For testing we set surface := reservoir - well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); + well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); ++well_control_iteration; if (!well_control_passed && well_control_iteration > max_well_control_iterations) { THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); @@ -829,7 +794,7 @@ main(int argc, char** argv) // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, - wells->c_wells(), well_perfrates, reorder_src); + wells->c_wells(), well_state.perfRates(), reorder_src); // Find inflow rate. @@ -851,19 +816,16 @@ main(int argc, char** argv) } for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { if (use_reorder) { - Opm::toWaterSat(state.saturation(), reorder_sat); reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0], stepsize, inflow_c, - &reorder_sat[0], &state.concentration()[0], &state.maxconcentration()[0]); - Opm::toBothSat(reorder_sat, state.saturation()); + state.saturation(), state.concentration(), state.maxconcentration()); Opm::computeInjectedProduced(*props, polyprop, state.saturation(), state.concentration(), state.maxconcentration(), reorder_src, simtimer.currentStepLength(), inflow_c, injected, produced, polyinj, polyprod); if (use_segregation_split) { if (use_column_solver) { if (use_gauss_seidel_gravity) { - reorder_model.solveGravity(columns, &porevol[0], stepsize, reorder_sat, + reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation(), state.concentration(), state.maxconcentration()); - Opm::toBothSat(reorder_sat, state.saturation()); } else { colsolver.solve(columns, stepsize, state.saturation(), state.concentration(), state.maxconcentration()); @@ -936,7 +898,7 @@ main(int argc, char** argv) if (wells->c_wells()) { wellreport.push(*props, *wells->c_wells(), state.saturation(), simtimer.currentTime() + simtimer.currentStepLength(), - well_bhp, well_perfrates); + well_state.bhp(), well_state.perfRates()); } } total_timer.stop();