20#ifndef OPM_TPFACOMPRESSIBLEASSEMBLER_HEADER_INCLUDED
21#define OPM_TPFACOMPRESSIBLEASSEMBLER_HEADER_INCLUDED
24#include <opm/core/pressure/tpfa/cfs_tpfa.h>
25#include <opm/core/pressure/tpfa/trans_tpfa.h>
26#include <opm/core/linalg/sparse_sys.h>
27#include <opm/core/pressure/flow_bc.h>
28#include <opm/core/pressure/legacy_well.h>
29#include <opm/core/pressure/tpfa/compr_quant.h>
30#include <dune/grid/common/GridAdapter.hpp>
44 : state_(Uninitialized), data_(0), bc_(0)
46 wells_.number_of_wells = 0;
56 flow_conditions_destroy(bc_);
57 cfs_tpfa_destroy(data_);
72 template <
class Gr
id,
class Wells>
73 void init(
const Grid& grid,
const Wells& wells,
const double* perm,
const double* porosity,
74 const typename Grid::Vector& gravity)
77 init(grid, perm, porosity, gravity);
87 void init(
const Grid& grid,
const double* perm,
const double* porosity,
const typename Grid::Vector& gravity)
95 if (wells_.number_of_wells != 0) {
98 data_ = cfs_tpfa_construct(grid_.c_grid(), w, num_phases);
100 throw std::runtime_error(
"Failed to initialize cfs_tpfa solver.");
104 int num_cells = grid.numCells();
105 int ngconn = grid_.c_grid()->cell_facepos[num_cells];
106 ncf_.resize(num_cells);
107 for (
int cell = 0; cell < num_cells; ++cell) {
108 int num_local_faces = grid.numCellFaces(cell);
109 ncf_[cell] = num_local_faces;
111 htrans_.resize(ngconn);
112 tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
115 trans_.resize(grid_.numFaces());
116 tpfa_trans_compute(grid_.c_grid(), &htrans_[0], &trans_[0]);
119 porevol_.resize(num_cells);
120 for (
int i = 0; i < num_cells; ++i) {
121 porevol_[i] = porosity[i]*grid.cellVolume(i);
125 if (Grid::dimension != 3) {
126 throw std::logic_error(
"Only 3 dimensions supported currently.");
128 std::copy(gravity.begin(), gravity.end(), gravity_);
130 state_ = Initialized;
153 const double* bcvalues,
155 const double* totcompr,
156 const double* voldiscr,
159 const double* wellperfA,
160 const double* phasemobf,
161 const double* phasemobwellperf,
162 const double* cell_pressure,
163 const double* gravcapf,
164 const double* wellperf_gpot,
165 const double* surf_dens)
167 if (state_ == Uninitialized) {
168 throw std::runtime_error(
"Error in TpfaCompressibleAssembler::assemble(): You must call init() prior to calling assemble().");
170 UnstructuredGrid* g = grid_.c_grid();
173 gather_boundary_conditions(bctypes, bcvalues);
176 double* src =
const_cast<double*
>(sources);
179 well_t* wells = NULL;
180 well_control_t* wctrl = NULL;
181 struct completion_data* wcompl = NULL;
182 if (wells_.number_of_wells != 0) {
187 std::copy(wellperf_gpot, wellperf_gpot + well_gpot_storage_.size(), well_gpot_storage_.begin());
188 std::copy(wellperfA, wellperfA + well_A_storage_.size(), well_A_storage_.begin());
189 std::copy(phasemobwellperf, phasemobwellperf + well_phasemob_storage_.size(), well_phasemob_storage_.begin());
193 compr_quantities cq = { 3 ,
194 const_cast<double *
>(totcompr ) ,
195 const_cast<double *
>(voldiscr ) ,
196 const_cast<double *
>(cellA ) ,
197 const_cast<double *
>(faceA ) ,
198 const_cast<double *
>(phasemobf) };
201 cfs_tpfa_assemble(g, dt, wells, bc_, src,
202 &cq, &trans_[0], gravcapf,
204 cell_pressure, &porevol_[0],
206 phasemobf_.assign(phasemobf, phasemobf + grid_.numFaces()*3);
207 gravcapf_.assign(gravcapf, gravcapf + grid_.numFaces()*3);
235 if (state_ != Assembled) {
236 throw std::runtime_error(
"Error in TpfaCompressibleAssembler::linearSystem(): "
237 "You must call assemble() prior to calling linearSystem().");
240 s.
nnz = data_->A->nnz;
259 std::vector<double>& face_pressures,
260 std::vector<double>& face_fluxes,
261 std::vector<double>& well_pressures,
262 std::vector<double>& well_fluxes)
264 if (state_ != Assembled) {
265 throw std::runtime_error(
"Error in TpfaCompressibleAssembler::computePressuresAndFluxes(): "
266 "You must call assemble() (and solve the linear system) "
267 "prior to calling computePressuresAndFluxes().");
269 int num_cells = grid_.c_grid()->number_of_cells;
270 int num_faces = grid_.c_grid()->number_of_faces;
271 cell_pressures.clear();
272 cell_pressures.resize(num_cells, 0.0);
273 face_pressures.clear();
274 face_pressures.resize(num_faces, 0.0);
276 face_fluxes.resize(num_faces, 0.0);
281 well_t* wells = NULL;
282 struct completion_data* wcompl = NULL;
285 if (wells_.number_of_wells != 0) {
288 well_pressures.resize(wells_.number_of_wells);
289 well_fluxes.resize(well_cells_storage_.size());
290 wpress = &well_pressures[0];
291 wflux = &well_fluxes[0];
294 cfs_tpfa_press_flux(grid_.c_grid(),
296 np, &trans_[0], &phasemobf_[0], &gravcapf_[0],
298 data_, &cell_pressures[0], &face_fluxes[0],
300 cfs_tpfa_fpress(grid_.c_grid(), bc_, np, &htrans_[0],
301 &phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0],
302 &face_fluxes[0], &face_pressures[0]);
311 const double* phasemobf,
312 const double* phasemobf_deriv,
313 const double* surf_dens)
315 compr_quantities cq = { 3,
319 const_cast<double *
>(faceA) ,
320 const_cast<double *
>(phasemobf) };
321 return cfs_tpfa_impes_maxtime(grid_.c_grid(), &cq, &trans_[0], &porevol_[0], data_,
322 phasemobf_deriv, surf_dens, gravity_);
330 double* cell_surfvols)
334 well_t* wells = NULL;
335 if (wells_.number_of_wells != 0) {
338 cfs_tpfa_expl_mass_transport(grid_.c_grid(), wells, &wcompl_,
339 np, dt, &porevol_[0],
340 data_, cell_surfvols);
359 std::vector<double>& cell_fluxes)
361 if (state_ != Assembled) {
362 throw std::runtime_error(
"Error in TpfaCompressibleAssembler::faceFluxToCellFlux(): "
363 "You must call assemble() (and solve the linear system) "
364 "prior to calling faceFluxToCellFlux().");
366 const UnstructuredGrid& g = *(grid_.c_grid());
367 int num_cells = g.number_of_cells;
368 cell_fluxes.resize(g.cell_facepos[num_cells]);
369 for (
int cell = 0; cell < num_cells; ++cell) {
370 for (
int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
371 int face = g.cell_faces[hface];
372 bool pos = (g.face_cells[2*face] == cell);
373 cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
399 enum State { Uninitialized, Initialized, Assembled };
403 cfs_tpfa_data* data_;
407 std::vector<int> ncf_;
409 std::vector<double> htrans_;
410 std::vector<double> trans_;
412 std::vector<double> porevol_;
414 std::vector<double> phasemobf_;
416 std::vector<double> gravcapf_;
421 FlowBoundaryConditions *bc_;
425 std::vector<int> well_connpos_storage_;
426 std::vector<int> well_cells_storage_;
427 well_control_t wctrl_;
428 std::vector<well_type> wctrl_type_storage_;
429 std::vector<well_control> wctrl_ctrl_storage_;
430 std::vector<double> wctrl_target_storage_;
431 struct completion_data wcompl_;
432 std::vector<double> well_prodind_storage_;
433 std::vector<double> well_gpot_storage_;
434 std::vector<double> well_A_storage_;
435 std::vector<double> well_phasemob_storage_;
444 template <
class Wells>
445 void initWells(
const Wells& w)
447 int num_wells = w.numWells();
448 if (num_wells == 0) {
449 wells_.number_of_wells = 0;
452 wctrl_type_storage_.resize(num_wells);
453 wctrl_ctrl_storage_.resize(num_wells);
454 wctrl_target_storage_.resize(num_wells);
455 for (
int i = 0; i < num_wells; ++i) {
456 wctrl_type_storage_[i] = (w.type(i) == Wells::Injector) ? INJECTOR : PRODUCER;
457 wctrl_ctrl_storage_[i] = (w.control(i) == Wells::Rate) ? RATE : BHP;
458 wctrl_target_storage_[i] = w.target(i);
459 int num_perf = w.numPerforations(i);
460 well_connpos_storage_.push_back(well_cells_storage_.size());
461 for (
int j = 0; j < num_perf; ++j) {
462 well_cells_storage_.push_back(w.wellCell(i, j));
463 well_prodind_storage_.push_back(w.wellIndex(i, j));
466 well_connpos_storage_.push_back(well_cells_storage_.size());
467 int tot_num_perf = well_prodind_storage_.size();
468 well_gpot_storage_.resize(tot_num_perf*3);
469 well_A_storage_.resize(3*3*tot_num_perf);
470 well_phasemob_storage_.resize(3*tot_num_perf);
472 wells_.number_of_wells = num_wells;
473 wells_.well_connpos = &well_connpos_storage_[0];
474 wells_.well_cells = &well_cells_storage_[0];
476 wctrl_.type = &wctrl_type_storage_[0];
477 wctrl_.ctrl = &wctrl_ctrl_storage_[0];
478 wctrl_.target = &wctrl_target_storage_[0];
480 wcompl_.WI = &well_prodind_storage_[0];
481 wcompl_.gpot = &well_gpot_storage_[0];
482 wcompl_.A = &well_A_storage_[0];
483 wcompl_.phasemob = &well_phasemob_storage_[0];
488 gather_boundary_conditions(
const FlowBCTypes* bctypes ,
489 const double* bcvalues)
492 bc_ = flow_conditions_construct(0);
495 flow_conditions_clear(bc_);
500 for (std::size_t i = 0, nf = grid_.numFaces(); ok && (i < nf); ++i) {
502 ok = flow_conditions_append(BC_PRESSURE,
507 else if (bctypes[ i ] ==
FBC_FLUX) {
508 ok = flow_conditions_append(BC_FLUX_TOTVOL,
516 flow_conditions_destroy(bc_);
Encapsulates the cfs_tpfa (= compressible flow solver two-point flux approximation) solver modules.
Definition: TpfaCompressibleAssembler.hpp:39
void init(const Grid &grid, 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:87
const std::vector< int > & numCellFaces() const
Access the number of connections (faces) per cell. Deprecated, will be removed.
Definition: TpfaCompressibleAssembler.hpp:383
void linearSystem(LinearSystem &s)
Access the linear system assembled. You must call assemble() prior to calling linearSystem().
Definition: TpfaCompressibleAssembler.hpp:232
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
void computePressuresAndFluxes(std::vector< double > &cell_pressures, std::vector< double > &face_pressures, std::vector< double > &face_fluxes, std::vector< double > &well_pressures, std::vector< double > &well_fluxes)
Compute cell pressures and face fluxes. You must call assemble() (and solve the linear system accesse...
Definition: TpfaCompressibleAssembler.hpp:258
~TpfaCompressibleAssembler()
Destructor.
Definition: TpfaCompressibleAssembler.hpp:54
FlowBCTypes
Boundary condition types.
Definition: TpfaCompressibleAssembler.hpp:136
@ 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
TpfaCompressibleAssembler()
Default constructor, does nothing.
Definition: TpfaCompressibleAssembler.hpp:43
void explicitTransport(const double dt, double *cell_surfvols)
Explicit IMPES transport.
Definition: TpfaCompressibleAssembler.hpp:329
void faceFluxToCellFlux(const std::vector< double > &face_fluxes, std::vector< double > &cell_fluxes)
Compute cell fluxes from face fluxes. You must call assemble() (and solve the linear system accessed ...
Definition: TpfaCompressibleAssembler.hpp:358
void assemble(const double *sources, const FlowBCTypes *bctypes, const double *bcvalues, const double dt, const double *totcompr, const double *voldiscr, const double *cellA, const double *faceA, const double *wellperfA, const double *phasemobf, const double *phasemobwellperf, const double *cell_pressure, const double *gravcapf, const double *wellperf_gpot, const double *surf_dens)
Assemble the sparse system. You must call init() prior to calling assemble().
Definition: TpfaCompressibleAssembler.hpp:151
const std::vector< double > & faceTransmissibilities() const
Definition: TpfaCompressibleAssembler.hpp:389
Encapsulate a sparse linear system in CSR format.
Definition: TpfaCompressibleAssembler.hpp:216
double * sa
Definition: TpfaCompressibleAssembler.hpp:221
double * b
Definition: TpfaCompressibleAssembler.hpp:222
int nnz
Definition: TpfaCompressibleAssembler.hpp:218
int * ja
Definition: TpfaCompressibleAssembler.hpp:220
int * ia
Definition: TpfaCompressibleAssembler.hpp:219
double * x
Definition: TpfaCompressibleAssembler.hpp:223
int n
Definition: TpfaCompressibleAssembler.hpp:217