From f7e2d88fd967a2c5e2b1b2a22c9a21741ace6b3f Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 9 May 2012 15:06:13 +0200 Subject: [PATCH] Added Newton compressible fluid solver. --- opm/core/pressure/IncompTpfa.cpp | 100 ++++++++++++++++++++++++++++++ opm/core/pressure/IncompTpfa.hpp | 22 +++++++ opm/core/pressure/tpfa/ifs_tpfa.c | 53 ++++++++++++++++ opm/core/pressure/tpfa/ifs_tpfa.h | 17 +++++ 4 files changed, 192 insertions(+) diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp index cd0a07b6b..36bb2fcb2 100644 --- a/opm/core/pressure/IncompTpfa.cpp +++ b/opm/core/pressure/IncompTpfa.cpp @@ -26,6 +26,7 @@ #include #include #include +#include namespace Opm { @@ -233,9 +234,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& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + const std::vector& porevol, + const std::vector& rock_comp, + const std::vector& prev_pressure, + const std::vector& initial_porevol, + const double dt, + std::vector& pressure_increment) + { + UnstructuredGrid* gg = const_cast(&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& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + std::vector& pressure, + std::vector& faceflux, + std::vector& well_bhp, + std::vector& well_rate) { + + UnstructuredGrid* gg = const_cast(&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 diff --git a/opm/core/pressure/IncompTpfa.hpp b/opm/core/pressure/IncompTpfa.hpp index f40d3e232..f14d9108f 100644 --- a/opm/core/pressure/IncompTpfa.hpp +++ b/opm/core/pressure/IncompTpfa.hpp @@ -121,6 +121,28 @@ namespace Opm std::vector& well_bhp, std::vector& well_rate); + void solveIncrement(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + const std::vector& porevol, + const std::vector& rock_comp, + const std::vector& prev_pressure, + const std::vector& initial_porevol, + const double dt, + std::vector& pressure_increment); + + void computeFaceFlux(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + std::vector& pressure, + std::vector& faceflux, + std::vector& well_bhp, + std::vector& well_rate); + /// Expose read-only reference to internal half-transmissibility. const ::std::vector& getHalfTrans() const { return htrans_; } diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index c3566c931..22a6e4826 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -11,6 +11,16 @@ #include +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 diff --git a/opm/core/pressure/tpfa/ifs_tpfa.h b/opm/core/pressure/tpfa/ifs_tpfa.h index 8f8aed757..c06c2f660 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.h +++ b/opm/core/pressure/tpfa/ifs_tpfa.h @@ -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 ,