merge.
This commit is contained in:
Kjetil Olsen Lye
commit
53dd538e99
@ -555,58 +555,69 @@ main(int argc, char** argv)
|
||||
if (rock_comp->isActive()) {
|
||||
rc.resize(num_cells);
|
||||
std::vector<double> initial_pressure = state.pressure();
|
||||
std::vector<double> initial_porevolume(num_cells);
|
||||
computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
|
||||
std::vector<double> pressure_increment(num_cells);
|
||||
std::vector<double> prev_pressure;
|
||||
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
|
||||
prev_pressure = state.pressure();
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
|
||||
}
|
||||
state.pressure() = initial_pressure;
|
||||
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
|
||||
state.pressure(), state.faceflux(), well_bhp, well_perfrates);
|
||||
double max_change = 0.0;
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
|
||||
}
|
||||
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
|
||||
if (max_change < nl_pressure_tolerance) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
|
||||
} 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();
|
||||
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
||||
ptime += pt;
|
||||
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, *rock_comp, state.pressure(), porevol);
|
||||
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);
|
||||
|
||||
if (check_well_controls) {
|
||||
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
|
||||
fractional_flows,
|
||||
well_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_iteration;
|
||||
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
|
||||
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
|
||||
}
|
||||
if (!well_control_passed) {
|
||||
std::cout << "Well controls not passed, solving again." << std::endl;
|
||||
} else {
|
||||
std::cout << "Well conditions met." << std::endl;
|
||||
}
|
||||
}
|
||||
} while (!well_control_passed);
|
||||
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]));
|
||||
}
|
||||
|
||||
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
|
||||
if (max_change < nl_pressure_tolerance) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
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();
|
||||
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
||||
ptime += pt;
|
||||
|
||||
if (check_well_controls) {
|
||||
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
|
||||
fractional_flows,
|
||||
well_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_iteration;
|
||||
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
|
||||
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
|
||||
}
|
||||
if (!well_control_passed) {
|
||||
std::cout << "Well controls not passed, solving again." << std::endl;
|
||||
} else {
|
||||
std::cout << "Well conditions met." << std::endl;
|
||||
}
|
||||
}
|
||||
} while (!well_control_passed);
|
||||
|
||||
// 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);
|
||||
if (!use_reorder) {
|
||||
// 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);
|
||||
if (!use_reorder) {
|
||||
clear_transport_source(tsrc);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
if (reorder_src[cell] > 0.0) {
|
||||
|
||||
@ -26,6 +26,7 @@
|
||||
#include <opm/core/linalg/sparse_sys.h>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/newwells.h>
|
||||
#include <string.h>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
@ -236,9 +237,108 @@ namespace Opm
|
||||
soln.well_flux = &well_rate[0];
|
||||
soln.well_press = &well_bhp[0];
|
||||
}
|
||||
|
||||
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
|
||||
}
|
||||
|
||||
|
||||
void IncompTpfa::solveIncrement(const std::vector<double>& totmob,
|
||||
const std::vector<double>& omega,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& wdp,
|
||||
const FlowBoundaryConditions* bcs,
|
||||
const std::vector<double>& porevol,
|
||||
const std::vector<double>& rock_comp,
|
||||
const std::vector<double>& prev_pressure,
|
||||
const std::vector<double>& initial_porevol,
|
||||
const double dt,
|
||||
std::vector<double>& pressure_increment)
|
||||
|
||||
{
|
||||
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
||||
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
|
||||
|
||||
if (!omega.empty()) {
|
||||
if (gpress_.empty()) {
|
||||
THROW("Nozero omega argument given, but gravity was null in constructor.");
|
||||
}
|
||||
mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
|
||||
&omega[0],
|
||||
&gpress_[0], &gpress_omegaweighted_[0]);
|
||||
} else {
|
||||
if (!gpress_.empty()) {
|
||||
THROW("Empty omega argument given, but gravity was non-null in constructor.");
|
||||
}
|
||||
}
|
||||
|
||||
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
|
||||
if (! src.empty()) { F.src = &src[0]; }
|
||||
F.bc = bcs;
|
||||
F.totmob = &totmob[0];
|
||||
F.wdp = &wdp[0];
|
||||
|
||||
if (rock_comp.empty()) {
|
||||
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
|
||||
} else {
|
||||
ifs_tpfa_assemble_comprock_increment(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
|
||||
&porevol[0], &rock_comp[0], dt, &prev_pressure[0],
|
||||
&initial_porevol[0], h_);
|
||||
}
|
||||
|
||||
linsolver_.solve(h_->A, h_->b, &pressure_increment[0]);
|
||||
|
||||
}
|
||||
|
||||
void IncompTpfa::computeFaceFlux(const std::vector<double>& totmob,
|
||||
const std::vector<double>& omega,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& wdp,
|
||||
const FlowBoundaryConditions* bcs,
|
||||
std::vector<double>& pressure,
|
||||
std::vector<double>& faceflux,
|
||||
std::vector<double>& well_bhp,
|
||||
std::vector<double>& well_rate) {
|
||||
|
||||
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
|
||||
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
|
||||
|
||||
if (!omega.empty()) {
|
||||
if (gpress_.empty()) {
|
||||
THROW("Nozero omega argument given, but gravity was null in constructor.");
|
||||
}
|
||||
mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
|
||||
&omega[0],
|
||||
&gpress_[0], &gpress_omegaweighted_[0]);
|
||||
} else {
|
||||
if (!gpress_.empty()) {
|
||||
THROW("Empty omega argument given, but gravity was non-null in constructor.");
|
||||
}
|
||||
}
|
||||
|
||||
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
|
||||
if (! src.empty()) { F.src = &src[0]; }
|
||||
F.bc = bcs;
|
||||
F.totmob = &totmob[0];
|
||||
F.wdp = &wdp[0];
|
||||
|
||||
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
|
||||
|
||||
faceflux.resize(grid_.number_of_faces);
|
||||
|
||||
ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
|
||||
soln.cell_press = &pressure[0];
|
||||
soln.face_flux = &faceflux[0];
|
||||
|
||||
if(wells_ != NULL) {
|
||||
well_bhp.resize(wells_->number_of_wells);
|
||||
well_rate.resize(wells_->well_connpos[ wells_->number_of_wells ]);
|
||||
soln.well_flux = &well_rate[0];
|
||||
soln.well_press = &well_bhp[0];
|
||||
}
|
||||
|
||||
memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
|
||||
|
||||
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
@ -121,6 +121,28 @@ namespace Opm
|
||||
std::vector<double>& well_bhp,
|
||||
std::vector<double>& well_rate);
|
||||
|
||||
void solveIncrement(const std::vector<double>& totmob,
|
||||
const std::vector<double>& omega,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& wdp,
|
||||
const FlowBoundaryConditions* bcs,
|
||||
const std::vector<double>& porevol,
|
||||
const std::vector<double>& rock_comp,
|
||||
const std::vector<double>& prev_pressure,
|
||||
const std::vector<double>& initial_porevol,
|
||||
const double dt,
|
||||
std::vector<double>& pressure_increment);
|
||||
|
||||
void computeFaceFlux(const std::vector<double>& totmob,
|
||||
const std::vector<double>& omega,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& wdp,
|
||||
const FlowBoundaryConditions* bcs,
|
||||
std::vector<double>& pressure,
|
||||
std::vector<double>& faceflux,
|
||||
std::vector<double>& well_bhp,
|
||||
std::vector<double>& well_rate);
|
||||
|
||||
/// Expose read-only reference to internal half-transmissibility.
|
||||
const ::std::vector<double>& getHalfTrans() const { return htrans_; }
|
||||
|
||||
|
||||
@ -11,6 +11,16 @@
|
||||
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
|
||||
|
||||
|
||||
void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v) {
|
||||
size_t i,j;
|
||||
for (j = 0; j < A->m; ++j) {
|
||||
v[j] = 0;
|
||||
for (i = (size_t) (A->ia[j]); i < (size_t) (A->ia[j+1]); ++i) {
|
||||
v[j] += A->sa[i]*u[A->ja[i]];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
struct ifs_tpfa_impl {
|
||||
double *fgrav; /* Accumulated grav contrib/face */
|
||||
|
||||
@ -725,6 +735,49 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G ,
|
||||
return ok;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int
|
||||
ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
|
||||
const struct ifs_tpfa_forces *F ,
|
||||
const double *trans ,
|
||||
const double *gpress ,
|
||||
const double *porevol ,
|
||||
const double *rock_comp,
|
||||
const double dt ,
|
||||
const double *prev_pressure,
|
||||
const double *initial_porevolume,
|
||||
struct ifs_tpfa_data *h )
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
int c;
|
||||
size_t j;
|
||||
int system_singular, ok;
|
||||
double *v;
|
||||
|
||||
assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
|
||||
|
||||
v = malloc(h->A->m * sizeof *v);
|
||||
|
||||
mult_csr_matrix(h->A, prev_pressure, v);
|
||||
|
||||
/* We want to solve a Newton step for the residual
|
||||
* (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible
|
||||
*
|
||||
*/
|
||||
|
||||
if (ok) {
|
||||
for (c = 0; c < G->number_of_cells; ++c) {
|
||||
j = csrmatrix_elm_index(c, c, h->A);
|
||||
h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
|
||||
h->b[c] += -(porevol[c] - initial_porevolume[c])/dt - v[c];
|
||||
}
|
||||
}
|
||||
|
||||
free(v);
|
||||
|
||||
return ok;
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void
|
||||
|
||||
@ -61,6 +61,10 @@ struct ifs_tpfa_data *
|
||||
ifs_tpfa_construct(struct UnstructuredGrid *G,
|
||||
struct Wells *W);
|
||||
|
||||
|
||||
void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v);
|
||||
|
||||
|
||||
int
|
||||
ifs_tpfa_assemble(struct UnstructuredGrid *G ,
|
||||
const struct ifs_tpfa_forces *F ,
|
||||
@ -78,6 +82,19 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G ,
|
||||
const double dt ,
|
||||
const double *pressure ,
|
||||
struct ifs_tpfa_data *h );
|
||||
int
|
||||
ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
|
||||
const struct ifs_tpfa_forces *F ,
|
||||
const double *trans ,
|
||||
const double *gpress ,
|
||||
const double *porevol ,
|
||||
const double *rock_comp,
|
||||
const double dt ,
|
||||
const double *prev_pressure ,
|
||||
const double *initial_porevolume,
|
||||
struct ifs_tpfa_data *h );
|
||||
|
||||
|
||||
void
|
||||
ifs_tpfa_press_flux(struct UnstructuredGrid *G ,
|
||||
const struct ifs_tpfa_forces *F ,
|
||||
|
||||
Loading…
Reference in New Issue
Block a user